diff --git a/Core/FEM/Optimization.py b/Core/FEM/Optimization.py
new file mode 100644
index 0000000..321bb43
--- /dev/null
+++ b/Core/FEM/Optimization.py
@@ -0,0 +1,94 @@
+import os
+import gc
+import numpy as np
+import pyvista as pv
+from datetime import datetime
+from scipy.optimize import differential_evolution, NonlinearConstraint, Bounds
+
+import FEM.Solver as solver
+import Meshing.modulation_envelope as mod_env
+
+
+class Optimization(solver.Solver):
+    def __init__(self,settings_file: dict, settings_header: str, electrode_system: str, units: str = 'mm') -> None:
+        super().__init__(settings_file, settings_header, electrode_system, units)
+        self.__region_volumes = {}
+        self.__AAL_regions = None
+        self.__settings_file = settings_file
+        self.__settings_header = settings_header
+
+    def initialization(self, model_name, max_solver_iterations=250, solver_relative_tol=1e-7, solver_absolute_tol=1e-3):
+        self.load_mesh(model_name)
+        self.solver_setup(max_solver_iterations, solver_relative_tol, solver_absolute_tol)
+
+        mesh = pv.UnstructuredGrid(self.__settings_file[self.__settings_header][model_name]['mesh_file'])
+        cell_volumes = mesh.compute_cell_sizes().cell_arrays['Volume']
+        self.__AAL_regions = mesh['AAL_regions']
+
+        for region in np.unique(self.__AAL_regions):
+            roi = np.where(self.__AAL_regions == region)[0]
+            self.__region_volumes[int(region)] = np.sum(cell_volumes[roi])
+        
+        del mesh
+
+    def objective_df(self, x, unwanted_regions=[], regions_of_interest=[]):
+        electrodes = np.round(x[:4]).astype(np.int32) # The first 4 indices are the electrode IDs
+        currents = np.round(x[4:], 3)
+        if (np.sum(currents) != 2) or (np.unique(np.round(x[:4]), return_counts=True)[1].size != 4):
+            return np.inf
+
+        self.essential_boundaries.clear()
+        self.fields.clear()
+        self.field_variables.clear()
+
+        self.define_field_variable('potential_base', 'voltage')
+        self.define_field_variable('potential_df', 'voltage')
+
+        electrode_names = [self._electrode_names[id] for id in electrodes]
+
+        self.define_essential_boundary(electrode_names[0], electrodes[0], 'potential_base', current=currents[0])
+        self.define_essential_boundary(electrode_names[1], electrodes[1], 'potential_base', current=-currents[0])
+        self.define_essential_boundary(electrode_names[2], electrodes[2], 'potential_df', current=currents[1])
+        self.define_essential_boundary(electrode_names[3], electrodes[3], 'potential_df', current=-currents[1])
+
+        solution = self.run_solver(save_results=False, post_process_calculation=True)
+
+        e_field_base = solution['e_field_(potential_base)'].data[:, 0, :, 0]
+        e_field_df = solution['e_field_(potential_df)'].data[:, 0, :, 0]
+        modulation_values = mod_env.modulation_envelope(e_field_base, e_field_df)
+
+        roi_region_sum = 0
+        non_roi_region_sum = 0
+
+        for region in np.unique(self.__AAL_regions):
+            if region in unwanted_regions:
+                continue
+            roi = np.where(self.__AAL_regions == region)[0]
+
+            if region in regions_of_interest:
+                roi_region_sum += np.sum(modulation_values[roi])/self.__region_volumes[region]
+            non_roi_region_sum += np.sum(modulation_values[roi])/self.__region_volumes[region]
+        
+        region_ratio = roi_region_sum/non_roi_region_sum
+        del solution
+        gc.collect()
+
+        return -np.round(region_ratio*10000, 1)
+
+    def run_optimization(self):
+        # TODO: Make variable count automatic
+        bounds = Bounds([10, 10, 10, 10, 0.05, 0.05], [70, 70, 70, 70, 2.0, 2.0])
+        unwanted_regions = np.append(np.arange(-90, -9), [-1, -2, -3, -7])
+        roi_ids = np.array([42])
+
+        # Constraints
+        unique_electrode_const = NonlinearConstraint(lambda x: np.unique(np.round(x[:4]), return_counts=True)[1].size, 4, 4) # Keep only unique combinations
+        current_sum_const = NonlinearConstraint(lambda x: np.round(x[4], 3) + np.round(x[5], 3), 2, 2)
+        current_step_const_1 = NonlinearConstraint(lambda x: (np.round(x[4], 3)/0.05).is_integer(), True, True)
+        current_step_const_2 = NonlinearConstraint(lambda x: (np.round(x[5], 3)/0.05).is_integer(), True, True)
+        consts = (unique_electrode_const, current_sum_const, current_step_const_1, current_step_const_2)
+
+        result = differential_evolution(self.objective_df, bounds, args=(unwanted_regions, roi_ids, ), constraints=consts, maxiter=100, disp=True)
+
+        return result
+
diff --git a/Core/FEM/Solver.py b/Core/FEM/Solver.py
new file mode 100644
index 0000000..dedb9b9
--- /dev/null
+++ b/Core/FEM/Solver.py
@@ -0,0 +1,266 @@
+from __future__ import absolute_import
+import os
+import gc
+import sys
+import numpy as np
+from numpy.core.defchararray import add
+import pyvista as pv
+import sfepy
+
+#### SfePy libraries
+from sfepy.base.base import Struct
+from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem, Function)
+from sfepy.discrete.fem import Mesh, FEDomain, Field
+from sfepy.terms import Term
+from sfepy.discrete.conditions import Conditions, EssentialBC
+from sfepy.solvers.ls import PETScKrylovSolver
+from sfepy.solvers.nls import Newton
+#### SfePy libraries
+
+class Solver:
+    def __init__(self, settings_file: dict, settings_header: str, electrode_system: str, units: str = 'mm'):
+        self.__settings: dict = settings_file
+        self.__settings_header: str = settings_header
+        if os.name == 'nt':
+            self.__extra_path = '_windows'
+        else:
+            self.__extra_path = ''
+        sys.path.append(os.path.realpath(settings_file[settings_header]['lib_path' + self.__extra_path]))
+
+        self.__linear_solver: sfepy.solvers.ls = None
+        self.__non_linear_solver: sfepy.solvers.nls = None
+        self.__overall_volume = None
+        self.__fields_to_calculate: list = []
+        self.__electrode_currents: dict = {}
+        self._electrode_names: dict = {}
+        self.conductivities: dict = {}
+        self.electrode_system: str = electrode_system
+        self.units: str = units
+        self.__out_of_range_assign_region: str = None
+        self.__out_of_range_group_threshold: int = None
+
+        # Read from settings
+        self.__material_conductivity = None
+        self.__selected_model: str = None
+        self.domain: sfepy.dicrete.fem.FEDomain = None
+        self.problem = None
+        self.essential_boundaries: list = []
+        self.field_variables: dict = {}
+        self.fields: dict = {}
+
+    def load_mesh(self, model: str=None, connectivity: str='3_4', id_array_name: str='cell_scalars') -> None:
+        if model is None:
+            raise AttributeError('No model was selected.')
+        mesh = pv.UnstructuredGrid(self.__settings[self.__settings_header][model]['mesh_file' + self.__extra_path])
+        self.__selected_model = model
+
+        vertices = mesh.points
+        vertex_groups = np.empty(vertices.shape[0])
+        cells = mesh.cell_connectivity.reshape((-1, 4)) # TODO: Generalize for higher order connectivity
+        cell_groups = mesh.cell_arrays[id_array_name]
+
+        for group in np.unique(cell_groups):
+            roi_cells = np.unique(cells[np.where(cell_groups == group)[0]])
+            vertex_groups[roi_cells] = group
+
+        loaded_mesh = Mesh.from_data('model_mesh', vertices, vertex_groups, [cells], [cell_groups], [connectivity])
+        self.domain = FEDomain('model_domain', loaded_mesh)
+
+        del mesh
+
+    def define_field_variable(self, field_var_name: str, field_name: str, out_of_range_assign_region: str = None, out_of_range_group_threshold: int = None) -> None:
+        self.__out_of_range_assign_region = out_of_range_assign_region
+        self.__out_of_range_group_threshold = out_of_range_group_threshold
+
+        if not self.__overall_volume:
+            self.__assign_regions(self.__out_of_range_assign_region, self.__out_of_range_group_threshold)
+        if field_name not in self.fields.keys():
+            self.fields[field_name] = Field.from_args(field_name, dtype=np.float64, shape=(1, ), region=self.__overall_volume, approx_order=1)
+            self.__electrode_currents.clear()
+
+        self.field_variables[field_var_name] = {
+            'unknown': FieldVariable(field_var_name, 'unknown', field=self.fields[field_name]),
+            'test': FieldVariable(field_var_name + '_test', 'test', field=self.fields[field_name], primary_var_name=field_var_name),
+        }
+
+    def define_essential_boundary(self, region_name: str, group_id: int, field_variable: str, potential: float = None, current: float = None) -> None:
+        assert field_variable in self.field_variables.keys(), 'The field variable {} is not defined'.format(field_variable)
+        assert (potential is not None) ^ (current is not None), 'Only potential or current value shall be provided.'
+
+        if current is not None:
+            try:
+                self.__electrode_currents[field_variable][region_name] = current
+            except KeyError:
+                self.__electrode_currents[field_variable] = {region_name: current}
+            potential = 1 if (current > 0) else -1
+
+        temporary_domain = self.domain.create_region(region_name, 'r.{} *v r.Skin'.format(region_name), 'facet', add_to_regions=False)
+        self.essential_boundaries.append(EssentialBC(region_name, temporary_domain, {field_variable + '.all' : potential}))
+
+    def solver_setup(self, max_iterations: int=250, relative_tol: float=1e-7, absolute_tol: float=1e-3, verbose: bool=False) -> None:
+        self.__linear_solver = PETScKrylovSolver({
+            'ksp_max_it': max_iterations,
+            'ksp_rtol': relative_tol,
+            'ksp_atol': absolute_tol,
+            'ksp_type': 'cg',
+            'pc_type': 'hypre',
+            'pc_hypre_type': 'boomeramg',
+            'pc_hypre_boomeramg_coarsen_type': 'HMIS',
+            'verbose': 2 if verbose else 0,
+        })
+
+        self.__non_linear_solver = Newton({
+            'i_max': 1,
+            'eps_a': absolute_tol,
+        }, lin_solver=self.__linear_solver)
+
+    def run_solver(self, save_results: bool, post_process_calculation: bool, field_calculation: list = ['E'], output_dir: str=None, output_file_name: str=None):
+        if not self.__non_linear_solver:
+            raise AttributeError('The solver is not setup. Please set it up before calling run.')
+        self.__material_definition()
+
+        if field_calculation:
+            self.__fields_to_calculate = field_calculation
+
+        self.problem = Problem('temporal_interference', equations=self.__generate_equations())
+        self.problem.set_bcs(ebcs=Conditions(self.essential_boundaries))
+        self.problem.set_solver(self.__non_linear_solver)
+        self.problem.setup_output(output_filename_trunk=output_file_name, output_dir=output_dir)
+
+        if post_process_calculation:
+            if save_results:
+                state = self.problem.solve(save_results=save_results, post_process_hook=self.__post_process)
+                return state
+            else:
+                state = self.problem.solve(save_results=save_results)
+            output_dictionary = state.create_output_dict()
+            output_dictionary = self.__post_process(output_dictionary, self.problem, state, extend=True)
+
+            del state
+            return output_dictionary
+        return self.problem.solve(save_results=save_results)
+
+    def set_custom_post_process(self, function) -> None:
+        self.__post_process = function
+
+    def clear_all(self) -> None:
+        del self.domain
+        del self.__overall_volume
+        del self.essential_boundaries
+        del self.field_variables
+        del self.fields
+        del self.problem
+        gc.collect()
+
+    def __generate_equations(self) -> sfepy.discrete.Equations:
+        # TODO: Add a check for the existence of the fields
+        integral = Integral('i1', order=2)
+
+        equations_list = []
+        for field_variable in self.field_variables.items():
+            term_arguments = {
+                'conductivity': self.__material_conductivity,
+                field_variable[0] + '_test': field_variable[1]['test'],
+                field_variable[0]: field_variable[1]['unknown']
+            }
+            equation_term = Term.new('dw_laplace(conductivity.val, ' + field_variable[0] + '_test, ' + field_variable[0] + ')', integral, self.__overall_volume, **term_arguments)
+            equations_list.append(Equation('equation_' + field_variable[0], equation_term))
+
+        return Equations(equations_list)
+
+    def __material_definition(self) -> None:
+        if not self.conductivities:
+            self.__assign_regions(self.__out_of_range_assign_region, self.__out_of_range_group_threshold)
+        self.__material_conductivity = Material('conductivity', function=Function('get_conductivity', lambda ts, coors, mode=None, equations=None, term=None, problem=None, **kwargs: self.__get_conductivity(ts, coors, mode, equations, term, problem, conductivities=self.conductivities)))
+
+    def __assign_regions(self, out_of_range_assign_region: str = None, out_of_range_group_threshold: int = None) -> None:
+        self.__overall_volume = self.domain.create_region('Omega', 'all')
+
+        if out_of_range_assign_region is not None:
+            unique_cell_groups = np.unique(self.domain.cmesh.cell_groups)
+            out_of_range_groups = unique_cell_groups[unique_cell_groups > out_of_range_group_threshold]
+
+            if out_of_range_groups.size > 0:
+                additions = ''
+                for group in out_of_range_groups:
+                    additions += 'cells of group {} +c '.format(group)
+                additions = ' +c c' + additions.strip('+c ')
+
+        for region in self.__settings[self.__settings_header][self.__selected_model]['regions'].items():
+            if (region[0] == out_of_range_assign_region) and (out_of_range_groups.size > 0):
+                self.domain.create_region(region[0], 'cells of group ' + str(region[1]['id']) + additions)
+            else:
+                self.domain.create_region(region[0], 'cells of group ' + str(region[1]['id']))
+            self.domain.create_region(region[0] + '_Gamma', 'vertices of group ' + str(region[1]['id']), 'facet')
+            self.conductivities[region[0]] = region[1]['conductivity']
+
+        for electrode in self.__settings[self.__settings_header]['electrodes'][self.electrode_system].items():
+            self.domain.create_region(electrode[0], 'cells of group ' + str(electrode[1]['id']))
+            self.domain.create_region(electrode[0] + '_Gamma', 'vertices of group ' + str(electrode[1]['id']), 'facet')
+            self.domain.create_region(electrode[0] + '_Gamma_cross', 'r.{} *v r.Skin'.format(electrode[0]), 'facet')
+            self.conductivities[electrode[0]] = self.__settings[self.__settings_header]['electrodes']['conductivity']
+            self._electrode_names[int(electrode[1]['id'])] = electrode[0]
+
+    def __get_conductivity(self, ts, coors, mode=None, equations=None, term=None, problem=None, conductivities=None):
+        # Execute only once at the initialization
+        if mode == 'qp':
+            values = np.empty(int(coors.shape[0]/4)) # Each element corresponds to one coordinate of the respective tetrahedral edge
+            num_nodes, dim = coors.shape
+            material_dict = {}
+
+            # Save the conductivity values
+            for domain in problem.domain.regions:
+                if domain.name in conductivities.keys():
+                    values[domain.entities[3]] = conductivities[domain.name]
+
+            # Values used in a matrix format for the material
+            tmp_fun = lambda x, dim: x*np.eye(dim) # Required for the diffusion velocity in the current density calculation
+
+            values = np.repeat(values, 4) # Account for the tetrahedral edges
+            if 'J' in np.char.upper(self.__fields_to_calculate):
+                mat_vec = np.array(list(map(tmp_fun, values, np.repeat(dim, num_nodes))))
+                material_dict['mat_vec'] = mat_vec
+
+            values.shape = (coors.shape[0], 1, 1)
+            material_dict['val'] = values
+
+            return material_dict
+
+    def __post_process(self, out, problem, state, extend=False):
+        electrode_conductivity = self.__settings[self.__settings_header]['electrodes']['conductivity']
+        electrode_material = Material('electrode', kind='stationary', values={'mat_vec': electrode_conductivity*np.eye(3)})
+        
+        # Save the output
+        for field_variable_name in self.field_variables.keys():
+            if self.units == 'mm':
+                distance_unit_multiplier = 1000
+            else:
+                distance_unit_multiplier = 1
+            
+            if self.__electrode_currents:
+                currents = list(self.__electrode_currents[field_variable_name].values())
+                regions = list(self.__electrode_currents[field_variable_name].keys())
+                assert np.sum(currents) == 0, 'The currents must sum to zero. The current sum is {}'.format(np.sum(currents))
+
+                surface_areas = []
+                fluxes = []
+                for region in regions:
+                    fluxes.append(1./(distance_unit_multiplier**2) * problem.evaluate('d_surface_flux.2.' + region + '_Gamma_cross(electrode.mat_vec, ' + field_variable_name + ')', mode='eval', copy_materials=False, electrode=electrode_material))
+                    surface_areas.append(problem.evaluate('d_surface.2.' + region + '_Gamma_cross(' + field_variable_name + ')', mode='eval'))
+                fluxes = np.abs(np.array(fluxes))
+                mean_current = np.average(fluxes*np.amax(surface_areas), weights=[surface_areas[0]/surface_areas[1], 1])
+                current_multiplier = np.abs(currents[0]*0.001)/mean_current # Current is always in mA
+
+                for potential in out.keys():
+                    out[potential].data *= current_multiplier
+
+            if 'E' in np.char.upper(self.__fields_to_calculate):
+                output_var_name = 'e_field_(' + field_variable_name + ')'
+                e_field = distance_unit_multiplier * current_multiplier * problem.evaluate('-ev_grad.2.Omega(' + field_variable_name + ')', mode='qp')
+                out[output_var_name] = Struct(name=output_var_name, mode='cell', data=e_field, dofs=None)
+            if 'J' in np.char.upper(self.__fields_to_calculate):
+                output_var_name = 'j_field_(' + field_variable_name + ')'
+                j_field = distance_unit_multiplier * current_multiplier * problem.evaluate('ev_diffusion_velocity.2.Omega(conductivity.mat_vec, ' + field_variable_name + ')', mode='qp', copy_materials=False)
+                out[output_var_name] = Struct(name=output_var_name, mode='cell', data=j_field, dofs=None)
+
+        return out
diff --git a/Core/FileOps/FileOperations.py b/Core/FileOps/FileOperations.py
new file mode 100644
index 0000000..596c9d5
--- /dev/null
+++ b/Core/FileOps/FileOperations.py
@@ -0,0 +1,141 @@
+import numpy as np
+
+class FileOperations():
+    def __init__(self):
+        print("here")
+
+    @staticmethod
+    def poly_write(file_name: str, nodes, faces, regions: dict, boundaries=None):
+        with open(file_name, "wb") as m_file:
+            # Nodes
+            m_file.write("{} 3\n".format(len(nodes)).encode("utf-8"))
+            for index in range(0, len(nodes)):
+                m_file.write("{} {} {} {}\n".format(index + 1, *nodes[index]).encode("utf-8"))
+
+            # Faces
+            if boundaries is not None:
+                m_file.write("{} 1\n".format(len(faces)).encode("utf-8"))
+                for index in range(0, len(faces)):
+                    m_file.write("1 0 {}\n".format(int(boundaries[index] + 1)).encode("utf-8"))
+                    m_file.write("3 {} {} {}\n".format(*faces[index] + 1).encode("utf-8"))
+            else:
+                m_file.write("{} 0\n".format(len(faces)).encode("utf-8"))
+                for index in range(0, len(faces)):
+                    m_file.write("1 0\n".encode("utf-8"))
+                    m_file.write("3 {} {} {}\n".format(*faces[index] + 1).encode("utf-8"))
+                
+
+            # Holes NOT IMPLEMENTED YET
+            m_file.write("0\n".encode("utf-8"))
+
+            # Regions
+            m_file.write("{}\n".format(len(list(regions.keys()))).encode("utf-8"))
+            for region_id in regions.keys():
+                m_file.write("{} {} {} {} {} {}\n".format(region_id, *regions[region_id]['coordinates'], region_id, regions[region_id]['max_volume']).encode("utf-8"))
+
+    @staticmethod
+    def gmsh_write(file_name: str, surfaces, domains, physical_tags, bounding_surface_tag):
+        c_size_t = np.dtype("P")
+
+        with open(file_name, "wb") as m_file:
+            # MeshFormat (Check)
+            m_file.write(b"$MeshFormat\n")
+            m_file.write("4.1 {} {}\n".format(1, c_size_t.itemsize).encode("utf-8"))
+
+            m_file.write(b"$EndMeshFormat\n")
+
+            # PhysicalNames (Check)
+            m_file.write(b"$PhysicalNames\n")
+            m_file.write("{}\n".format(len(physical_tags)).encode("utf-8"))
+            for physical_tag in physical_tags:
+                m_file.write('3 {} "{}"\n'.format(*physical_tag).encode("utf-8"))
+            m_file.write(b"$EndPhysicalNames\n")
+
+            # Entities (Check)
+            m_file.write(b"$Entities\n")
+            
+            m_file.write("0 0 {} {}\n".format(len(surfaces), len(domains)).encode("utf-8"))
+            for i in range(0, len(surfaces)):
+                m_file.write("{} {} {} {} {} {} {} {} {}\n".format(i + 1, np.amin(surfaces[i].vertices[:, 0]), np.amin(surfaces[i].vertices[:, 1]), np.amin(surfaces[i].vertices[:, 2]), np.amax(surfaces[i].vertices[:, 0]), np.amax(surfaces[i].vertices[:, 1]), np.amax(surfaces[i].vertices[:, 2]), 0, 0).encode("utf-8"))
+            for i in range(0, len(domains)):
+                m_file.write("{} {} {} {} {} {} {} {} {} {} {}\n".format(i + 1, np.amin(domains[i].vertices[:, 0]), np.amin(domains[i].vertices[:, 1]), np.amin(domains[i].vertices[:, 2]), np.amax(domains[i].vertices[:, 0]), np.amax(domains[i].vertices[:, 1]), np.amax(domains[i].vertices[:, 2]), 1, physical_tags[i][0], 1, bounding_surface_tag[i]).encode("utf-8"))
+            m_file.write(b"$EndEntities\n")
+
+            # Nodes (Checked?)
+            float_fmt=".16e"
+
+            surfs = 0
+            voxels = 0
+            for surf in surfaces:
+                surfs = surfs + surf.num_vertices
+            for domain in domains:
+                voxels = voxels + domain.num_vertices
+
+            m_file.write(b"$Nodes\n")
+            num_blocks = len(surfaces) + len(domains)
+            min_tag = 1
+            max_tag = surfs + voxels
+
+            m_file.write("{} {} {} {}\n".format(num_blocks, surfs + voxels, min_tag, max_tag).encode("utf-8"))
+
+            ta = 1
+            face_elements = []
+            for i in range(0, len(surfaces)):
+                # dim_entity, is_parametric
+                num_elem = surfaces[i].num_vertices
+
+                m_file.write("{} {} {} {}\n".format(2, i + 1, 0, num_elem).encode("utf-8"))
+                np.arange(ta, ta + num_elem).tofile(m_file, "\n", "%d")
+                m_file.write(b"\n")
+                np.savetxt(m_file, surfaces[i].vertices, delimiter=" ", fmt="%" + float_fmt)
+                #face_elements.append(surfaces[i].faces + 1 + (ta - 1))
+                face_elements.append(surfaces[i].faces.astype(c_size_t) + 1 + (ta - 1))
+                ta = ta + num_elem
+                m_file.write(b"\n")
+
+            voxel_elements = []
+            for i in range(0, len(domains)):
+                num_elem = domains[i].num_vertices
+
+                m_file.write("{} {} {} {}\n".format(3, i + 1, 0, num_elem).encode("utf-8"))
+                np.arange(ta, ta + num_elem).tofile(m_file, "\n", "%d")
+                m_file.write(b"\n")
+                np.savetxt(m_file, domains[i].vertices, delimiter=" ", fmt="%" + float_fmt)
+                #voxel_elements.append(domains[i].voxels + 1 + (ta - 1))
+                voxel_elements.append(domains[i].voxels.astype(c_size_t) + 1 + (ta - 1))
+                ta = ta + num_elem
+                m_file.write(b"\n")
+            m_file.write(b"$EndNodes\n")
+
+            # Elements
+            surfs = 0
+            voxels = 0
+            for surf in surfaces:
+                surfs = surfs + surf.num_faces
+            for domain in domains:
+                voxels = voxels + domain.num_voxels
+
+            m_file.write(b"$Elements\n")
+            num_blocks = len(face_elements) + len(voxel_elements)
+            min_tag = 1
+            max_tag = surfs + voxels
+
+            m_file.write("{} {} {} {}\n".format(num_blocks, surfs + voxels, min_tag, max_tag).encode("utf-8"))
+
+            ta = 1
+            for i in range(0, len(face_elements)):
+                num_elem = face_elements[i].shape[0]
+
+                m_file.write("{} {} {} {}\n".format(2, i + 1, 2, num_elem).encode("utf-8"))
+                np.savetxt(m_file, np.column_stack([np.arange(ta, ta + num_elem), face_elements[i]]), "%d", " ")
+                ta = ta + num_elem
+                m_file.write(b"\n")
+
+            for i in range(0, len(voxel_elements)):
+                num_elem = voxel_elements[i].shape[0]
+
+                m_file.write("{} {} {} {}\n".format(3, i + 1, 4, num_elem).encode("utf-8"))
+                np.savetxt(m_file, np.column_stack([np.arange(ta, ta + num_elem), voxel_elements[i]]), "%d", " ")
+                ta = ta + num_elem
+                m_file.write(b"\n")
+            m_file.write(b"$EndElements\n")
\ No newline at end of file
diff --git a/Core/Meshing/ElectrodeOperations.py b/Core/Meshing/ElectrodeOperations.py
new file mode 100644
index 0000000..77237f6
--- /dev/null
+++ b/Core/Meshing/ElectrodeOperations.py
@@ -0,0 +1,113 @@
+import pymesh
+import numpy as np
+
+class ElectrodeOperations:
+    def __init__(self, surface_mesh: pymesh.Mesh, electrode_attributes: dict):
+        self._surface_mesh: pymesh.Mesh = surface_mesh
+        self._surface_mesh.enable_connectivity() # Required by some functions below
+        self._electrode_attributes: dict = electrode_attributes
+        self.electrode_array: dict = {}
+
+    def add_electrodes_on_skin(self) -> list:
+        if not self.electrode_array:
+            raise AttributeError('The electrodes shall be positioned first before added on the surface. Please call the positioning function first.')
+        electrode_mesh = self.get_electrode_single_mesh()
+
+        # Get the surface outline including the electrode
+        model = pymesh.merge_meshes((self._surface_mesh, electrode_mesh))
+        outer_hull = pymesh.compute_outer_hull(model)
+
+        # Create the surface with the electrode mesh imprinted
+        electrode_tan_mesh = pymesh.boolean(electrode_mesh, self._surface_mesh, 'difference')
+        outer_diff = pymesh.boolean(outer_hull, electrode_tan_mesh, 'difference')
+        conditioned_surface = pymesh.merge_meshes((outer_diff, electrode_tan_mesh))
+
+        # Generate the surface with the electrode on
+        face_id = np.arange(conditioned_surface.num_faces)
+        conditioned_surface = pymesh.remove_duplicated_vertices(conditioned_surface)[0]
+        non_intersecting_surface = np.isin(face_id, pymesh.detect_self_intersection(conditioned_surface)[:, 0], invert=True)
+
+        return [pymesh.submesh(conditioned_surface, non_intersecting_surface, 0), outer_diff]
+
+    def standard_electrode_positioning(self, electrodes_to_omit: list=None) -> None:
+        width = self._electrode_attributes['width']
+        radius = self._electrode_attributes['radius']
+        elements = self._electrode_attributes['elements']
+
+        closest_points = pymesh.distance_to_mesh(self._surface_mesh, self._electrode_attributes['coordinates'])[1] # Get the closest point to the one provided
+
+        for electrode_name, closest_point in zip(self._electrode_attributes['names'], closest_points):
+            if electrode_name in electrodes_to_omit:
+                continue
+            p_i = self._surface_mesh.vertices[self._surface_mesh.faces[closest_point]][0] # Get the surface vertex coordinates
+            electrode_orientation = self.__orient_electrode(p_i) # Orient the electrode perpendicular to the surface
+
+            # Generate the electrode cylinder and save to the output dictionary
+            electrode_cylinder = pymesh.generate_cylinder(p_i - (width * electrode_orientation)/4., p_i + (width * electrode_orientation)/4., radius, radius, elements)
+            self.electrode_array[electrode_name] = electrode_cylinder
+
+    def sphere_electrode_positioning(self) -> None:
+        cyl_radius = self._electrode_attributes['cylinder_radius']
+        #skin_radius = self._electrode_attributes['skin_radius']
+        skin_radius = np.amax(self._surface_mesh.vertices[:, 0])
+        cyl_width = self._electrode_attributes['cylinder_width']
+        elements = self._electrode_attributes['elements']
+
+        for electrode in self._electrode_attributes['electrodes'].items():
+            print(electrode)
+            p_i = self.electrode_position_sphere(skin_radius, electrode[1]['theta'], electrode[1]['phi'])
+            electrode_orientation = self.__orient_electrode_sphere(p_i, np.array([0, 0, -cyl_width]))
+
+            # Generate the electrode cylinder and save to the output dictionary
+            electrode_cylinder = pymesh.generate_cylinder(p_i - (cyl_width * electrode_orientation)/4., p_i + (cyl_width * electrode_orientation)/4., cyl_radius, cyl_radius, elements)
+            self.electrode_array[electrode[0]] = electrode_cylinder
+
+    def get_electrode_array(self) -> dict:
+        if not self.electrode_array:
+            raise AttributeError('Electrodes are not positioned yet. Please call the positioning function.')
+        return self.electrode_array
+
+    def get_electrode_single_mesh(self) -> pymesh.Mesh:
+        assert self.electrode_array, 'Electrodes are not positioned yet. Please call the positioning function.'
+
+        merged_meshes = pymesh.merge_meshes((e_mesh for e_mesh in self.electrode_array.values()))
+        merged_ids = merged_meshes.get_attribute('face_sources')
+        new_attribute_ids = np.empty(merged_meshes.num_faces, dtype=int)
+
+        for electrode_face_source, electrode_name in zip(np.unique(merged_ids), self.electrode_array.keys()):
+            id_array_positions = np.where(merged_ids == electrode_face_source)[0]
+            new_attribute_ids[id_array_positions] = self._electrode_attributes['ids'][electrode_name]['id']
+        merged_meshes.set_attribute('face_sources', new_attribute_ids)
+        return merged_meshes
+
+    def __orient_electrode(self, init_point: list):
+        point_id = np.where(np.sum(self._surface_mesh.vertices == init_point, axis=1))[0][0] # Unique point assumed
+        face = self._surface_mesh.get_vertex_adjacent_faces(point_id)[0]
+
+        points = []
+        for point in self._surface_mesh.vertices[self._surface_mesh.faces[face]]:
+            if np.sum(point != init_point):
+                points.append(point)
+        p_1 = points[0] - init_point
+        p_2 = points[1] - init_point
+
+        normal = np.cross(p_1, p_2)
+        return normal/np.linalg.norm(normal)
+
+    def __orient_electrode_sphere(self, init_point: list, delta_point: list):
+        dist = pymesh.signed_distance_to_mesh(self._surface_mesh, init_point)
+        face = dist[1][0]
+
+        # Create a vector with the direction of the face
+        p_1 = self._surface_mesh.vertices[self._surface_mesh.faces[face][0]]
+        p_2 = self._surface_mesh.vertices[self._surface_mesh.faces[face][1]]
+        dir_vector = p_1 - p_2
+        dir_vector = dir_vector/np.linalg.norm(dir_vector)
+
+        normal = np.cross(delta_point, dir_vector)
+        return normal/np.linalg.norm(normal)
+
+    @staticmethod
+    def electrode_position_sphere(radius: float, theta: float, phi: float=0):
+        return np.array([radius*np.cos(np.deg2rad(phi))*np.cos(np.deg2rad(theta)), radius*np.cos(np.deg2rad(phi))*np.sin(np.deg2rad(theta)), radius*np.sin(np.deg2rad(phi))])
+
diff --git a/Core/Meshing/MeshOperations.py b/Core/Meshing/MeshOperations.py
new file mode 100644
index 0000000..d65b8c6
--- /dev/null
+++ b/Core/Meshing/MeshOperations.py
@@ -0,0 +1,124 @@
+import os
+import pymesh
+import numpy as np
+import itertools
+import warnings
+
+import Meshing.ElectrodeOperations as ElecOps
+import FileOps.FileOperations as FileOps
+
+class MeshOperations(ElecOps.ElectrodeOperations, FileOps.FileOperations):
+    def __init__(self, surface_mesh: pymesh.Mesh, electrode_attributes: dict):
+        self.merged_meshes: pymesh.Mesh = None
+        self.skin_with_electrodes: list = None
+        self.electrode_mesh: pymesh.Mesh = None
+        self.surface_meshes: list = []
+        ElecOps.ElectrodeOperations.__init__(self, surface_mesh, electrode_attributes)
+
+    def phm_model_meshing(self, mesh_filename: str, resolve_intersections: bool=True, electrodes_to_omit: list=None) -> None:
+        """Generate a .poly file of the provided model mesh.
+
+        Args:
+            mesh_filename (str): The path where to save the meshed model file
+            resolve_intersections (bool, optional): Whether to resolve any self-intersections the model might have. Defaults to True.
+
+        Raises:
+            AttributeError: If the class variable `surface_meshes` does not contain the model meshes, throw the error.
+        """
+        if not self.surface_meshes:
+            raise AttributeError('Meshes must be loaded first. Please load the meshes.')
+        # self.electrode_meshing(electrodes_to_omit=electrodes_to_omit)
+
+        # self.merged_meshes = pymesh.merge_meshes((self.skin_with_electrodes, *self.surface_meshes[1:]))
+        self.merged_meshes = pymesh.merge_meshes((*self.surface_meshes,))
+        regions = self.region_points(self.surface_meshes, 0.1, electrode_mesh=None)
+
+        self_intersection = pymesh.detect_self_intersection(self.merged_meshes)
+        if self_intersection.size != 0:
+            if resolve_intersections:
+                warnings.warn('Self intersections were detected. Trying to resolve', UserWarning)
+            else:
+                raise Warning('Self intersections were detected. User selected not to resolve.')
+
+            temp_meshes = []
+            combination_pairs = 2
+            surface_meshes = np.array(self.surface_meshes)
+            # surface_meshes[0] = self.skin_with_electrodes
+
+            unique_intersections = np.unique(self.merged_meshes.get_attribute('face_sources')[self_intersection]).astype(np.int32)
+            intersection_combinations = np.array(list(itertools.combinations(unique_intersections, combination_pairs))).astype(np.int32)
+
+            for combination in intersection_combinations:
+                if pymesh.detect_self_intersection(pymesh.merge_meshes((*surface_meshes[combination], ))).size != 0:
+                    temp_meshes.append(pymesh.resolve_self_intersection(pymesh.merge_meshes((*surface_meshes[combination], ))))
+
+            mask_intersecting = np.ones(len(surface_meshes), np.bool)
+            mask_intersecting[unique_intersections] = False
+
+            temp_meshes = temp_meshes + [surface_meshes[j] for j, i in enumerate(mask_intersecting) if mask_intersecting[j]]
+            temp_merged_meshes = pymesh.merge_meshes((*temp_meshes, ))
+
+            self.poly_write(mesh_filename, temp_merged_meshes.vertices, temp_merged_meshes.faces, regions)
+        else:
+            self.poly_write(mesh_filename, self.merged_meshes.vertices, self.merged_meshes.faces, regions)
+
+    def sphere_model_meshing(self, mesh_filename: str) -> None:
+        if not self.surface_meshes:
+            raise AttributeError('Meshes must be loaded first. Please load the meshes.')
+        self.electrode_meshing(sphere=True)
+
+        self.merged_meshes = pymesh.merge_meshes((self.skin_with_electrodes, *self.surface_meshes[1:]))
+        regions = self.region_points(self.surface_meshes, 0.1, electrode_mesh=self.electrode_mesh)
+
+        self.poly_write(mesh_filename, self.merged_meshes.vertices, self.merged_meshes.faces, regions)
+
+    def load_surface_meshes(self, base_path: str, file_names: list) -> None:
+        # TODO: Order matters
+        for file_name in file_names:
+            self.surface_meshes.append(pymesh.load_mesh(os.path.join(base_path, file_name)))
+
+    def electrode_meshing(self, sphere: bool=False, electrodes_to_omit: list=None):
+        print("Placing Electrodes") # INFO log
+        if sphere:
+            self.sphere_electrode_positioning()
+        else:
+            self.standard_electrode_positioning(electrodes_to_omit)
+        self.skin_with_electrodes = self.add_electrodes_on_skin()[0]
+
+        electrodes_single_mesh = self.get_electrode_single_mesh()
+        self.electrode_mesh = pymesh.boolean(electrodes_single_mesh, self._surface_mesh, 'difference')
+        pymesh.map_face_attribute(electrodes_single_mesh, self.electrode_mesh, 'face_sources')
+
+        return self.electrode_mesh
+
+    @staticmethod
+    def region_points(boundary_surfaces: list, shift: float, electrode_mesh: pymesh.Mesh=None, max_volume=30):
+        points = {}
+        for index, surface in enumerate(boundary_surfaces):
+            vertices = surface.vertices
+            max_z_point_id = np.where(vertices[:, 2] == np.amax(vertices[:, 2]))[0][0]
+            max_z_point_unit_vector = vertices[max_z_point_id]/np.linalg.norm(vertices[max_z_point_id])
+
+            points[str(index + 1)] = {
+                'coordinates': np.round(vertices[max_z_point_id] - np.absolute(np.multiply(max_z_point_unit_vector, [0, 0, shift])), 6),
+                'max_volume': max_volume,
+            }
+
+        if electrode_mesh is not None:
+            electrode_regions = electrode_mesh.get_attribute('face_sources')
+            electrode_mesh.add_attribute('vertex_valance')
+            vertex_valance = electrode_mesh.get_attribute('vertex_valance')
+
+            for region in np.unique(electrode_regions):
+                faces = electrode_mesh.faces[np.where(electrode_regions == region)[0]]
+                vertices = electrode_mesh.vertices[np.unique(faces)]
+
+                max_valance_point = np.where(vertex_valance[np.unique(faces)] == np.amax(vertex_valance[np.unique(faces)]))[0][0]
+                max_valance_unit_vector = vertices[max_valance_point]/np.linalg.norm(vertices[max_valance_point])
+
+                points[str(region)] = {
+                    'coordinates': np.round(vertices[max_valance_point] - np.multiply(max_valance_unit_vector, shift), 6),
+                    'max_volume': max_volume,
+                }
+
+        return points
diff --git a/Core/NiiPHM/NiiMesh.py b/Core/NiiPHM/NiiMesh.py
new file mode 100644
index 0000000..d4230ce
--- /dev/null
+++ b/Core/NiiPHM/NiiMesh.py
@@ -0,0 +1,102 @@
+from typing import Union, Tuple
+import os
+import PVGeo
+import numpy as np
+import pandas as pd
+import pyvista as pv
+import nibabel as nib
+
+class NiiMesh:
+    def __init__(self, mesh_path: str=None, index_array_name: str='cell_scalars', roi_ids: list=[1, 8]) -> None:
+        self._mesh_pts: pv.UnstructuredGrid = None
+        self.__index_array_name: str = index_array_name
+        self.__intial_bounds: dict = {}
+        self.voxels: pv.UniformGrid = None
+        self.affine: np.array = None
+        self.intensities: np.array = None
+        self._mesh: pv.UnstructuredGrid = self.load_mesh(os.path.realpath(mesh_path), index_array_name, roi_ids) if mesh_path else None
+
+    def load_mesh(self, mesh_path: str, index_array_name: str='cell_scalars', roi_ids: list=[1, 8]) -> None:
+        self.__index_array_name = index_array_name
+        self._mesh = pv.UnstructuredGrid(os.path.realpath(mesh_path))
+        self._mesh_pts = self._mesh.threshold(value=roi_ids, scalars=self.__index_array_name, preference='cell').cell_centers()
+
+        return self._mesh
+
+    def assign_intensities(self, assign_values_per_region: bool, assign_index: bool=False, values_to_assign=None, unwanted_region_ids: list=None):
+        roi_pts_loc = np.isin(self._mesh_pts[self.__index_array_name], unwanted_region_ids, invert=True)
+        roi_pts_ids = self._mesh_pts[self.__index_array_name][roi_pts_loc]
+        roi_intensities = np.zeros(roi_pts_ids.shape)
+
+        self.intensities = np.zeros(self._mesh_pts[self.__index_array_name].shape)
+        assert (assign_index == False) and (values_to_assign is not None), 'To assign values a list of those shall be given.'
+
+        if assign_values_per_region == False:
+            assert len(roi_intensities) == len(values_to_assign), 'The provided intensity values shall be the same length as the points.'
+
+        for index, id in enumerate(np.unique(roi_pts_ids)):
+            pts_loc = np.where(roi_pts_ids == id)[0]
+            if assign_index:
+                roi_intensities[pts_loc] = 0 if (id in unwanted_region_ids) else id
+            else:
+                if assign_values_per_region:
+                    roi_intensities[pts_loc] = 0 if (id in unwanted_region_ids) else values_to_assign[index]
+                else:
+                    roi_intensities[pts_loc] = 0 if (id in unwanted_region_ids) else values_to_assign[pts_loc]
+        self.intensities[roi_pts_loc] = roi_intensities
+
+        return self.intensities
+
+    def generate_uniform_grid(self, distance_margin: float, voxel_size: float, interp_radius_multiplier: int=4, interp_sharpness: int=8, **kwargs) -> Tuple[pv.UniformGrid, np.array]:
+        if self.intensities is None:
+            _ = self.assign_intensities(kwargs)
+        data = {'x': self._mesh_pts.points[:, 0], 'y': self._mesh_pts.points[:, 1], 'z': self._mesh_pts.points[:, 2], 'intensities': self.intensities}
+        df = pd.DataFrame(data, columns=['x', 'y', 'z', 'intensities'])
+        grid_points = PVGeo.points_to_poly_data(df)
+
+        bounds = self._mesh.bounds
+        x_bounds_init = np.ceil(np.sum(np.abs(bounds[0:2]) + distance_margin)/voxel_size).astype(np.int32)
+        y_bounds_init = np.ceil(np.sum(np.abs(bounds[2:4]) + distance_margin)/voxel_size).astype(np.int32)
+        z_bounds_init = np.ceil(np.sum(np.abs(bounds[4:6]) + distance_margin)/voxel_size).astype(np.int32)
+        self.__intial_bounds = {
+                'x_bounds': x_bounds_init, 
+                'y_bounds': y_bounds_init, 
+                'z_bounds': z_bounds_init
+            }
+
+        grid = pv.UniformGrid((x_bounds_init, y_bounds_init, z_bounds_init))
+        grid.origin = -1*np.array([np.sum(np.abs(bounds[0:2]))/2., np.sum(np.abs(bounds[2:4]))/2., np.sum(np.abs(bounds[4:6]))/2.]) - distance_margin
+        grid.spacing = [voxel_size]*3
+
+        self.voxels = grid.interpolate(grid_points, radius=voxel_size*interp_radius_multiplier, sharpness=interp_sharpness)
+
+        ## Affine calculation
+        self.affine = self.__calculate_affine(bounds, distance_margin, voxel_size)
+
+        return self.voxels, self.affine
+
+    def __calculate_affine(self, bounds: list, distance_margin: float, voxel_size: float):
+        # Get the lower end bounds. The bounds list contains elements as [x-, x+, y-, y+, z-, z+]
+        x_minus_bound = bounds[0] - distance_margin
+        y_minus_bound = bounds[2] - distance_margin
+        z_minus_bound = bounds[4] - distance_margin
+        affine = np.array([[voxel_size, 0, 0, x_minus_bound], [0, voxel_size, 0, y_minus_bound], [0, 0, voxel_size, z_minus_bound], [0, 0, 0, 1]])
+
+        return affine
+
+    def generate_nifti(self, image_units: tuple=('mm', 'sec'), save_image: bool=True, image_path: str=None, **knwargs):
+        if self.voxels is None:
+            _ = self.generate_uniform_grid(knwargs)
+        img_header = nib.Nifti1Header()
+        img_header.set_xyzt_units(*image_units)
+        voxel_data = self.voxels['intensities'].reshape((self.__intial_bounds['z_bounds'], self.__intial_bounds['y_bounds'], self.__intial_bounds['x_bounds'])).transpose()
+
+        img = nib.Nifti1Image(voxel_data, self.affine, img_header)
+
+        if save_image:
+            self.save_nifti('converted.nii' if image_path is None else image_path, img)
+
+        return img
+
+    def save_nifti(self, image_path: str, image_object: Union[nib.Nifti1Image, nib.Nifti2Image]):
+        nib.save(image_object, os.path.realpath(image_path))
diff --git a/Core/modulation_envelope.py b/Core/modulation_envelope.py
new file mode 100644
index 0000000..e75d232
--- /dev/null
+++ b/Core/modulation_envelope.py
@@ -0,0 +1,133 @@
+import itertools
+import numpy as np
+
+try:
+    import cupy as cp
+
+    def modulation_envelope_gpu(e_field_1, e_field_2, dir_vector=None):
+        if dir_vector is None:
+            envelope = cp.zeros(e_field_1.shape[0])
+
+            # Calculate the angles between the two fields for each vector
+            dot_angle = cp.einsum('ij,ij->i', e_field_1, e_field_2)
+            cross_angle = cp.linalg.norm(cp.cross(e_field_1, e_field_2), axis=1)
+            angles = cp.arctan2(cross_angle, dot_angle)
+
+            # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
+            e_field_2 = cp.where(cp.broadcast_to(angles >= cp.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
+
+            # Recalculate the angles
+            dot_angle = cp.einsum('ij,ij->i', e_field_1, e_field_2)
+            cross_angle = cp.linalg.norm(cp.cross(e_field_1, e_field_2), axis=1)
+            angles = cp.arctan2(cross_angle, dot_angle)
+            E_minus = cp.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+
+            # Condition to have two times the E2 field amplitude
+            max_condition_1 = cp.linalg.norm(e_field_2, axis=1) < cp.linalg.norm(e_field_1, axis=1)*cp.cos(angles)
+            e1_gr_e2 = cp.where(cp.linalg.norm(e_field_1, axis=1) > cp.linalg.norm(e_field_2, axis=1), max_condition_1, False)
+
+            # Condition to have two times the E1 field amplitude
+            max_condition_2 = cp.linalg.norm(e_field_1, axis=1) < cp.linalg.norm(e_field_2, axis=1)*cp.cos(angles)
+            e2_gr_e1 = cp.where(cp.linalg.norm(e_field_2, axis=1) > cp.linalg.norm(e_field_1, axis=1), max_condition_2, False)
+
+            # Double magnitudes
+            envelope = cp.where(e1_gr_e2, 2.0*cp.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
+            envelope = cp.where(e2_gr_e1, 2.0*cp.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
+
+            # Calculate the complement area to the previous calculation
+            e1_gr_e2 = cp.where(cp.linalg.norm(e_field_1, axis=1) > cp.linalg.norm(e_field_2, axis=1), cp.logical_not(max_condition_1), False)
+            e2_gr_e1 = cp.where(cp.linalg.norm(e_field_2, axis=1) > cp.linalg.norm(e_field_1, axis=1), cp.logical_not(max_condition_2), False)
+
+            # Cross product
+            envelope = cp.where(e1_gr_e2, 2.0*(cp.linalg.norm(cp.cross(e_field_2, E_minus), axis=1)/cp.linalg.norm(E_minus, axis=1)), envelope) # (First case)
+            envelope = cp.where(e2_gr_e1, 2.0*(cp.linalg.norm(cp.cross(e_field_1, -E_minus), axis=1)/cp.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
+        else:
+            # Calculate the values of the E field modulation envelope along the desired n direction
+            E_plus = cp.add(e_field_1, e_field_2) # Create the sum of the E fields
+            E_minus = cp.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+            envelope = cp.abs(cp.abs(cp.dot(E_plus, dir_vector)) - cp.abs(cp.dot(E_minus, dir_vector)))
+
+        return cp.nan_to_num(envelope)
+
+    def modulation_envelope_gpu_multi(*args):
+        assert len(args)%2 == 0, 'The electric fields shall be in pairs.'
+        unique_combinations = np.array(list(itertools.combinations(np.arange(0, len(args)), 2))).astype(np.int32)
+        unique_combinations_gpu = cp.array(unique_combinations)
+        e_fields = []
+        
+        for index in unique_combinations_gpu:
+            e_fields.append(modulation_envelope_gpu(args[int(index[0])], args[int(index[1])]))
+        
+        e_fields_gpu = cp.array(e_fields)
+        return cp.array(list(map(cp.amin, e_fields_gpu.transpose())))
+except ImportError:
+    pass
+
+def modulation_envelope(e_field_1, e_field_2, dir_vector=None):
+    """[summary]
+
+    Args:
+        e_field_1 ([type]): [description]
+        e_field_2 ([type]): [description]
+        dir_vector ([type], optional): [description]. Defaults to None.
+
+    Returns:
+        [type]: [description]
+    """
+    if dir_vector is None:
+        envelope = np.zeros(e_field_1.shape[0])
+
+        # Calculate the angles between the two fields for each vector
+        dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
+        cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
+        angles = np.arctan2(cross_angle, dot_angle)
+
+        # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
+        e_field_2 = np.where(np.broadcast_to(angles >= np.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
+
+        # Recalculate the angles
+        dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
+        cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
+        angles = np.arctan2(cross_angle, dot_angle)
+
+        E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+
+        # Condition to have two times the E2 field amplitude
+        max_condition_1 = np.linalg.norm(e_field_2, axis=1) < np.linalg.norm(e_field_1, axis=1)*np.cos(angles)
+        e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), max_condition_1, False)
+
+        # Condition to have two times the E1 field amplitude
+        max_condition_2 = np.linalg.norm(e_field_1, axis=1) < np.linalg.norm(e_field_2, axis=1)*np.cos(angles)
+        e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), max_condition_2, False)
+
+        # Double magnitudes
+        envelope = np.where(e1_gr_e2, 2.0*np.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
+        envelope = np.where(e2_gr_e1, 2.0*np.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
+
+        # Calculate the complement area to the previous calculation
+        e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), np.logical_not(max_condition_1), False)
+        e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), np.logical_not(max_condition_2), False)
+
+        # Cross product
+        envelope = np.where(e1_gr_e2, 2.0*(np.linalg.norm(np.cross(e_field_2, E_minus), axis=1)/np.linalg.norm(E_minus, axis=1)), envelope) # (First case)
+        envelope = np.where(e2_gr_e1, 2.0*(np.linalg.norm(np.cross(e_field_1, -E_minus), axis=1)/np.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
+    else:
+        # Calculate the values of the E field modulation envelope along the desired n direction
+        E_plus = np.add(e_field_1, e_field_2) # Create the sum of the E fields
+        E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+        envelope = np.abs(np.abs(np.dot(E_plus, dir_vector)) - np.abs(np.dot(E_minus, dir_vector)))
+    return np.nan_to_num(envelope)
+
+def modulation_envelope_multi(*args):
+    assert len(args)%2 == 0, 'The electric fields shall be in pairs.'
+    unique_combinations = np.array(list(itertools.combinations(np.arange(0, len(args)), 2))).astype(np.int32)
+    e_fields = []
+    
+    for index in unique_combinations:
+        e_fields.append(modulation_envelope(args[int(index[0])], args[int(index[1])]))
+    
+    e_fields = np.array(e_fields)
+    return np.array(list(map(np.amin, e_fields.transpose())))
+
+def mod_env_multi_multiprocessing(x):
+    return modulation_envelope(x[0], x[1])
diff --git a/Scripts/Data Analysis/[Optimization]_Result_Validation.py b/Scripts/Data Analysis/[Optimization]_Result_Validation.py
new file mode 100644
index 0000000..f3d3712
--- /dev/null
+++ b/Scripts/Data Analysis/[Optimization]_Result_Validation.py	
@@ -0,0 +1,137 @@
+import os
+import gc
+import numpy as np
+import cupy as cp
+import pandas as pd
+
+def modulation_envelope_gpu(e_field_1, e_field_2, dir_vector=None):
+    if dir_vector is None:
+        envelope = cp.zeros(e_field_1.shape[0])
+        
+        # Calculate the angles between the two fields for each vector
+        dot_angle = cp.einsum('ij,ij->i', e_field_1, e_field_2)
+        cross_angle = cp.linalg.norm(cp.cross(e_field_1, e_field_2), axis=1)
+        angles = cp.arctan2(cross_angle, dot_angle)
+        
+        # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
+        e_field_2 = cp.where(cp.broadcast_to(angles >= cp.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
+        
+        # Recalculate the angles
+        dot_angle = cp.einsum('ij,ij->i', e_field_1, e_field_2)
+        cross_angle = cp.linalg.norm(cp.cross(e_field_1, e_field_2), axis=1)
+        angles = cp.arctan2(cross_angle, dot_angle)
+        E_minus = cp.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+        
+        # Condition to have two times the E2 field amplitude
+        max_condition_1 = cp.linalg.norm(e_field_2, axis=1) < cp.linalg.norm(e_field_1, axis=1)*cp.cos(angles)
+        e1_gr_e2 = cp.where(cp.linalg.norm(e_field_1, axis=1) > cp.linalg.norm(e_field_2, axis=1), max_condition_1, False)
+        
+        # Condition to have two times the E1 field amplitude
+        max_condition_2 = cp.linalg.norm(e_field_1, axis=1) < cp.linalg.norm(e_field_2, axis=1)*cp.cos(angles)
+        e2_gr_e1 = cp.where(cp.linalg.norm(e_field_2, axis=1) > cp.linalg.norm(e_field_1, axis=1), max_condition_2, False)
+        
+        # Double magnitudes
+        envelope = cp.where(e1_gr_e2, 2.0*cp.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
+        envelope = cp.where(e2_gr_e1, 2.0*cp.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
+        
+        # Calculate the complement area to the previous calculation
+        e1_gr_e2 = cp.where(cp.linalg.norm(e_field_1, axis=1) > cp.linalg.norm(e_field_2, axis=1), cp.logical_not(max_condition_1), False)
+        e2_gr_e1 = cp.where(cp.linalg.norm(e_field_2, axis=1) > cp.linalg.norm(e_field_1, axis=1), cp.logical_not(max_condition_2), False)
+        
+        # Cross product
+        envelope = cp.where(e1_gr_e2, 2.0*(cp.linalg.norm(cp.cross(e_field_2, E_minus), axis=1)/cp.linalg.norm(E_minus, axis=1)), envelope) # (First case)
+        envelope = cp.where(e2_gr_e1, 2.0*(cp.linalg.norm(cp.cross(e_field_1, -E_minus), axis=1)/cp.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
+    else:
+        # Calculate the values of the E field modulation envelope along the desired n direction
+        E_plus = cp.add(e_field_1, e_field_2) # Create the sum of the E fields
+        E_minus = cp.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+        envelope = cp.abs(cp.abs(cp.dot(E_plus, dir_vector)) - cp.abs(cp.dot(E_minus, dir_vector)))
+    
+    return cp.nan_to_num(envelope)
+
+
+def objective_df(x, field_data, regions_of_interest, aal_regions, region_volumes, currents, threshold):  
+    if np.unique(np.round(x[:4]), return_counts=True)[1].size != 4:
+        return 100*(np.abs((np.unique(np.round(x[:4]), return_counts=True)[1].size - 4))**2) + 10000
+    
+    penalty = 0
+    electrodes = np.round(x[:4]).astype(np.int32) # The first 4 indices are the electrode IDs
+    roi = cp.isin(aal_regions, cp.array(regions_of_interest))
+    max_vals = []
+    fitness_vals = []
+    
+    for current in currents:
+        e_field_base = current[0]*field_data[electrodes[0]] - current[0]*field_data[electrodes[1]]
+        e_field_df = current[1]*field_data[electrodes[2]] - current[1]*field_data[electrodes[3]]
+
+        e_field_base_gpu = cp.array(e_field_base)
+        e_field_df_gpu = cp.array(e_field_df)
+        modulation_values = modulation_envelope_gpu(e_field_base_gpu, e_field_df_gpu)
+        max_vals.append(float(cp.amax(modulation_values[roi])))
+    
+        roi_region_sum = 0
+        non_roi_region_sum = 0
+
+        for region in cp.unique(aal_regions):
+            roi = cp.where(aal_regions == region)[0]
+        
+            if int(region) in regions_of_interest:
+                roi_region_sum += cp.sum(modulation_values[roi])/region_volumes[int(region)]
+            else:
+                non_roi_region_sum += cp.sum(modulation_values[roi])/region_volumes[int(region)]
+        
+        region_ratio = cp.nan_to_num(roi_region_sum/non_roi_region_sum)
+        fitness_measure = region_ratio*10000
+        fitness_vals.append(float(fitness_measure))
+    
+    max_vals = np.array(max_vals)
+    max_val_curr = np.where(max_vals >= threshold)[0]
+    fitness_vals = np.array(fitness_vals)
+
+    return_fitness = 0
+    if max_val_curr.size == 0:
+        penalty += 100*((threshold - np.mean(max_vals))**2) + 1000
+        return_fitness = np.amin(fitness_vals)
+    else:
+        fitness_candidate = np.amax(fitness_vals[max_val_curr])
+        return_fitness = fitness_candidate
+    
+    return -float(np.round(return_fitness - penalty, 2))
+
+if __name__ == "__main__":
+    npz_dir = 'directory of the NPZ model files'
+    csv_dir = 'directory of the optimized electrode pairs'
+    files = next(os.walk(csv_dir))[2]
+    
+    for model in files:
+        model = model.split('.')[0].split('_')[-1]
+        
+        npz_arrays = np.load(os.path.join(npz_dir, model + '_fields_brain_reg.npz'), allow_pickle=True)
+        electrode_vals = pd.read_csv(os.path.join(csv_dir, 'optimized_electrodes_' + model + '.csv'))
+        
+        field_data = npz_arrays['e_field']
+        aal_regions = npz_arrays['aal_ids']
+        region_volumes = {}
+        for region in np.unique(aal_regions):
+            region_volumes[region] = np.where(aal_regions == region)[0].size
+        roi_ids = np.array([42])
+        
+        ideal_case = None
+        
+        cur_potential_values = np.arange(0.5, 1.55, 0.05)
+        cur_x, cur_y = np.meshgrid(cur_potential_values, cur_potential_values)
+        cur_all_combinations = np.hstack((cur_x.reshape((-1, 1)), cur_y.reshape((-1, 1))))
+        usable_currents = cur_all_combinations[np.where(np.sum(np.round(cur_all_combinations, 2), axis=1) == 2)[0]]
+        
+        print(model)
+        electrodes_model = electrode_vals['electrodes'].to_numpy().astype(int)
+        aal_regions_gpu = cp.array(aal_regions)
+        value, cur_id = objective_df(electrodes_model, field_data, roi_ids, aal_regions=aal_regions_gpu, region_volumes=region_volumes, currents=usable_currents, ideal_case=ideal_case)
+        
+        electrode_vals['currents'] = np.round([usable_currents[cur_id][0], usable_currents[cur_id][0], usable_currents[cur_id][1], usable_currents[cur_id][1]], 2)
+        electrode_vals.to_csv(os.path.join(csv_dir, 'optimized_electrodes_' + model + '.csv'))
+        print('Value: {}, Current: {}, {}'.format(value, usable_currents[cur_id][0], usable_currents[cur_id][1]))
+        
+        del field_data
+        del npz_arrays
+        gc.collect()
diff --git a/Scripts/Gmsh/spheres.geo b/Scripts/Gmsh/spheres.geo
new file mode 100644
index 0000000..7a78a68
--- /dev/null
+++ b/Scripts/Gmsh/spheres.geo
@@ -0,0 +1,38 @@
+// Gmsh project created on Tue Oct 13 06:02:07 2020
+SetFactory("OpenCASCADE");
+
+head_radius = 87.53;
+brain_ratio = 0.83;
+csf_ratio = 0.023;
+skull_ratio = 0.085;
+scalp_ratio = 0.05;
+
+//Sphere(1) = {0, 0, 0, head_radius};
+//Sphere(2) = {0, 0, 0, head_radius*(brain_ratio + csf_ratio + skull_ratio + scalp_ratio)};
+//Sphere(3) = {0, 0, 0, head_radius*(brain_ratio + csf_ratio + skull_ratio)};
+//Sphere(4) = {0, 0, 0, head_radius*(brain_ratio + csf_ratio)};
+Sphere(5) = {0, 0, 0, head_radius*(brain_ratio)};
+
+//Physical Volume("Outer Boundary", 1) = {1};
+
+/*
+Physical Volume("Skin", 2) = {2};
+Physical Volume("Skull", 3) = {3};
+Physical Volume("CSF", 4) = {4};
+Physical Volume("Brain", 5) = {5};
+
+Physical Surface("Outer Boundary", 1) = {1};
+Physical Surface("Skin", 2) = {2};
+Physical Surface("Skull", 3) = {3};
+Physical Surface("CSF", 4) = {4};
+Physical Surface("Brain", 5) = {5};
+*/
+
+Mesh 2;
+RefineMesh;
+RefineMesh;
+RefineMesh;
+Mesh 3;
+RefineMesh;
+
+
diff --git a/Scripts/Nifti Operations/[Post-process]_Electrodes_to_Nifti.py b/Scripts/Nifti Operations/[Post-process]_Electrodes_to_Nifti.py
new file mode 100644
index 0000000..e3dedab
--- /dev/null
+++ b/Scripts/Nifti Operations/[Post-process]_Electrodes_to_Nifti.py	
@@ -0,0 +1,82 @@
+import os
+import gc
+import yaml
+import numpy as np
+import pyvista as pv
+import pandas as pd
+import nibabel as nib
+import PVGeo
+
+model_path_vtk = 'full path of the directory containing the VTK solved models'
+save_path = 'output directory to save the Nifti files'
+npz_path = ''
+
+files = np.sort(next(os.walk(npz_path))[2])
+
+with open(os.path.realpath('settings file path')) as stream:
+    settings = yaml.safe_load(stream)
+
+distance_margin = 15
+voxel_size = 1
+
+for fl in files:
+    model_id = os.path.splitext(fl)[0].split('_')[0]
+    print(model_id)
+    msh = pv.UnstructuredGrid(os.path.join(model_path_vtk, 'Meshed_AAL_10-10_' + model_id + '.vtk'))
+
+    mesh_pts = msh.threshold(value=[1, 8], scalars='cell_scalars', preference='cell').cell_centers()
+    pts_loc = np.isin(mesh_pts['cell_scalars'], [4, 5, 6])
+    mesh_points = mesh_pts.points[pts_loc]
+
+    bounds = msh.bounds
+    x_bounds_init = np.ceil(np.sum(np.abs(bounds[0:2]) + distance_margin)/voxel_size).astype(np.int32)
+    y_bounds_init = np.ceil(np.sum(np.abs(bounds[2:4]) + distance_margin)/voxel_size).astype(np.int32)
+    z_bounds_init = np.ceil(np.sum(np.abs(bounds[4:6]) + distance_margin)/voxel_size).astype(np.int32)
+
+    grid = pv.UniformGrid((x_bounds_init, y_bounds_init, z_bounds_init))
+    grid.origin = -1*np.array([np.sum(np.abs(bounds[0:2]))/2., np.sum(np.abs(bounds[2:4]))/2., np.sum(np.abs(bounds[4:6]))/2.]) - distance_margin
+    grid.spacing = [voxel_size]*3
+
+    img_header = nib.Nifti1Header()
+    img_header.set_xyzt_units('mm', 'sec')
+
+    npz_arrays = np.load(os.path.join(npz_path, fl), allow_pickle=True)
+    fields = npz_arrays['e_field']
+    electrodes = list(settings['SfePy']['electrodes']['10-10-mod'].items())
+    elec_interest = ['PO7', 'F7', 'P8', 'FC6']
+
+    for electrode_id, electrode_name in enumerate(elec_interest):
+        if electrode_name == 'P9':
+            continue
+
+        data = {'x': mesh_points[:, 0], 'y': mesh_points[:, 1], 'z': mesh_points[:, 2], 'int_x': fields[electrode_id, :, 0], 'int_y': fields[electrode_id, :, 1], 'int_z': fields[electrode_id, :, 2]}
+        df = pd.DataFrame(data, columns=['x', 'y', 'z', 'int_x', 'int_y', 'int_z'])
+        grid_points = PVGeo.points_to_poly_data(df)
+
+        voxels = grid.interpolate(grid_points, radius=voxel_size*4, sharpness=8)
+
+        for e_field_array in voxels.array_names:
+            voxel_data = voxels[e_field_array].reshape((z_bounds_init, y_bounds_init, x_bounds_init)).transpose()
+            direction = e_field_array.split('_')[1]
+
+            x_minus_bound = bounds[0] - distance_margin
+            y_minus_bound = bounds[2] - distance_margin
+            z_minus_bound = bounds[4] - distance_margin
+            affine = np.array([[voxel_size, 0, 0, x_minus_bound], [0, voxel_size, 0, y_minus_bound], [0, 0, voxel_size, z_minus_bound], [0, 0, 0, 1]])
+
+            try:
+                os.makedirs(os.path.join(save_path, model_id, electrode_name))
+            except OSError:
+                pass
+
+            img = nib.Nifti1Image(voxel_data, affine, img_header)
+            nib.save(img, os.path.join(save_path, model_id, electrode_name, 'Efield_{}_{}_{}.nii'.format(direction, electrode_name, model_id)))
+
+        del voxels
+        del voxel_data
+        gc.collect()
+
+    del msh
+    del npz_arrays
+    del fields
+    gc.collect()
diff --git a/Scripts/Nifti Operations/[Post-process]_FEM_to_Nifti.py b/Scripts/Nifti Operations/[Post-process]_FEM_to_Nifti.py
new file mode 100644
index 0000000..88d7314
--- /dev/null
+++ b/Scripts/Nifti Operations/[Post-process]_FEM_to_Nifti.py	
@@ -0,0 +1,36 @@
+import os
+import yaml
+import multiprocessing
+
+import NiiPHM.NiiMesh as NiiMesh
+
+
+with open(os.path.realpath('settings file path')) as stream:
+    settings = yaml.safe_load(stream)
+
+index_array_name = 'cell_scalars'
+
+model_path_vtk = 'full path of the directory containing the VTK solved models'
+files = next(os.walk(model_path_vtk))[2]
+
+conductivities = [1 for a in settings['SfePy']['real_brain']['regions'].values()]
+
+def nfts(fl):
+    global conductivities
+    model_path_vtk = 'full path of the directory containing the VTK solved models'
+    index_array_name = 'cell_scalars'
+    structures_to_remove = [1, 2, 3, 5, 6, 7]
+    
+    output_dir = 'output directory to save the Nifti files'
+    model_id = os.path.splitext(fl)[0].split('_')[-1].split('-')[0].split('.')[0]
+    print("Converting model {}".format(model_id))
+    
+    nifti = NiiMesh.NiiMesh(os.path.join(model_path_vtk, fl), index_array_name, [1, 8])
+    _ = nifti.assign_intensities(assign_index=False, unwanted_region_ids=structures_to_remove, assign_values_per_region=True, values_to_assign=conductivities)
+    _ = nifti.generate_uniform_grid(15, 0.75, )
+    nifti.generate_nifti(image_path=os.path.join(output_dir, 'PHM_' + model_id + '.nii'))
+
+
+if __name__ == '__main__':
+    pool = multiprocessing.Pool(processes=4)
+    pool.map(nfts, files)
diff --git a/Scripts/Nifti Operations/[Post-process]_Nifti_to_NPZ.py b/Scripts/Nifti Operations/[Post-process]_Nifti_to_NPZ.py
new file mode 100644
index 0000000..6b3fb0b
--- /dev/null
+++ b/Scripts/Nifti Operations/[Post-process]_Nifti_to_NPZ.py	
@@ -0,0 +1,84 @@
+import os
+import gc
+import yaml
+import numpy as np
+import nibabel as nib
+
+from argparse import ArgumentParser, RawDescriptionHelpFormatter
+
+#### Argument parsing
+helps = {
+    'settings-file' : "File having the settings to be loaded",
+    'aal-atlas' : "Full path of the AAL atlas Nifti",
+    'model-path' : "Full path of the dirctory which includes the models",
+}
+
+parser = ArgumentParser(description=__doc__,
+                        formatter_class=RawDescriptionHelpFormatter)
+parser.add_argument('--version', action='version', version='%(prog)s')
+parser.add_argument('--settings-file', metavar='str', type=str,
+                    action='store', dest='settings_file',
+                    default=None, help=helps['settings-file'], required=True)
+parser.add_argument('--aal-atlas', metavar='str', type=str,
+                    action='store', dest='atlas_path',
+                    default=None, help=helps['aal-atlas'], required=True)
+parser.add_argument('--model-path', metavar='str', type=str,
+                    action='store', dest='model_path',
+                    default='real_brain', help=helps['model-path'], required=True)
+parser.add_argument('--save-dir', metavar='str', type=str,
+                    action='store', dest='save_dir',
+                    default=None, required=True)
+options = parser.parse_args()
+#### Argument parsing
+
+folders = np.sort(next(os.walk(options.model_path))[1])
+e_field_size = []
+
+with open(os.path.realpath(options.settings_file)) as stream:
+    settings = yaml.safe_load(stream)
+
+atlas_nifti = nib.load(options.atlas_path)
+aal_regions_loc = np.where(atlas_nifti.get_fdata().flatten() != 0)[0].astype(int)
+aal_regions_ids = atlas_nifti.get_fdata().flatten()[aal_regions_loc]
+
+settings_electrodes = settings['SfePy']['electrodes']['10-10-mod'].keys()
+electrode_constant = 10
+
+for fld in folders:
+    e_field_values = []
+    model_id = fld
+    print(model_id)
+    
+    for electrode in settings_electrodes:
+        electrode_name = electrode
+        if electrode_name == 'P9':
+            continue
+        print(electrode_name)
+        files = np.sort(next(os.walk(os.path.join(options.model_path, model_id, electrode_name)))[2])
+        
+        nifti_images_dict = {}
+        for fl in files:
+            if not fl.startswith('w'):
+                continue
+            field_montage = fl.split('.')[0].split('_')[2]
+            field_direction = fl.split('.')[0].split('_')[1]
+            dict_key = '{}_{}'.format(field_montage, field_direction)
+            
+            nifti_images_dict[dict_key] = nib.load(os.path.join(options.model_path, model_id, electrode_name, fl))
+        
+        e_field = np.nan_to_num(np.vstack((nifti_images_dict[electrode_name + '_x'].get_fdata().flatten(), nifti_images_dict[electrode_name + '_y'].get_fdata().flatten(), nifti_images_dict[electrode_name + '_z'].get_fdata().flatten())).transpose())
+        e_field = e_field[aal_regions_loc].reshape((1, -1, 3))
+        
+        if isinstance(e_field_values, list):
+            e_field_values = e_field
+        else:
+            e_field_values = np.append(e_field_values, e_field, axis=0)
+        
+        del nifti_images_dict
+        del e_field
+        gc.collect()
+
+    np.savez_compressed(os.path.join(options.save_dir, model_id + '_fields_brain_reg'), e_field=e_field_values.reshape((4, -1, 3)), aal_ids=aal_regions_ids, aal_loc=aal_regions_loc)
+
+    del e_field_values
+    gc.collect()
diff --git a/Scripts/[Meshing]_Model_POLY.py b/Scripts/[Meshing]_Model_POLY.py
new file mode 100644
index 0000000..b57806f
--- /dev/null
+++ b/Scripts/[Meshing]_Model_POLY.py
@@ -0,0 +1,65 @@
+import os
+import yaml
+import pymesh
+import Meshing.MeshOperations as MeshOps
+
+import scipy.io as sio
+import numpy as np
+
+from argparse import ArgumentParser, RawDescriptionHelpFormatter
+
+#### Argument parsing
+helps = {
+    'settings-file' : "File having the settings to be loaded",
+    'models-dir' : "Full path of the model directory",
+    'electrode-attributes-path' : "Full path of the Mesh2EEG generated electrodes file",
+}
+
+parser = ArgumentParser(description=__doc__,
+                        formatter_class=RawDescriptionHelpFormatter)
+parser.add_argument('--version', action='version', version='%(prog)s')
+parser.add_argument('--settings-file', metavar='str', type=str,
+                    action='store', dest='settings_file',
+                    default=None, help=helps['settings-file'], required=True)
+parser.add_argument('--models-dir', metavar='str', type=str,
+                    action='store', dest='stl_models_path',
+                    default=None, help=helps['models-dir'], required=True)
+parser.add_argument('--electrode-attributes-path', metavar='str', type=str,
+                    action='store', dest='electrode_attributes_path',
+                    default='real_brain', help=helps['electrode-attributes-path'], required=True)
+parser.add_argument('--save-dir', metavar='str', type=str,
+                    action='store', dest='save_dir',
+                    default=None, required=True)
+options = parser.parse_args()
+#### Argument parsing
+
+folders = np.sort(next(os.walk(options.stl_models_path))[1])
+
+# Electrodes to ommit based on the imported settings file
+electrodes_to_omit=['Nz', 'N2', 'AF10', 'F10', 'FT10', 'T10(M2)', 'TP10', 'PO10', 'I2', 'Iz', 'I1', 'PO9', 'TP9', 'T9(M1)', 'FT9', 'F9', 'AF9', 'N1', 'P10']
+
+with open(os.path.realpath(options.settings_file)) as stream:
+    settings = yaml.safe_load(stream)
+
+if __name__ == '__main__':
+    for folder in folders:
+        if folder == 'meshed':
+            continue
+        print("############")
+        print("Model " + folder)
+        print("############\n")
+        standard_electrodes = sio.loadmat(os.path.join(options.electrode_attributes_path, '10-10_elec_' + folder + '.mat'))
+        elec_attributes = {
+            'names': [name[0][0] for name in standard_electrodes['ElectrodeNames']],
+            'coordinates': standard_electrodes['ElectrodePts'],
+            'ids': settings['SfePy']['electrodes']['10-10-mod'],
+            'width': 4,
+            'radius': 4,
+            'elements': 200,
+        }
+
+        skin_stl = pymesh.load_mesh(os.path.join(options.stl_models_path, folder, 'skin_fixed.stl'))
+
+        meshing = MeshOps.MeshOperations(skin_stl, elec_attributes)
+        meshing.load_surface_meshes(os.path.join(options.stl_models_path, folder), ['skin_fixed.stl', 'skull_fixed.stl', 'csf_fixed.stl', 'gm_fixed.stl', 'wm_fixed.stl', 'cerebellum_fixed.stl', 'ventricles_fixed.stl'])
+        meshing.phm_model_meshing(os.path.join(options.save_dir, 'meshed_model_10-10_' + folder + '.poly'), electrodes_to_omit=electrodes_to_omit)
diff --git a/Scripts/[Optimization]_Data_Analysis.py b/Scripts/[Optimization]_Data_Analysis.py
new file mode 100644
index 0000000..9e3333b
--- /dev/null
+++ b/Scripts/[Optimization]_Data_Analysis.py
@@ -0,0 +1,146 @@
+import os
+import gc
+import numpy as np
+import pandas as pd
+
+def modulation_envelope(e_field_1, e_field_2):
+    envelope = np.zeros(e_field_1.shape[0])
+    # Calculate the angles between the two fields for each vector
+    dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
+    cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
+    angles = np.arctan2(cross_angle, dot_angle)
+    # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
+    e_field_2 = np.where(np.broadcast_to(angles >= np.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
+    # Recalculate the angles
+    dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
+    cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
+    angles = np.arctan2(cross_angle, dot_angle)
+    E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
+    # Condition to have two times the E2 field amplitude
+    max_condition_1 = np.linalg.norm(e_field_2, axis=1) < np.linalg.norm(e_field_1, axis=1)*np.cos(angles)
+    e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), max_condition_1, False)
+    # Condition to have two times the E1 field amplitude
+    max_condition_2 = np.linalg.norm(e_field_1, axis=1) < np.linalg.norm(e_field_2, axis=1)*np.cos(angles)
+    e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), max_condition_2, False)
+    # Double magnitudes
+    envelope = np.where(e1_gr_e2, 2.0*np.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
+    envelope = np.where(e2_gr_e1, 2.0*np.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
+    # Calculate the complement area to the previous calculation
+    e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), np.logical_not(max_condition_1), False)
+    e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), np.logical_not(max_condition_2), False)
+    # Cross product
+    envelope = np.where(e1_gr_e2, 2.0*(np.linalg.norm(np.cross(e_field_2, E_minus), axis=1)/np.linalg.norm(E_minus, axis=1)), envelope) # (First case)
+    envelope = np.where(e2_gr_e1, 2.0*(np.linalg.norm(np.cross(e_field_1, -E_minus), axis=1)/np.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
+    return np.nan_to_num(envelope)
+
+UNOPTIMIZED = False
+MAX_E_FIELD_CALC = False
+
+if MAX_E_FIELD_CALC:
+    def region_calc(modulation_values, aal_regions, region_volumes, regions_of_interest, ommit=None):
+        roi = np.isin(aal_regions, regions_of_interest)    
+        return np.amax(modulation_values[roi])
+else:
+    def region_calc(modulation_values, aal_regions, region_volumes, regions_of_interest, ommit=None):
+        roi_region_sum = 0
+        non_roi_region_sum = 0
+        
+        for region in np.unique(aal_regions):
+            roi = np.where(aal_regions == region)[0]
+            #if region in ommit:
+            #        continue
+            if int(region) in regions_of_interest:
+                roi_region_sum += np.sum(modulation_values[roi])/region_volumes[int(region)]
+            else:
+                non_roi_region_sum += np.sum(modulation_values[roi])/region_volumes[int(region)]
+        
+        return np.nan_to_num(roi_region_sum/non_roi_region_sum)
+
+
+if __name__ == "__main__":
+    npz_dir = 'directory of the NPZ model files'
+    csv_dir = 'directory of the optimized electrode pairs'
+    files = next(os.walk(csv_dir))[2]
+    
+    models = np.sort(next(os.walk(csv_dir))[2])
+    sum_values_df = {'TS': []}
+    focal_values = {'TS': []}
+    focal_values_df = {'TS': []}
+    
+    if MAX_E_FIELD_CALC:
+        mult_factor = 1
+    else:
+        mult_factor = 10000
+    
+    for model in models:
+        model = model.split('.')[0].split('_')[-1]
+        npz_arrays = np.load(os.path.join(npz_dir, model + '_fields_brain_reg.npz'), allow_pickle=True)
+        electrode_vals = pd.read_csv(os.path.join(csv_dir, 'optimized_electrodes_' + model + '.csv'))
+        
+        field_data = npz_arrays['e_field']
+        aal_regions = npz_arrays['aal_ids']
+        aal_regions[np.isin(aal_regions, np.arange(122, 150, 2))] = 122
+        aal_regions[np.isin(aal_regions, np.arange(121, 150, 2))] = 121
+        region_volumes = {}
+        for region in np.unique(aal_regions):
+            region_volumes[region] = np.where(aal_regions == region)[0].size
+        
+        print(model)
+        
+        if UNOPTIMIZED:
+            electrodes_model = [53, 8, 52, 24]
+            currents = [1.25, 1.25, 0.75, 0.75]
+        else:
+            electrodes_model = electrode_vals['electrodes'].to_numpy().astype(int)
+            currents = electrode_vals['currents'].to_numpy()
+        
+        base_e_field = currents[0]*(field_data[electrodes_model[0]] - field_data[electrodes_model[1]])
+        df_e_field = currents[2]*(field_data[electrodes_model[2]] - field_data[electrodes_model[3]])
+        
+        max_mod = modulation_envelope(base_e_field, df_e_field)
+        
+        sum_values_df[model] = {}
+        for region in np.unique(aal_regions):
+            roi = np.where(aal_regions == region)[0]
+            sum_values_df[model][str(int(region))] = np.sum(max_mod[roi])/float(roi.size)
+        
+        ## Right area  
+        thal_right = region_calc(max_mod, aal_regions, region_volumes, np.arange(122, 150, 2), [42])
+        hip_right = region_calc(max_mod, aal_regions, region_volumes, [42], np.arange(122, 150, 2))
+        focal_values[model] = {'Thal_Right': thal_right*mult_factor}
+        focal_values[model]['Hipp_Right'] = hip_right*mult_factor
+        ## Right area
+        
+        ## Left area      
+        thal_left = region_calc(max_mod, aal_regions, region_volumes, np.arange(121, 150, 2), [41])
+        hip_left = region_calc(max_mod, aal_regions, region_volumes, [41], np.arange(121, 150, 2))
+        focal_values[model]['Thal_Left'] = thal_left*mult_factor
+        focal_values[model]['Hipp_Left'] = hip_left*mult_factor
+        ## Left area
+        
+        ## Both areas
+        thal_both = region_calc(max_mod, aal_regions, region_volumes, np.arange(121, 151, 1), [41, 42])
+        hip_both = region_calc(max_mod, aal_regions, region_volumes, [41, 42], np.arange(121, 151, 1))
+        focal_values[model]['Thal_Both'] = thal_both*mult_factor
+        focal_values[model]['Hipp_Both'] = hip_both*mult_factor
+        ## Both areas
+        
+        ## Cortex
+        cortex = region_calc(max_mod, aal_regions, region_volumes, [-4], [-5])
+        wm = region_calc(max_mod, aal_regions, region_volumes, [-5], [-4])
+        focal_values[model]['Cortex'] = cortex*mult_factor
+        focal_values[model]['White_Matter'] = wm*mult_factor
+        ## Cortex
+        
+        del npz_arrays
+        del field_data
+        del base_e_field
+        del df_e_field
+        del max_mod
+        gc.collect()
+
+sum_values_df.pop('TS')
+focal_values.pop('TS')
+
+np.savez('input the save path', sum_values_df=sum_values_df, focal_values=focal_values)
+
diff --git a/Scripts/[Optimization]_Gentic_Algorithm.py b/Scripts/[Optimization]_Gentic_Algorithm.py
new file mode 100644
index 0000000..59bb257
--- /dev/null
+++ b/Scripts/[Optimization]_Gentic_Algorithm.py
@@ -0,0 +1,136 @@
+import os
+import gc
+import sys
+import yaml
+import time
+import numpy as np
+import pandas as pd
+import cupy as cp
+from geneticalgorithm import geneticalgorithm as ga
+
+from argparse import ArgumentParser, RawDescriptionHelpFormatter
+
+#### Argument parsing
+helps = {
+    'settings-file' : "File having the settings to be loaded",
+    'model' : "Name of the model. Selection from the settings file",
+}
+
+parser = ArgumentParser(description=__doc__,
+                        formatter_class=RawDescriptionHelpFormatter)
+parser.add_argument('--version', action='version', version='%(prog)s')
+parser.add_argument('--settings-file', metavar='str', type=str,
+                    action='store', dest='settings_file',
+                    default=None, help=helps['settings-file'], required=True)
+parser.add_argument('--npz-file', metavar='str', type=str,
+                    action='store', dest='npz_file',
+                    default='real_brain', help=helps['model'], required=True)
+parser.add_argument('--save-dir', metavar='str', type=str,
+                    action='store', dest='save_dir',
+                    default=None, required=True)
+options = parser.parse_args()
+#### Argument parsing
+
+with open(os.path.realpath(options.settings_file)) as stream:
+    settings = yaml.safe_load(stream)
+
+if os.name == 'nt':
+    extra_path = '_windows'
+else:
+    extra_path = ''
+
+sys.path.append(os.path.realpath(settings['SfePy']['lib_path' + extra_path]))
+
+from modulation_envelope import modulation_envelope_gpu
+
+def objective_df(x, field_data, regions_of_interest, aal_regions, region_volumes, currents, threshold):  
+    if np.unique(np.round(x[:4]), return_counts=True)[1].size != 4:
+        return 100*(np.abs((np.unique(np.round(x[:4]), return_counts=True)[1].size - 4))**2) + 10000
+    
+    penalty = 0
+    electrodes = np.round(x[:4]).astype(np.int32) # The first 4 indices are the electrode IDs
+    roi = cp.isin(aal_regions, cp.array(regions_of_interest))
+    max_vals = []
+    fitness_vals = []
+    
+    for current in currents:
+        e_field_base = current[0]*field_data[electrodes[0]] - current[0]*field_data[electrodes[1]]
+        e_field_df = current[1]*field_data[electrodes[2]] - current[1]*field_data[electrodes[3]]
+
+        e_field_base_gpu = cp.array(e_field_base)
+        e_field_df_gpu = cp.array(e_field_df)
+        modulation_values = modulation_envelope_gpu(e_field_base_gpu, e_field_df_gpu)
+        max_vals.append(float(cp.amax(modulation_values[roi])))
+    
+        roi_region_sum = 0
+        non_roi_region_sum = 0
+
+        for region in cp.unique(aal_regions):
+            roi = cp.where(aal_regions == region)[0]
+        
+            if int(region) in regions_of_interest:
+                roi_region_sum += cp.sum(modulation_values[roi])/region_volumes[int(region)]
+            else:
+                non_roi_region_sum += cp.sum(modulation_values[roi])/region_volumes[int(region)]
+        
+        region_ratio = cp.nan_to_num(roi_region_sum/non_roi_region_sum)
+        fitness_measure = region_ratio*10000
+        fitness_vals.append(float(fitness_measure))
+    
+    max_vals = np.array(max_vals)
+    max_val_curr = np.where(max_vals >= threshold)[0]
+    fitness_vals = np.array(fitness_vals)
+
+    return_fitness = 0
+    if max_val_curr.size == 0:
+        penalty += 100*((threshold - np.mean(max_vals))**2) + 1000
+        return_fitness = np.amin(fitness_vals)
+    else:
+        fitness_candidate = np.amax(fitness_vals[max_val_curr])
+        return_fitness = fitness_candidate
+    
+    return -float(np.round(return_fitness - penalty, 2))
+
+if __name__ == "__main__":
+    OPTIMIZATION_THRESHOLD = 0.2
+    
+    npz_arrays = np.load(options.npz_file, allow_pickle=True)
+    
+    field_data = npz_arrays['e_field']
+    aal_regions = npz_arrays['aal_ids']
+    region_volumes = {}
+    for region in np.unique(aal_regions):
+        region_volumes[region] = np.where(aal_regions == region)[0].size
+    roi_ids = np.array([42])
+    
+    algorithm_param = {'max_num_iteration': 25,
+                       'population_size': 100,
+                       'mutation_probability': 0.4,
+                       'elit_ratio': 0.01,
+                       'crossover_probability': 0.5,
+                       'parents_portion': 0.1,
+                       'crossover_type': 'uniform',
+                       'max_iteration_without_improv': None
+                    }
+    
+    cur_potential_values = np.arange(0.5, 1.55, 0.05)
+    cur_x, cur_y = np.meshgrid(cur_potential_values, cur_potential_values)
+    cur_all_combinations = np.hstack((cur_x.reshape((-1, 1)), cur_y.reshape((-1, 1))))
+    usable_currents = cur_all_combinations[np.where(np.sum(np.round(cur_all_combinations, 2), axis=1) == 2)[0]]
+    
+    aal_regions_gpu = cp.array(aal_regions)
+    ga_objective_df = lambda x, **kwargs: objective_df(x, field_data, roi_ids, aal_regions=aal_regions_gpu, region_volumes=region_volumes, currents=usable_currents, threshold=OPTIMIZATION_THRESHOLD)
+    
+    var_type = np.array([['int']]*4)
+    bounds = np.array([[0, 60]]*4)
+    result = ga(function=ga_objective_df, dimension=bounds.shape[0], variable_type_mixed=var_type, variable_boundaries=bounds, algorithm_parameters=algorithm_param, function_timeout=120., convergence_curve=False)
+    result.run()
+
+    convergence = result.report
+    solution = result.output_dict
+
+    model_id = os.path.basename(options.npz_file).split('.')[0].split('_')[0]
+    df_dict = {'electrodes': solution['variable'], 'value': solution['function']}
+    df = pd.DataFrame(df_dict)
+
+    df.to_csv(os.path.join(options.save_dir, 'optimized_electrodes_' + model_id + '.csv'))
diff --git a/Scripts/[Post-process]_Electrode_Combinations.py b/Scripts/[Post-process]_Electrode_Combinations.py
new file mode 100644
index 0000000..7c01fce
--- /dev/null
+++ b/Scripts/[Post-process]_Electrode_Combinations.py
@@ -0,0 +1,102 @@
+from __future__ import absolute_import
+import os
+import gc
+import sys
+import yaml
+import numpy as np
+import pyvista as pv
+
+from argparse import ArgumentParser, RawDescriptionHelpFormatter
+
+#### Argument parsing
+helps = {
+    'settings-file' : "File having the settings to be loaded",
+    'model' : "Name of the model. Selection from the settings file",
+}
+
+parser = ArgumentParser(description=__doc__,
+                        formatter_class=RawDescriptionHelpFormatter)
+parser.add_argument('--version', action='version', version='%(prog)s')
+parser.add_argument('--settings-file', metavar='str', type=str,
+                    action='store', dest='settings_file',
+                    default=None, help=helps['settings-file'], required=True)
+parser.add_argument('--meshf', metavar='str', type=str,
+                    action='store', dest='meshf',
+                    default=None, required=True)
+parser.add_argument('--model', metavar='str', type=str,
+                    action='store', dest='model',
+                    default='real_brain', help=helps['model'], required=True)
+parser.add_argument('--csv_save_dir', metavar='str', type=str,
+                    action='store', dest='csv_save_dir',
+                    default=None, required=True)
+parser.add_argument('--job_id', metavar='str', type=str,
+                    action='store', dest='job_id',
+                    default='', required=False)
+options = parser.parse_args()
+#### Argument parsing
+
+with open(os.path.realpath(options.settings_file)) as stream:
+    settings = yaml.safe_load(stream)
+
+if os.name == 'nt':
+    extra_path = '_windows'
+else:
+    extra_path = ''
+
+sys.path.append(os.path.realpath(settings['SfePy']['lib_path' + extra_path]))
+
+import FEM.Solver as slv
+
+settings['SfePy'][options.model]['mesh_file' + extra_path] = options.meshf
+model_id = os.path.basename(options.meshf).split('.')[0].split('_')[-1]
+
+## Read the mesh with pyvista to get the area ids and AAL regions
+msh = pv.UnstructuredGrid(options.meshf)
+brain_regions_mask = np.isin(msh['cell_scalars'], [4, 5, 6])
+
+cell_ids_brain = msh['cell_scalars'][brain_regions_mask]
+aal_regions = msh['AAL_regions'][brain_regions_mask]
+cell_volumes_brain = msh.compute_cell_sizes().cell_arrays['Volume'][brain_regions_mask]
+
+region_volumes_brain = {}
+for region in np.unique(aal_regions):
+    roi = np.where(aal_regions == region)[0]
+    region_volumes_brain[int(region)] = np.sum(cell_volumes_brain[roi])
+del msh
+region_volumes_brain = np.array(region_volumes_brain)
+## Read the mesh with pyvista to get the area ids and AAL regions
+
+electrodes = settings['SfePy']['electrodes']['10-10-mod']
+e_field_values_brain = []
+
+solve = slv.Solver(settings, 'SfePy', '10-10-mod')
+solve.load_mesh(options.model)
+
+for electrode in electrodes.items():
+    if electrode[0] == 'P9':
+        continue
+    solve.essential_boundaries.clear()
+    solve.fields.clear()
+    solve.field_variables.clear()
+
+    solve.define_field_variable('potential', 'voltage', out_of_range_assign_region='Skin', out_of_range_group_threshold=71)
+
+    solve.define_essential_boundary(electrode[0], electrode[1]['id'], 'potential', current=1)
+    solve.define_essential_boundary('P9', 71, 'potential', current=-1)
+
+    solve.solver_setup(600, 1e-12, 5e-12, verbose=True)
+    solution = solve.run_solver(save_results=False, post_process_calculation=True)
+
+    e_field_base = solution['e_field_(potential)'].data[:, 0, :, 0]
+    if isinstance(e_field_values_brain, list):
+        e_field_values_brain = e_field_base[brain_regions_mask]
+    else:
+        e_field_values_brain = np.append(e_field_values_brain, e_field_base[brain_regions_mask], axis=0)
+
+    del solution
+    gc.collect()
+
+del solve
+gc.collect
+
+np.savez_compressed(os.path.join(options.csv_save_dir, model_id + '_fields_brain'), e_field=e_field_values_brain.reshape((61, -1, 3)), cell_ids=cell_ids_brain, aal_regions=aal_regions, volumes=region_volumes_brain)
diff --git a/Utilities/mesh_fixing.py b/Utilities/mesh_fixing.py
new file mode 100644
index 0000000..6d2b68f
--- /dev/null
+++ b/Utilities/mesh_fixing.py
@@ -0,0 +1,24 @@
+import os
+import gc
+import multiprocessing
+import meshio
+
+base_path = "full path of the model's directory"
+mesh_fix_path = "full path of meshfix executable"
+folders = next(os.walk(base_path))[1]
+
+def fixing(folder):
+    global base_path
+    for file_n in next(os.walk(os.path.join(base_path, folder)))[2]:
+        print("\nConverting to ASCII file")
+        temp_file = meshio.read(os.path.join(base_path, folder, file_n))
+        temp_file.write(os.path.join(base_path, folder, file_n))
+        del temp_file
+
+        print("Fixing")
+        os.system(mesh_fix_path + " " + os.path.join(base_path, folder, file_n) + r" " + os.path.join(base_path, folder, file_n.split('.')[0] + '_fixed.stl') + r" -j")
+    return 0
+
+if __name__ == '__main__':
+    pool = multiprocessing.Pool(processes=8)
+    pool.map_async(fixing, folders)
diff --git a/Utilities/tetgen_meshing.sh b/Utilities/tetgen_meshing.sh
new file mode 100644
index 0000000..125b14e
--- /dev/null
+++ b/Utilities/tetgen_meshing.sh
@@ -0,0 +1,11 @@
+#!/usr/bin/env bash
+
+set -e
+
+cd 'to model directory'
+
+for file in *.poly; do 
+    if [ -f "$file" ]; then 
+        time 'tetgen executable path' -zpq1.8/0O4a30kNEFAV "$file"
+    fi 
+done
diff --git a/sim_settings.yml b/sim_settings.yml
new file mode 100644
index 0000000..876f588
--- /dev/null
+++ b/sim_settings.yml
@@ -0,0 +1,368 @@
+SfePy:
+  lib_path: 'full path of the python code root directory on linux'
+  lib_path_windows: 'full path of the python code root directory on windows'
+  electrodes:
+    conductivity: 1.4
+    10-20:
+      Fp1:
+        id: 10
+      Fp2:
+        id: 11
+      F7:
+        id: 12
+      F3:
+        id: 13
+      Fz:
+        id: 14
+      F4:
+        id: 15
+      F8:
+        id: 16
+      T7:
+        id: 17
+      C3:
+        id: 18
+      Cz:
+        id: 19
+      C4:
+        id: 20
+      T8:
+        id: 21
+      P7:
+        id: 22
+      P3:
+        id: 23
+      Pz:
+        id: 24
+      P4:
+        id: 25
+      P8:
+        id: 26
+      O1:
+        id: 27
+      O2:
+        id: 28
+    10-10:
+      Nz:
+        id: 10
+      N1:
+        id: 11
+      Fp1:    
+        id: 12
+      Fpz:    
+        id: 13
+      Fp2:    
+        id: 14
+      N2:     
+        id: 15
+      AF9:    
+        id: 16
+      AF7:    
+        id: 17
+      AF3:    
+        id: 18
+      AFz:    
+        id: 19
+      AF4:    
+        id: 20
+      AF8:    
+        id: 21
+      AF10:   
+        id: 22
+      F9:     
+        id: 23
+      F7:     
+        id: 24
+      F5:     
+        id: 25
+      F3:     
+        id: 26
+      F1:     
+        id: 27
+      Fz:     
+        id: 28
+      F2:     
+        id: 29
+      F4:     
+        id: 30
+      F6:     
+        id: 31
+      F8:     
+        id: 32
+      F10:    
+        id: 33
+      FT9:    
+        id: 34
+      FT7:    
+        id: 35
+      FC5:    
+        id: 36
+      FC3:    
+        id: 37
+      FC1:    
+        id: 38
+      FCz:    
+        id: 39
+      FC2:    
+        id: 40
+      FC4:    
+        id: 41
+      FC6:    
+        id: 42
+      FT8:    
+        id: 43
+      FT10:   
+        id: 44
+      T9(M1): 
+        id: 45
+      T7:     
+        id: 46
+      C5:     
+        id: 47
+      C3:     
+        id: 48
+      C1:     
+        id: 49
+      Cz:     
+        id: 50
+      C2:     
+        id: 51
+      C4:     
+        id: 52
+      C6:     
+        id: 53
+      T8:     
+        id: 54
+      T10(M2):
+        id: 55
+      TP9:    
+        id: 56
+      TP7:    
+        id: 57
+      CP5:    
+        id: 58
+      CP3:    
+        id: 59
+      CP1:    
+        id: 60
+      CPz:    
+        id: 61
+      CP2:    
+        id: 62
+      CP4:    
+        id: 63
+      CP6:    
+        id: 64
+      TP8:    
+        id: 65
+      TP10:   
+        id: 66
+      P9:     
+        id: 67
+      P7:     
+        id: 68
+      P5:     
+        id: 69
+      P3:     
+        id: 70
+      P1:     
+        id: 71
+      Pz:     
+        id: 72
+      P2:     
+        id: 73
+      P4:     
+        id: 74
+      P6:     
+        id: 75
+      P8:     
+        id: 76
+      P10:    
+        id: 77
+      PO9:    
+        id: 78
+      PO7:    
+        id: 79
+      PO3:    
+        id: 80
+      POz:    
+        id: 81
+      PO4:    
+        id: 82
+      PO8:    
+        id: 83
+      PO10:   
+        id: 84
+      I1:     
+        id: 85
+      O1:     
+        id: 86
+      Oz:     
+        id: 87
+      O2:     
+        id: 88
+      I2:     
+        id: 89
+      Iz:     
+        id: 90
+    10-10-mod:
+      Fp1:    
+        id: 10
+      Fpz:    
+        id: 11
+      Fp2:    
+        id: 12
+      AF7:    
+        id: 13
+      AF3:    
+        id: 14
+      AFz:    
+        id: 15
+      AF4:    
+        id: 16
+      AF8:    
+        id: 17
+      F7:     
+        id: 18
+      F5:     
+        id: 19
+      F3:     
+        id: 20
+      F1:     
+        id: 21
+      Fz:     
+        id: 22
+      F2:     
+        id: 23
+      F4:     
+        id: 24
+      F6:     
+        id: 25
+      F8:     
+        id: 26
+      FT7:    
+        id: 27
+      FC5:    
+        id: 28
+      FC3:    
+        id: 29
+      FC1:    
+        id: 30
+      FCz:    
+        id: 31
+      FC2:    
+        id: 32
+      FC4:    
+        id: 33
+      FC6:    
+        id: 34
+      FT8:    
+        id: 35
+      T7:     
+        id: 36
+      C5:     
+        id: 37
+      C3:     
+        id: 38
+      C1:     
+        id: 39
+      Cz:     
+        id: 40
+      C2:     
+        id: 41
+      C4:     
+        id: 42
+      C6:     
+        id: 43
+      T8:     
+        id: 44
+      TP7:    
+        id: 45
+      CP5:    
+        id: 46
+      CP3:    
+        id: 47
+      CP1:    
+        id: 48
+      CPz:    
+        id: 49
+      CP2:    
+        id: 50
+      CP4:    
+        id: 51
+      CP6:    
+        id: 52
+      TP8:    
+        id: 53
+      P7:     
+        id: 54
+      P5:     
+        id: 55
+      P3:     
+        id: 56
+      P1:     
+        id: 57
+      Pz:     
+        id: 58
+      P2:     
+        id: 59
+      P4:     
+        id: 60
+      P6:     
+        id: 61
+      P8:     
+        id: 62
+      PO7:    
+        id: 63
+      PO3:    
+        id: 64
+      POz:    
+        id: 65
+      PO4:    
+        id: 66
+      PO8:    
+        id: 67
+      O1:     
+        id: 68
+      Oz:     
+        id: 69
+      O2:     
+        id: 70
+      P9:
+        id: 71
+    sphere:
+      B_VCC: 
+        id: 12
+      B_GND: 
+        id: 13
+      D_VCC: 
+        id: 10
+      D_GND: 
+        id: 11
+  real_brain:
+    mesh_file: 'full path of the mesh file on linux'
+    output_dir: 'full path of the output directory on linux'
+    mesh_file_windows: 'full path of the mesh file on windows'
+    output_dir_windows: 'full path of the mesh file on windows'
+    regions:
+      Skin:
+        id: 1
+        conductivity: 1 # Conductivity values in S/m
+      Skull:
+        id: 2
+        conductivity: 1
+      CSF:
+        id: 3
+        conductivity: 1
+      GM:
+        id: 4
+        conductivity: 1
+      WM:
+        id: 5
+        conductivity: 1
+      Cerebellum:
+        id: 6
+        conductivity: 1
+      Ventricles:
+        id: 7
+        conductivity: 1