"""
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._parameters import parameters
from .fenics_utils._FEMSimulation import FEMSimulation
from ._FENICS_model import FENICS_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]
class FENICS_lumped_impedance_model(FENICS_model):
"""
A class for FENICS Finite Element Models, inherits from FENICS_model.
"""
[docs]
def __init__(
self,
fname=None,
Ncore=None,
handle_server=False,
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__(
fname=fname,
Ncore=Ncore,
handle_server=handle_server,
elem=elem,
comm=comm,
rank=rank,
)
self.axons = {}
self.axons_gmem = {}
self.Myelin_mat = "endoneurium_ranck"
self.Axoplasmic_mat = 1 / 70
[docs]
def reshape_axon(
self,
d,
y=0,
z=0,
ID=None,
myelinated=False,
res="default",
res_node="default",
**kwargs
):
if not self.mesh_file_status:
self.mesh.reshape_axon(
d=d,
y=y,
z=z,
ID=ID,
myelinated=myelinated,
res=res,
res_node=res_node,
**kwargs
)
self.__update_parameters()
######################
## setup simulation ##
######################
def __set_domains(self):
super().__set_domains()
for i in self.axons.keys():
self.sim.add_domain(
mesh_domain=ENT_DOM_offset["Axons"] + (2 * i),
mat_pty=self.Endoneurium_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],
)
[docs]
def set_materials(
self,
Outer_mat=None,
Epineurium_mat=None,
Endoneurium_mat=None,
Perineurium_mat=None,
Electrodes_mat=None,
Myelin_mat=None,
Axoplasmic_mat=None,
):
super().set_materials(
Outer_mat, Epineurium_mat, Endoneurium_mat, Perineurium_mat, Electrodes_mat
)
self.Myelin_mat = Myelin_mat or self.Myelin_mat
self.Axoplasmic_mat = Axoplasmic_mat or self.Axoplasmic_mat
[docs]
def set_axon_membrane(id, gmem_pty):
pass
###################
## 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):
self.j_electrode[i_elec] = self.jstims[E]
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