Source code for nrv.fmod.FEM._FENICS_model

"""
NRV-:class:`.FENICS_model` handling.
"""

import configparser
import os
import time

import numpy as np
from mpi4py import MPI

from ...backend._file_handler import rmv_ext
from ...utils._units import V, mm
from ...utils._misc import get_perineurial_thickness
from ...backend._MCore import MCH, synchronize_processes
from ...backend._log_interface import rise_warning
from ...backend._parameters import parameters
from .fenics_utils._FEMSimulation import FEMSimulation

from ._FEM import FEM_model
from .fenics_utils._FEMResults import save_sim_res_list
from .mesh_creator._NerveMshCreator import NerveMshCreator, ENT_DOM_offset, pi

# built in FENICS models
dir_path = parameters.nrv_path + "/_misc"
# material_library = os.listdir(dir_path+"/fenics_templates/")

###############
## Constants ##
###############
machine_config = configparser.ConfigParser()
config_fname = dir_path + "/NRV.ini"
machine_config.read(config_fname)

FENICS_Ncores = int(machine_config.get("FENICS", "FENICS_CPU"))
FENICS_Status = machine_config.get("FENICS", "FENICS_STATUS") == "True"

fem_verbose = True


[docs] def check_sim_dom(sim: FEMSimulation, mesh: NerveMshCreator): """ Check if all the physical domains of a mesh are linked to a material property. This function can either prevent errors or help to debbug Paramters --------- sim: FEMSimulation Simulation to check. mesh: NerveMshCreator Mesh in used for the silmulation. Returns ------- bool True if all domains are linked to a material, False otherwise. """ uset_dom = [] mdom = mesh.domains_3D if sim.D == 2: mdom = mesh.domains_2D for i in mdom: if not i in sim.domains_list: uset_dom += [i] if len(uset_dom): rise_warning(f" the following domains are not set {uset_dom}") return False return True
[docs] class FENICS_model(FEM_model): """ A class for FENICS Finite Element Models, inherits from FEM_model. """
[docs] def __init__( self, fname=None, Ncore=None, elem=None, comm="default", rank=0, ): """ Creates a instance of a FENICS Finite Element Model object. Parameters ---------- fname : str path to the mesh file (.msh) file Ncore : int number of FENICS cores for computation. If None is specified, this number is taken from the NRV2.ini configuration file. Byu default set to None handle_server : bool if True, the instantiation creates the server, else a external server is used. Usefull for multiple cells sharing the same model """ super().__init__(Ncore=Ncore) self.type = "FENICS" # Default paramerters self.L = 5000 self.y_c = 0 self.z_c = 0 self.Outer_D = 5 # mm self.Nerve_D = 500 # um self.fascicles = {} self.Perineurium_thickness = {} self.electrodes = {} self.i_stim = 1e-3 # A self.jstims = [] self.j_electrode = {} self.Outer_mat = "saline" self.Epineurium_mat = "epineurium" self.Endoneurium_mat = "endoneurium_ranck" self.Perineurium_mat = "perineurium" self.Electrodes_mat = 1 # "platinum" self.default_fascicle = {"d": 250, "y_c": 0, "z_c": 0, "res": 250 / 3} self.default_electrode = { "elec_type": "LIFE", "x_c": self.L / 2, "y_c": 0, "z_c": 0, "length": 1000, "d": 25, "res": 3, } # Mesh self.mesh_file = fname self.mesh = None self.elem = elem if elem is None: self.elem = ("Lagrange", 2) # FEM self.sim = None self.sim_res = [] # Mcore attributes if self.Ncore is None: self.Ncore = FENICS_Ncores self.comm = None self.set_Ncore(self.Ncore) if comm != "default": self.comm = comm self.rank = rank # Status self.mesh_file_status = not self.mesh_file is None self.is_sim_ready = False self.is_meshed = False self.is_computed = False if not self.mesh_file_status: self.mesh = NerveMshCreator( Length=self.L, Outer_D=self.Outer_D, Nerve_D=self.Nerve_D, y_c=self.y_c, z_c=self.z_c, )
@property def N_fascicle(self): return len(self.fascicles) @property def N_electrode(self): return len(self.electrodes)
[docs] def save(self, save=False, fname="Fenics_model.json", blacklist=[], **kwargs): """ Save the changes to the model file. (Avoid for the overal weight of the package) """ bl = [i for i in blacklist] bl += ["comm", "rank", "mesh", "sim", "sim_res"] return super().save(save=save, fname=fname, blacklist=bl, **kwargs)
[docs] def save_results(self, save=False, fname="Fenics_model.json"): """ Save the changes to the model file. (Avoid for the overal weight of the package) """ if self.is_computed: save_sim_res_list(self.sim_res, fname) elif self.is_meshed: self.mesh_file = rmv_ext(fname) self.mesh.save(fname)
[docs] def load(self, data="extracel_context.json", **kwargs): """ Load all extracellular context properties from a dictionary or a json file Parameters ---------- data : str or dict json file path or dictionary containing extracel_context information """ bl = ["comm", "rank", "mesh", "sim", "sim_res"] super().load(data=data, blacklist=bl, **kwargs) param = self.get_parameters() self.load_from_parameters(**param)
[docs] def load_from_parameters(self, **param): """ """ self.reset_parameters() self.reshape_outerBox(param["Outer_D"]) self.reshape_nerve( Nerve_D=param["Nerve_D"], Length=param["L"], y_c=param["y_c"], z_c=param["z_c"], ) for id in param["fascicles"]: fasc = param["fascicles"][id] per_th = param["Perineurium_thickness"][id] self.reshape_fascicle( fasc["d"], y_c=fasc["y_c"], z_c=fasc["z_c"], ID=int(id), Perineurium_thickness=per_th, res=fasc["res"], ) for id in param["electrodes"]: elec = param["electrodes"][id] self.add_electrode( elec_type=elec["type"], ID=int(id), res=elec["res"], **elec["kwargs"] )
############################# ## Access model parameters ## #############################
[docs] def set_Ncore(self, N): """ Set the number of cores to use for the FEM Parameters ---------- N : int Number of cores to set """ super().set_Ncore(N) if self.is_multi_proc: self.comm = MPI.COMM_WORLD else: self.comm = MPI.COMM_SELF
[docs] def get_parameters(self): """ Get the all the parameters in the model as a python dictionary. Returns ------- param : dict all parameters as dictionnaries in a list, names a keys, with corresponding values """ param = {} param["L"] = self.L param["y_c"] = self.y_c param["z_c"] = self.z_c param["Outer_D"] = self.Outer_D / mm param["Nerve_D"] = self.Nerve_D param["fascicles"] = self.fascicles param["Perineurium_thickness"] = self.Perineurium_thickness param["N_electrode"] = self.N_electrode param["electrodes"] = self.electrodes param["Outer_mat"] = self.Outer_mat param["Epineurium_mat"] = self.Epineurium_mat param["Endoneurium_mat"] = self.Endoneurium_mat param["Perineurium_mat"] = self.Perineurium_mat param["i_stim"] = self.i_stim return param
def __update_parameters(self, bcast=False): """ Internal use only: updates all the parameters from the mesh Parameters ---------- bcast : bool if true the parameters are updated on all process """ if MCH.do_master_only_work(): if bcast: self.__update_parameters(bcast=False) param = self.get_parameters() else: param = self.mesh.get_parameters() else: param = None if bcast: param = MCH.master_broadcasts_array_to_all(param) if MCH.do_master_only_work() ^ bcast: for key in param: self.__dict__[key] = param[key]
[docs] def reset_parameters(self): """ reset model parameters """ self.L = 5000 self.y_c = 0 self.z_c = 0 self.Outer_D = 5 # mm self.Nerve_D = 250 # um self.fascicles = {} self.Perineurium_thickness = {} self.electrodes = {} self.i_stim = 1e-3 # A self.jstims = [] self.j_electrode = {} self.mesh_file_status = self.mesh_file is not None self.is_sim_ready = False self.is_meshed = False self.is_computed = False del self.mesh if not self.mesh_file_status: if MCH.do_master_only_work: self.mesh = NerveMshCreator( Length=self.L, Outer_D=self.Outer_D, Nerve_D=self.Nerve_D, y_c=self.y_c, z_c=self.z_c, ) else: self.mesh = None
##################### ## customize model ## #####################
[docs] def reshape_outerBox(self, Outer_D, res="default"): """ Reshape the size of the FEM simulation outer box Parameters ---------- outer_D : float FEM simulation outer box diameter, in mm, WARNING, this is the only parameter in mm ! """ if not self.mesh_file_status: self.mesh.reshape_outerBox(Outer_D=Outer_D, res=res) self.__update_parameters()
[docs] def reshape_nerve(self, Nerve_D=None, Length=None, y_c=0, z_c=0, res="default"): """ Reshape the nerve of the FEM simulation Parameters ---------- Nerve_D : float Nerve diameter, in um Length : float Nerve length, in um y_c : float Nerve center y-coordinate in um, 0 by default z_c : float Nerve z-coordinate center in um, 0 by default Perineurium_thickness :float Thickness of the Perineurium sheet surounding the fascicles in um, 5 by default """ if not self.mesh_file_status: self.L = Length or self.L self.Nerve_D = Nerve_D or self.Nerve_D self.mesh.reshape_nerve( Nerve_D=self.Nerve_D, Length=self.L, y_c=y_c, z_c=z_c, res=res ) self.__update_parameters()
[docs] def reshape_fascicle( self, Fascicle_D, y_c=0, z_c=0, ID=None, Perineurium_thickness=None, res="default", ): """ Reshape a fascicle of the FEM simulation Parameters ---------- Fascicle_D : float Fascicle diameter, in um y_c : float Fascicle center y-coodinate in um, 0 by default z_c : float Fascicle center y-coodinate in um, 0 by default ID : int If the simulation contains more than one fascicles, ID number of the fascicle to reshape as in FENICS """ if not self.mesh_file_status: self.mesh.reshape_fascicle(Fascicle_D, y_c, z_c, ID, res) self.__update_parameters() if ID is None: if self.Perineurium_thickness == {}: ID = 0 else: ## To check when not all ID are fixed ID = max(self.Perineurium_thickness) + 1 p_th = Perineurium_thickness or get_perineurial_thickness(Fascicle_D) self.Perineurium_thickness[ID] = p_th
[docs] def remove_fascicles(self, ID=None): """ remove a fascicle of the FEM simulation Parameters ---------- ID : int, None ID number of the fascicle to remove, if None, remove all fascicles, by default None """ if ID is None: self.fascicles = {} self.Perineurium_thickness = {} elif ID in self.fascicles: del self.fascicles[ID] del self.Perineurium_thickness[ID] self.mesh.remove_fascicles(ID=ID)
[docs] def add_electrode(self, elec_type, ID=None, res="default", **kwargs): """ TO BE WRITTEN """ if not self.mesh_file_status: self.mesh.add_electrode(elec_type=elec_type, ID=ID, res=res, **kwargs)
###################### ## setup simulation ## ######################
[docs] def setup_simulations(self, comm=None, rank=None): """ TO BE WRITTEN """ if comm is not None: self.comm = comm if rank is not None: self.rank = rank if not self.is_meshed: self.build_and_mesh() if not self.is_sim_ready and (MCH.do_master_only_work() or self.is_multi_proc): t0 = time.time() # SETTING DOMAINS # del self.mesh self.sim = FEMSimulation( mesh_file=self.mesh_file, mesh=self.mesh, elem=self.elem, comm=self.comm, rank=self.rank, ) # Outerbox domain self.__set_domains() # SETTING INTERNAL BOUNDARY CONDITION (for perineuriums) self.__set_iboundaries() # SETTING BOUNDARY CONDITION # Ground (to the external ring of Outerbox) self.sim.add_boundary(mesh_domain=1, btype="Dirichlet", value=0) # Injected current from active electrode # For EIT change in for E in elec_patren: for _, (E, active_elec) in enumerate(self.electrodes.items()): E_var = "E" + str(E) mesh_domain_3D = self.__find_elec_subdomain(active_elec) self.j_electrode[E_var] = 0 e_dom = ( ENT_DOM_offset["Surface"] + ENT_DOM_offset["Electrode"] + (2 * E) ) self.sim.add_boundary( mesh_domain=e_dom, btype="Neuman", variable=E_var, mesh_domain_3D=mesh_domain_3D, ) # set a parallelizable preconditionner if sim solve on multiple precesses if self.is_multi_proc: self.sim.set_solver_opt(pc_type="hypre") self.sim.setup_sim(**self.j_electrode) self.is_sim_ready = True self.setup_timer += time.time() - t0
[docs] def set_materials( self, Outer_mat=None, Epineurium_mat=None, Endoneurium_mat=None, Perineurium_mat=None, Electrodes_mat=None, ): """ Set material files for any domain Parameters ---------- Outer_mat :str Outer box material fname if None not changed, by default None Epineurium_mat :str Epineurium material fname if None not changed, by default None Endoneurium_mat :str Endoneurium material fname if None not changed, by default None Perineurium_mat :str Outer material fname if None not changed, by default None Electrodes_mat :str Electrodes material fname if None not changed, by default None """ self.Outer_mat = Outer_mat or self.Outer_mat self.Epineurium_mat = Epineurium_mat or self.Epineurium_mat self.Endoneurium_mat = Endoneurium_mat or self.Endoneurium_mat self.Perineurium_mat = Perineurium_mat or self.Perineurium_mat self.Electrodes_mat = Electrodes_mat or self.Electrodes_mat
def __set_domains(self): """ Internal use only: set the material properties of physical domains """ self.sim.add_domain( mesh_domain=ENT_DOM_offset["Outerbox"], mat_pty=self.Outer_mat ) # Nerve domain self.sim.add_domain( mesh_domain=ENT_DOM_offset["Nerve"], mat_pty=self.Epineurium_mat ) for i in self.fascicles.keys(): self.sim.add_domain( mesh_domain=ENT_DOM_offset["Fascicle"] + (2 * i), mat_pty=self.Endoneurium_mat, ) for _, (i, elec) in enumerate(self.electrodes.items()): if not elec["type"] == "LIFE": self.sim.add_domain( mesh_domain=ENT_DOM_offset["Electrode"] + (2 * i), mat_pty=self.Electrodes_mat, ) def __set_iboundaries(self): """ Internal use only: set internam boundaries """ for i in self.fascicles: thickness = self.Perineurium_thickness[i] f_dom = ENT_DOM_offset["Fascicle"] + (2 * i) self.sim.add_inboundary( mesh_domain=ENT_DOM_offset["Surface"] + f_dom, mat_pty=self.Perineurium_mat, thickness=thickness, in_domains=[f_dom], ) def __find_elec_subdomain(self, elec) -> int: """ Internal use only: """ if elec["type"] == "LIFE": y_e, z_e = elec["kwargs"]["y_c"], elec["kwargs"]["z_c"] for i in self.fascicles: fascicle = self.fascicles[i] d_f, y_f, z_f = fascicle["d"], fascicle["y_c"], fascicle["z_c"] if (y_e - y_f) ** 2 + (z_e - z_f) ** 2 < (d_f / 2) ** 2: return 10 + 2 * i return 0 def __find_elec_jstim(self, elec, I=None) -> float: """ Internal use only: return electrical current density in electrode from current valeu and electrode geometry """ # Unitary stimulation if I is not None: self.i_stim = I """if elec["type"] == "CUFF": d_e = self.Nerve_D l_e = elec["kwargs"]["contact_length"] elif elec["type"] == "LIFE": d_e = elec["kwargs"]["d"] l_e = elec["kwargs"]["length"]""" S = self.sim jstim = self.i_stim / S return jstim ################### ## Use the model ## ###################
[docs] def get_meshes(self): """ Visualize the model mesh WARNING: debbug use only + might not work on all os """ if self.is_meshed: self.mesh.visualize()
[docs] def build_and_mesh(self): """ Build the geometry and perform meshing process """ if not self.mesh_file_status and not self.is_meshed: t0 = time.time() self.__update_parameters() if self.N_fascicle == 0: self.reshape_fascicle( Fascicle_D=self.default_fascicle["d"], y_c=self.default_fascicle["y_c"], z_c=self.default_fascicle["z_c"], res=self.default_fascicle["res"], ) if self.N_electrode == 0: self.add_electrode( elec_type=self.default_electrode["elec_type"], x_c=self.default_electrode["x_c"], y_c=self.default_electrode["y_c"], z_c=self.default_electrode["z_c"], length=self.default_electrode["length"], d=self.default_electrode["d"], res=self.default_electrode["res"], ) self.__update_parameters(bcast=self.is_multi_proc) if MCH.do_master_only_work(): self.mesh.compute_mesh() if self.is_multi_proc: synchronize_processes() self.is_meshed = True self.mesh.get_info(verbose=True) self.meshing_timer += time.time() - t0
[docs] def solve(self, comm=None, rank=None): """ Solve the model """ if not self.is_sim_ready: self.setup_simulations() if not self.is_computed: t0 = time.time() # For EIT change in for E in elec_patren: for E in range(self.N_electrode): for i_elec in self.j_electrode: if i_elec == "E" + str(E): e_dom = ( ENT_DOM_offset["Surface"] + ENT_DOM_offset["Electrode"] + (2 * E) ) self.j_electrode[i_elec] = self.i_stim / self.sim.get_surface( e_dom ) else: self.j_electrode[i_elec] = 0 self.sim.setup_sim(**self.j_electrode) self.sim_res += [self.sim.solve()] self.sim_res[-1].vout.label = "E" + str(E) self.is_computed = True self.solving_timer += time.time() - t0
###################### ## results handling ## ######################
[docs] def get_potentials(self, x, y, z, E=-1): """ Get the potential on a line to get extracellular potential for axons stimulation. Parameters ---------- x : np.array array of x coordinates in the model y : float y-coordinate of the axon z : float z-coordinate of the axon E : int ID of the electrode from witch the potentials should be evaluated Returns: -------- array All potential for all paramtric sweeps (all electrodes in NRV2 models)\ (line: electrode selection, column: potential) """ if self.is_computed: t0 = time.time() line = [[x_, y, z] for x_ in x] if self.N_electrode == 1 or (E < self.N_electrode and E >= 0): potentials = self.sim_res[E].eval(line, self.is_multi_proc) * V else: potentials = [] for E in range(self.N_electrode): potentials += [self.sim_res[E].eval(line, self.is_multi_proc)] potentials = np.transpose(np.array(potentials) * V) self.access_res_timer += time.time() - t0 return potentials