Source code for nrv.fmod.FEM.mesh_creator._NerveMshCreator

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

from cmath import phase
import numpy as np
from math import log10, floor

from ....backend._NRV_Class import NRV_class, load_any
from ....backend._log_interface import rise_error, rise_warning
from ....backend._MCore import MCH

# from ....nmod.myelinated import myelinated
from ....utils._units import mm
from ....backend._file_handler import rmv_ext
from ._MshCreator import MshCreator, pi

ENT_DOM_offset = {
    "Volume": 0,
    "Surface": 1,
    "Outerbox": 0,
    "Nerve": 2,
    "Fascicle": 10,
    "Electrode": 100,
    "Axon": 1000,
}

ELEC_TYPES = ["CUFF MP", "CUFF", "LIFE"]


[docs] def is_NerveMshCreator(object): """ check if an object is a NerveMshCreator, return True if yes, else False Parameters ---------- result : object object to test Returns ------- bool True it the type is a NerveMshCreator object """ return isinstance(object, NerveMshCreator)
[docs] class NerveMshCreator(MshCreator): """ Class allowing to generate Nerve shape 3D gmsh mesh with labeled physical domain Contains methodes dealing with the mesh geometries, physical domains and feilds Inherit from MshCreator class. see MshCreator for further detail """
[docs] def __init__( self, Length=10000, Outer_D=5, Nerve_D=4000, y_c=0, z_c=0, ver_level=2, ): """ initialisation of the NerveMshCreator Parameters ---------- Length : float Nerve length in um, by default 10000 Outer_D : float Outer box diameter in mm, by default 5 Nerve_D : float Nerve diameter in um, by default 4000 y_c : float y-axis position of the Nerve center, by default 0 z_c : float z-axis position of the Nerve center, by default 0 ver_level : int(0,1,2,3,4,5,99) verbosity level of gmsh (see MshCreator.set_verbosity), by default 2 """ super().__init__(D=3, ver_level=ver_level) self.L = Length self.y_c = y_c self.z_c = z_c self.surf_c = [self.L / 2] self.gnd_facet = {"outfacet": True, "lfacet": False, "rfacet": False} if Outer_D is None: self.is_outer = False self.gnd_facet["outfacet"] = False self.gnd_facet["rfacet"] = True self.Outer_D = 0 self.Outer_entities = {"face": [], "volume": []} else: self.is_outer = True self.Outer_D = Outer_D * mm self.Outer_entities = {"face": [], "volume": []} self.Nerve_D = Nerve_D self.Nerve_entities = {"face": [], "volume": []} self.N_fascicle = 0 self.fascicles = {} self.N_axon = 0 self.axons = {} self.N_electrode = 0 self.electrodes = {} self.default_res = { "Outerbox": self.Outer_D / 5, "Outerbox_tresholded": False, "Nerve": self.Nerve_D / 5, "Fascicle": 500, "Axon": 10, "Electrode": 1000, } self.Ox = None # add an Ox line (need for thresholded resolution) self.res = self.default_res["Outerbox"] self.alpha_Outerboxres = 0.05 self.is_geo = False self.is_dom = False self.is_refined = False
[docs] def get_parameters(self): param = {} param["res"] = self.res param["L"] = self.L param["y_c"] = self.y_c param["z_c"] = self.z_c param["surf_c"] = self.surf_c param["Outer_D"] = self.Outer_D param["Outer_entities"] = self.Outer_entities param["Nerve_D"] = self.Nerve_D param["Nerve_entities"] = self.Nerve_entities param["N_fascicle"] = self.N_fascicle param["fascicles"] = self.fascicles param["N_axon"] = self.N_axon param["axons"] = self.axons param["N_electrode"] = self.N_electrode param["electrodes"] = self.electrodes return param
[docs] def save( self, fname="nervemshcreator.json", save=True, mshfname=None, blacklist=[], **kwargs, ): """ Return extracellular context as dictionary and eventually save it as json file NB: caution, first argument is fname to match with MshCreator.save arguments Parameters ---------- save : bool if True, save in json files fname : str Path and Name of the saving file, by default "extracel_context.json" Returns ------- context_dic : dict dictionary containing all information """ blacklist += [ "model", ] if self.is_generated and save: if mshfname is None: mshfname = rmv_ext(fname) + ".msh" super().save(fname=mshfname, generate=False) fname = rmv_ext(fname) + ".json" return NRV_class.save( self, save=save, fname=fname, blacklist=blacklist, **kwargs )
[docs] def load(self, data, mshfname=None, blacklist={}, **kwargs): """ Generic loading method for NRV class instance Parameters ---------- data : dict Dictionary containing the NRV object blacklist : dict, optional Dictionary containing the keys to be excluded from the load **kwargs : dict, optional Additional arguments to be passed to the load method of the NRV object """ NRV_class.load(self, data, **kwargs) self.electrodes = {int(k): v for k, v in self.electrodes.items()} self.fascicles = {int(k): v for k, v in self.fascicles.items()} self.axons = {int(k): v for k, v in self.axons.items()} if self.is_generated and self.file != "": super().load(self.file)
[docs] def compute_mesh(self, master_only=True): """ Compute mesh geometry, domains and resolution and then generate the mesh """ if MCH.do_master_only_work() or not master_only: self.compute_geo() self.compute_domains() self.compute_res() self.generate()
[docs] def set_gnd_facet(self, outfacet=None, lfacet=None, rfacet=None): """ Set which of the outer facet should be the ground (element 1) Parameters ---------- outfacet : bool or None if true outer ring facet is included to the element 1, if None not modified by default None lfacet : bool or None if true left facet is included to the element 1, if None not modified by default None rfacet : bool or None if true right facet is included to the element 1, if None not modified by default None """ if outfacet is not None: self.gnd_facet["outfacet"] = outfacet if lfacet is not None: self.gnd_facet["lfacet"] = lfacet if rfacet is not None: self.gnd_facet["rfacet"] = rfacet
############################################################################################### ##################################### geometry definition ################################## ###############################################################################################
[docs] def compute_geo(self): """ Compute the mesh geometry """ if not self.is_geo: if self.is_outer: self.add_cylinder(0, self.y_c, self.z_c, self.L, self.Outer_D / 2) if self.default_res["Outerbox_tresholded"]: self.Ox = self.add_line((0, 0, 0), (self.L, 0, 0)) self.add_cylinder(0, self.y_c, self.z_c, self.L, self.Nerve_D / 2) for fascicle in self.fascicles.values(): self.add_cylinder( 0, fascicle["y_c"], fascicle["z_c"], self.L, fascicle["d"] / 2 ) for i in self.axons: self.add_axon(i) for i in self.electrodes: electrode = self.electrodes[i] if "CUFF MP" in electrode["type"]: self.add_CUFF_MP(ID=i, **electrode["kwargs"]) elif "CUFF MEA" in electrode["type"]: self.add_CUFF_MEA(ID=i, **electrode["kwargs"]) elif "CUFF" in electrode["type"]: self.add_CUFF(ID=i, **electrode["kwargs"]) elif "LIFE" in electrode["type"]: self.add_LIFE(ID=i, **electrode["kwargs"]) self.fragment(verbose=False) self.is_geo = True
[docs] def reshape_outerBox(self, Outer_D=None, res="default", tresholded_res=None): """ 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 Outer_D is not None: """ ## See how if it should be added if self.default_res["Outerbox"] == self.Outer_D/5: self.default_res["Outerbox"] = Outer_D/5 """ self.Outer_D = Outer_D * mm if tresholded_res is not None: self.default_res["Outerbox_tresholded"] = tresholded_res if not res == "default": self.default_res["Outerbox"] = res if self.default_res["Outerbox"] > self.Outer_D / 3: self.default_res["Outerbox"] = self.Outer_D / 3
[docs] def reshape_nerve( self, Nerve_D=None, Length=None, y_c=None, z_c=None, 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 """ if Length is not None: self.L = Length if Nerve_D is not None: if self.default_res["Nerve"] == self.Nerve_D / 5: self.default_res["Nerve"] = Nerve_D / 5 self.Nerve_D = Nerve_D if y_c is not None: self.y_c = y_c if z_c is not None: self.z_c = z_c if not res == "default": self.default_res["Nerve"] = res if self.default_res["Nerve"] > self.Nerve_D / 5: self.default_res["Nerve"] = self.Nerve_D / 5
[docs] def reshape_fascicle(self, d, y_c=0, z_c=0, ID=None, res="default"): """ Reshape a fascicle of the FEM simulation Parameters ---------- 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 COMSOL """ if ID not in self.fascicles: if ID is None: ID = 0 while ID in self.fascicles: ID += 1 self.N_fascicle += 1 if res == "default": res = self.default_res["Fascicle"] if d / 5 < res: res = d / 5 self.fascicles[ID] = { "y_c": y_c, "z_c": z_c, "d": d, "res": res, "face": None, "volume": None, }
[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.N_fascicle = 0 elif ID in self.fascicles: del self.fascicles[ID] self.N_fascicle -= 1
[docs] def reshape_axon( self, d, y=0, z=0, ID=None, myelinated=False, res="default", res_node="default", **kwargs, ): """ Reshape a axon of the FEM simulation Parameters ---------- 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 COMSOL """ if ID not in self.axons: if ID is None: ID = 0 while ID in self.axons: ID += 1 self.N_axon += 1 if res == "default": res = self.default_res["Axon"] if d / 3 < res: res = d / 3 if res_node == "default": res_node = res / 3 self.axons[ID] = kwargs self.axons[ID].update( { "y": y, "z": z, "d": d, "myelinated": myelinated, "res": res, "res_node": res_node, "face": [], "volume": [], } )
[docs] def add_electrode(self, elec_type, ID=None, res="default", **kwargs): """ """ if elec_type not in ELEC_TYPES + ["CUFF MEA"]: rise_warning( elec_type + " not implemented, electrode will not be added\nList of Electrode types: " + str(ELEC_TYPES) ) if ID not in self.electrodes: if ID is None: ID = 0 while ID in self.electrodes: if ( "CUFF MEA" in self.electrodes[ID]["type"] or "CUFF MP" in self.electrodes[ID]["type"] ): ID += self.electrodes[ID]["kwargs"]["N"] else: ID += 1 self.N_electrode += 1 if res == "default": res = self.default_res["Electrode"] if "LIFE" in elec_type: if "d" in kwargs: d = kwargs["d"] else: d = 25 if d / 3 < res: res = d / 3 self.electrodes[ID] = {"type": elec_type, "res": res, "kwargs": kwargs}
#################################################################################################### ##################################### domains definition ######################################## ####################################################################################################
[docs] def compute_domains(self): if not self.is_geo: rise_error("compute geometry before domain") elif not self.is_dom: self.__link_entity_domains(dim=2) self.__link_entity_domains(dim=3) self.compute_entity_domain() self.is_dom = True
def __is_outerbox(self, dx, dy, dz, com, dim_key): """ Internal use only: check if volume is the box or face is external face of box Parameters ---------- dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox """ outer_ring_test = ( np.allclose([dx, dy, dz], [self.L, self.Outer_D, self.Outer_D]) and self.is_outer ) if dim_key == "volume": return outer_ring_test outer_ring_test = outer_ring_test and self.gnd_facet["outfacet"] left_facet_test = np.allclose(com[0], [0]) and self.gnd_facet["lfacet"] rigth_facet_test = np.allclose(com[0], [self.L]) and self.gnd_facet["rfacet"] return outer_ring_test or left_facet_test or rigth_facet_test def __is_nerve(self, dx, dy, dz): """ Internal use only: check if volume is nerve or face is external face of nerve Parameters ---------- dx : float length along x of the tested entity dy : float length along y of the tested entity dz : float length along z of the tested entity """ return np.allclose([dy, dz], [self.Nerve_D, self.Nerve_D]) and not np.isclose( dx, 0 ) def __is_fascicle(self, ID, dx, dy, dz, com, dim_key): """ Internal use only: check if volume is fascicle ID or face is external face of fascicle ID Parameters ---------- ID :int ID of fascicle to test dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass dim_key :"face" or "volume" element type """ fascicle = self.fascicles[ID] # Already computed status_test = fascicle[dim_key] is None # test good diameter size_test = np.allclose([dx, dy, dz], [self.L, fascicle["d"], fascicle["d"]]) # test center of mass in fascicle com_test = np.allclose( com, (self.L / 2, fascicle["y_c"], fascicle["z_c"]), rtol=1, atol=fascicle["d"] / 2, ) return status_test and size_test and com_test def __is_axon_node(self, ID, dx, dy, dz, com, dim_key): """ Internal use only: check if volume is axon ID or face is external face of axon ID Parameters ---------- ID :int ID of axon to test dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass dim_key :"face" or "volume" element type """ axon = self.axons[ID] # Already computed if not axon["myelinated"]: return False # test good diameter size_test = np.allclose( [dx, dy, dz], [axon["node_l"], axon["node_d"], axon["node_d"]] ) # test first and last node if "first_node_l" in axon and not size_test: size_test = np.allclose( [dx, dy, dz], [axon["first_node_l"], axon["node_d"], axon["node_d"]] ) if "last_node_l" in axon and not size_test: size_test = np.allclose( [dx, dy, dz], [axon["last_node_l"], axon["node_d"], axon["node_d"]] ) # test center of mass in axon com_test = np.allclose(com[1:], (axon["y"], axon["z"])) return size_test and com_test def __is_axon(self, ID, dx, dy, dz, com, dim_key): """ Internal use only: check if volume is axon ID or face is external face of axon ID Parameters ---------- ID :int ID of axon to test dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass dim_key :"face" or "volume" element type """ axon = self.axons[ID] # Already computed status_test = True # axon[dim_key] is None # test good diameter size_test = np.allclose([dy, dz], [axon["d"], axon["d"]]) size_test &= not np.isclose(dx, 0) # test center of mass in axon com_test = np.allclose(com[1:], (axon["y"], axon["z"]), atol=axon["d"] / 2) return status_test and size_test and com_test def __is_CUFF_MP_electrode(self, ID, rc, dx, teta, com, dim_key): """ Internal use only: check if volume is electrode ID or face is external face of fascicle ID Parameters ---------- dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass """ elec_kwargs = self.electrodes[ID]["kwargs"] if self.electrodes[ID]["type"] != "CUFF MP": return False # test good length size_test = np.allclose([dx], [elec_kwargs["contact_length"]]) # test good x-location com_test = np.isclose(com[0], elec_kwargs["x_c"]) # test good polar coordinate if elec_kwargs["is_volume"] or dim_key == "volume": N = elec_kwargs["N"] tetas = [teta for i in range(N)] angle_test = np.isclose(tetas, self.electrodes[ID]["angles"]).any() radius_test = rc > self.Nerve_D / 2 polar_test = angle_test and radius_test else: polar_test = rc <= self.Nerve_D / 2 return size_test and com_test and polar_test def __is_CUFF_electrode(self, ID, dx, dy, dz, com, dim_key): """ Internal use only: check if volume is electrode ID or face is external face of fascicle ID Parameters ---------- dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass """ elec_kwargs = self.electrodes[ID]["kwargs"] if self.electrodes[ID]["type"] != "CUFF": return False # test good diameter Syz = self.Nerve_D if elec_kwargs["is_volume"] or dim_key == "volume": Syz += 2 * elec_kwargs["contact_thickness"] size_test = np.allclose([dx, dy, dz], [elec_kwargs["contact_length"], Syz, Syz]) # test center of mass in CUFF com_test = np.allclose(com, (elec_kwargs["x_c"], self.y_c, self.z_c)) return size_test and com_test def __is_LIFE_electrode(self, ID, dx, dy, dz, com): """ Internal use only: check if volume is electrode ID or face is external face of fascicle ID Parameters ---------- dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass """ elec_kwargs = self.electrodes[ID]["kwargs"] if self.electrodes[ID]["type"] != "LIFE": return False # test good diameter size_test = np.allclose( [dx, dy, dz], [elec_kwargs["length"], elec_kwargs["d"], elec_kwargs["d"]] ) # test center of mass in LIFE com_test = np.allclose( com, (elec_kwargs["x_c"], elec_kwargs["y_c"], elec_kwargs["z_c"]), atol=elec_kwargs["d"] / 2, ) return size_test and com_test def __link_entity_domains(self, dim): """ Internal use only: link all entities from Parameters ---------- dx : float length along x of the entity boundbox dy : float length along y of the entity boundbox dz : float length along z of the entity boundbox com : tupple(float) entity Center of Mass """ if dim == 2: entities, ent_com, ent_bd = self.get_faces(com=True, bd=True) key = "face" offset = 1 else: entities, ent_com, ent_bd = self.get_volumes(com=True, bd=True) key = "volume" offset = 0 for i in range(len(entities)): bd_x = abs(ent_bd[i][3] - ent_bd[i][0]) bd_y = abs(ent_bd[i][4] - ent_bd[i][1]) bd_z = abs(ent_bd[i][5] - ent_bd[i][2]) if self.__is_outerbox(bd_x, bd_y, bd_z, ent_com[i], key): self.Outer_entities[key] += [entities[i][1]] elif self.__is_nerve(bd_x, bd_y, bd_z): self.Nerve_entities[key] += [entities[i][1]] for j in self.fascicles: if self.__is_fascicle(j, bd_x, bd_y, bd_z, ent_com[i], key): self.fascicles[j][key] = entities[i][1] for j in self.axons: if self.__is_axon_node(j, bd_x, bd_y, bd_z, ent_com[i], key): node_id = int(ent_com[i][0] // self.axons[j]["deltax"]) self.axons[j]["nodes_" + key][node_id] = [entities[i][1]] elif self.__is_axon(j, bd_x, bd_y, bd_z, ent_com[i], key): self.axons[j][key] += [entities[i][1]] for j in self.electrodes: r_c = ( (ent_com[i][1] - self.y_c) ** 2 + (ent_com[i][2] - self.z_c) ** 2 ) ** 0.5 teta = phase( complex(ent_com[i][2] - self.z_c, ent_com[i][1] - self.y_c) ) % (2 * pi) if self.__is_CUFF_MP_electrode(j, r_c, bd_x, teta, ent_com[i], key): N = self.electrodes[j]["kwargs"]["N"] ID_EA = round(teta * N / (2 * pi)) % N if dim == 3: if self.electrodes[j][key][ID_EA] is None: self.electrodes[j][key][ID_EA] = entities[i][1] else: rise_warning( "two volumes can be same electrode only the first one is kept" ) else: if self.electrodes[j][key][ID_EA] is None: self.electrodes[j][key][ID_EA] = entities[i][1] else: self.electrodes[j][key][ID_EA] = [ entities[i][1], self.electrodes[j][key][ID_EA], ] elif self.__is_CUFF_electrode(j, bd_x, bd_y, bd_z, ent_com[i], key): self.electrodes[j][key] = entities[i][1] if not self.electrodes[j]["kwargs"]["is_volume"] and key == "face": self.Nerve_entities[key].remove(entities[i][1]) elif self.__is_LIFE_electrode(j, bd_x, bd_y, bd_z, ent_com[i]): self.electrodes[j][key] = entities[i][1]
[docs] def compute_entity_domain(self): """ link all geometrical entities to corresponding physical domains """ # Outer box domain if self.is_outer: self.add_domains(obj_IDs=self.Outer_entities["volume"], phys_ID=0, dim=3) self.add_domains(obj_IDs=self.Outer_entities["face"], phys_ID=1, dim=2) else: self.add_domains(obj_IDs=self.Outer_entities["face"], phys_ID=1, dim=2) # nerve domain self.add_domains(obj_IDs=self.Nerve_entities["volume"], phys_ID=2, dim=3) self.add_domains(obj_IDs=self.Nerve_entities["face"], phys_ID=3, dim=2) for j in self.fascicles: id_ph = ENT_DOM_offset["Fascicle"] + (2 * j) if id_ph < ENT_DOM_offset["Electrode"]: self.add_domains( obj_IDs=self.fascicles[j]["volume"], phys_ID=id_ph, dim=3 ) id_ph += ENT_DOM_offset["Surface"] self.add_domains( obj_IDs=self.fascicles[j]["face"], phys_ID=id_ph, dim=2 ) else: rise_warning("Too much Fascicles: " + str(j) + " not added") for j, ax in self.axons.items(): id_ph = ENT_DOM_offset["Axon"] + (2 * j) self.add_domains(obj_IDs=ax["volume"], phys_ID=id_ph, dim=3) id_ph += ENT_DOM_offset["Surface"] self.add_domains(obj_IDs=ax["face"], phys_ID=id_ph, dim=2) # adding one physical domain by node if ax["myelinated"]: for i_node in range(len(ax["nodes_volume"])): node_id = 2 * i_node n = floor(log10(node_id + 1) + 1) ax_id = (ENT_DOM_offset["Axon"] + (2 * j)) * 10**n id_ph = node_id + ax_id self.add_domains( obj_IDs=ax["nodes_volume"][i_node], phys_ID=id_ph, dim=3 ) id_ph += ENT_DOM_offset["Surface"] self.add_domains( obj_IDs=ax["nodes_face"][i_node], phys_ID=id_ph, dim=2 ) for j in self.electrodes: id_ph = ENT_DOM_offset["Electrode"] + (2 * j) if id_ph < ENT_DOM_offset["Axon"]: if "CUFF MP" in self.electrodes[j]["type"]: for ID_EA in range(self.electrodes[j]["kwargs"]["N"]): if self.electrodes[j]["kwargs"]["is_volume"]: self.add_domains( obj_IDs=self.electrodes[j]["volume"][ID_EA], phys_ID=id_ph, dim=3, ) id_ph += ENT_DOM_offset["Surface"] self.add_domains( obj_IDs=self.electrodes[j]["face"][ID_EA], phys_ID=id_ph, dim=2, ) id_ph += 1 else: if self.electrodes[j]["kwargs"]["is_volume"]: self.add_domains( obj_IDs=self.electrodes[j]["volume"], phys_ID=id_ph, dim=3 ) id_ph += ENT_DOM_offset["Surface"] self.add_domains( obj_IDs=self.electrodes[j]["face"], phys_ID=id_ph, dim=2 ) else: rise_warning("Too much Electrodes: " + str(j) + " not added")
#################################################################################################### ################################### field (res) definition ###################################### ####################################################################################################
[docs] def compute_res(self): if not self.is_dom and self.is_geo: rise_error("compute geometry and domain before resolution") elif not self.is_refined: self.res = max(self.default_res.values()) fields = [] # Outerbox field if self.is_outer: fields += self.__refine_Outer_box(alpha=self.alpha_Outerboxres) # Nerve field fields += [ self.refine_entities( ent_ID=self.Nerve_entities["volume"], res_in=self.default_res["Nerve"], dim=3, res_out=None, IncludeBoundary=True, ) ] # Facsicle fields for fascicle in self.fascicles.values(): fields += [ self.refine_entities( ent_ID=fascicle["volume"], res_in=fascicle["res"], dim=3, res_out=None, IncludeBoundary=True, ) ] # Axon fields for axon in self.axons.values(): fields += self.__refine_axon(axon) # Electrodes fields for electrode in self.electrodes.values(): if electrode["type"] in ["CUFF MEA", "CUFF MP"]: for ID_EA in range(electrode["kwargs"]["N"]): fields += [ self.refine_entities( ent_ID=electrode["volume"][ID_EA], res_in=electrode["res"], dim=3, res_out=None, IncludeBoundary=True, ) ] else: fields += [ self.refine_entities( ent_ID=electrode["volume"], res_in=electrode["res"], dim=3, res_out=None, IncludeBoundary=True, ) ] self.refine_min(fields) self.is_refined = True
def __refine_axon(self, axon): """ """ # First set the field for unmyelinated axon or myelin fields = [ self.refine_entities( ent_ID=axon["volume"], res_in=axon["res"], dim=3, res_out=None, IncludeBoundary=True, ) ] # Then set the field for myelinated axon nodes if axon["myelinated"]: for i in range(len(axon["nodes_volume"])): fields += [ self.refine_entities( ent_ID=axon["nodes_volume"][i], res_in=axon["res_node"], dim=3, res_out=None, IncludeBoundary=True, ) ] return fields def __refine_Outer_box(self, alpha=0.1): """ """ if not self.default_res["Outerbox_tresholded"]: field = [ self.refine_entities( ent_ID=self.Outer_entities["volume"], res_in=self.default_res["Outerbox"], dim=3, res_out=None, IncludeBoundary=True, ) ] else: dmin = self.Nerve_D dmax = dmin + alpha * self.Outer_D field = [ self.refine_threshold( ent_ID=self.Ox, dim=1, res_min=self.default_res["Nerve"], res_max=self.default_res["Outerbox"], dist_min=self.Nerve_D, dist_max=dmax, ) ] return field ################################################################ ############## complex volumes adding methods ################## ################################################################
[docs] def add_axon(self, ID): """ Add an axon to the mesh. Note ---- - if the axon is unmyelinated this method add only a cylinder to the mesh Parameters ---------- """ axon = self.axons[ID] res_min = axon["res_node"] if not axon["myelinated"]: self.add_cylinder(0, axon["y"], axon["z"], self.L, axon["d"] / 2) else: axon.update( { "L": self.L, "rec": "all", "N_seg_per_sec": True, "__NRVObject__": True, "nrv_type": "myelinated", } ) # Ideally it would be easier to use: ax = myelinated(**axon) # But myelinated cannot be imported without causing an import loop ax = load_any(data=axon) x = ax.first_section_size init_l = 0 n_nodes = 0 for sec in ax.axon_path_type[1:-1]: if sec == "node": r_sec = ax.nodeD / 2 l_sec = ax.nodelength n_nodes += 1 else: r_sec = ax.d / 2 if sec == "MYSA": l_sec = ax.paralength1 elif sec == "FLUT": l_sec = ax.paralength2 elif sec == "STIN": l_sec = ax.interlength # mrege the first (last) section with the next (previous) when it is too small if abs(x - l_sec) < res_min: init_l += l_sec ax.axon_path_type.pop(0) x += l_sec elif abs(self.L - (l_sec + x)) < res_min: ax.axon_path_type.pop(-1) else: self.add_cylinder(x, ax.y, ax.z, l_sec, r_sec) x += l_sec # adding first section if ax.axon_path_type[0] == "node": r_sec = ax.nodeD / 2 n_nodes += 1 self.axons[ID]["first_node_l"] = ax.first_section_size + init_l else: r_sec = ax.d / 2 self.add_cylinder(0, ax.y, ax.z, ax.first_section_size + init_l, r_sec) # adding first section if ax.axon_path_type[-1] == "node": r_sec = ax.nodeD / 2 n_nodes += 1 self.axons[ID]["last_node_l"] = self.L - x else: r_sec = ax.d / 2 self.add_cylinder(x, ax.y, ax.z, self.L - x, r_sec) self.axons[ID]["node_d"] = ax.nodeD self.axons[ID]["node_l"] = ax.nodelength self.axons[ID]["deltax"] = ax.deltax self.axons[ID]["n_nodes"] = n_nodes self.axons[ID]["nodes_volume"] = [None for _ in range(n_nodes)] self.axons[ID]["nodes_face"] = [None for _ in range(n_nodes)] del ax
[docs] def add_LIFE( self, ID=None, x_c=0, y_c=0, z_c=0, length=1000, d=25, is_volume=False ): """ Add LIFE electrode to the mesh Parameters ---------- ID :int electrod ID, ,by defalt None x_c :float x-position of the LIFE center in um, by default 0 y_c :float y-position of the LIFE center in um, by default 0 z_c :float z-position of the LIFE center in um, by default 0 length :float length of the LIFE electrod in um, by default 1000 d :float diameter of the LIFE electrod in um, by default 25 is_volume : bool if True in cylinder of the LIFE is kept on the mesh """ if ID is not None: self.electrodes[ID]["volume"] = [] self.electrodes[ID]["face"] = [] self.electrodes[ID]["kwargs"]["is_volume"] = is_volume x_active = x_c - length / 2 y_active = y_c z_active = z_c self.add_cylinder(x_active, y_active, z_active, length, d / 2)
[docs] def add_CUFF( self, ID=None, x_c=0, contact_length=100, is_volume=True, contact_thickness=None, insulator=True, insulator_thickness=None, insulator_length=None, insulator_offset=0, ): """ Add CUFF electrode to the mesh Parameters ---------- ID :int if not none and ID exist, change electrod ID,by defalt None x_c :float x-position of the CUFF center in um, by default 0 length of the CUFF electrod in um, by default 100 contact_length :float length along x of the contact site in um, by default 100 is_volume : bool if True the contact is kept on the mesh as a volume, by default True contact_thickness :float thickness of the contact site in um, by default 5 insulator :bool remove insulator ring from the mesh (no conductivity), by default True insulator_thickness :float thickness of the insulator ring in um, by default 20 insulator_length :float length along x of the insulator ring in um, by default 1000 """ if contact_thickness is None: contact_thickness = 0.1 * (self.Outer_D - self.Nerve_D) / 2 if contact_thickness < 0: contact_thickness = self.Nerve_D * 0.1 if insulator_thickness is None: insulator_thickness = min( 5 * contact_thickness, 0.4 * (self.Outer_D - self.Nerve_D) / 2 ) if insulator_length is None: insulator_length = 2 * contact_length # self.reshape_outerBox(tresholded_res=True) if ID is not None: self.electrodes[ID]["kwargs"]["contact_length"] = contact_length self.electrodes[ID]["kwargs"]["contact_thickness"] = contact_thickness self.electrodes[ID]["kwargs"]["is_volume"] = is_volume self.electrodes[ID]["kwargs"]["insulator_offset"] = insulator_offset self.electrodes[ID]["volume"] = [] self.electrodes[ID]["face"] = [] if ( self.electrodes[ID]["res"] > min(contact_length, insulator_thickness) / 2 ): self.electrodes[ID]["res"] = ( min(contact_length, insulator_thickness) / 2 ) x_active = x_c - contact_length / 2 y_active = self.y_c z_active = self.z_c cyl_act = self.add_cylinder( x_active, y_active, z_active, contact_length, self.Nerve_D / 2 + contact_thickness, ) cyl_ner2 = self.add_cylinder( x_active, y_active, z_active, contact_length, self.Nerve_D / 2 ) self.model.occ.cut([(3, cyl_act)], [(3, cyl_ner2)]) if insulator: x_insulator = x_c + insulator_offset - insulator_length / 2 y_insulator = self.y_c z_insulator = self.z_c cyl_ina = self.add_cylinder( x_insulator, y_insulator, z_insulator, insulator_length, self.Nerve_D / 2 + insulator_thickness, ) cyl_ner = self.add_cylinder( x_insulator, y_insulator, z_insulator, insulator_length, self.Nerve_D / 2, ) self.model.occ.cut([(3, cyl_ina)], [(3, cyl_ner), (3, cyl_act)])
[docs] def add_CUFF_MP( self, ID=None, N=4, x_c=0, contact_width=None, contact_length=100, is_volume=True, contact_thickness=None, insulator=True, insulator_thickness=None, insulator_length=None, insulator_offset=0, ): """ Add MultiPolar CUFF electrodes to the mesh Parameters ---------- ID :int if not none and ID exist, change electrod ID,by defalt None N :int Number of active site, by default 4 x_c :float x-position of the CUFF center in um, by default 0 length of the CUFF electrod in um, by default 100 contact_width :float or None width of the active sites around the nerve, if None set to cover 4/5 of the perimeter with active sites contact_length :float length along x of the contact site in um, by default 100 is_volume : bool if True the contact is kept on the mesh as a volume, by default True contact_thickness :float thickness of the contact site in um, by default 5 insulator :bool remove insulator ring from the mesh (no conductivity), by default True insulator_thickness :float thickness of the insulator ring in um, by default 20 insulator_length :float length along x of the insulator ring in um, by default 1000 """ if contact_width is None: contact_width = pi * self.Nerve_D * 4 / (5 * N) if contact_thickness is None: contact_thickness = 0.1 * (self.Outer_D - self.Nerve_D) / 2 if contact_thickness < 0: contact_thickness = self.Nerve_D * 0.1 if insulator_thickness is None: insulator_thickness = min( 5 * contact_thickness, 0.4 * (self.Outer_D - self.Nerve_D) / 2 ) if insulator_length is None: insulator_length = 2 * contact_length if ID is not None: self.electrodes[ID]["kwargs"]["contact_width"] = contact_width self.electrodes[ID]["kwargs"]["contact_length"] = contact_length self.electrodes[ID]["kwargs"]["contact_thickness"] = contact_thickness self.electrodes[ID]["kwargs"]["insulator_thickness"] = insulator_thickness self.electrodes[ID]["kwargs"]["insulator_length"] = insulator_length self.electrodes[ID]["kwargs"]["insulator_offset"] = insulator_offset self.electrodes[ID]["kwargs"]["is_volume"] = is_volume self.electrodes[ID]["volume"] = [None for i in range(N)] self.electrodes[ID]["face"] = [None for i in range(N)] if contact_length / 3 < self.electrodes[ID]["res"]: self.electrodes[ID]["res"] = contact_length / 2 angles = [] bar_fus = [] x_active = x_c - contact_length / 2 y_active = self.y_c - contact_width / 2 y_c = self.y_c z_c = self.z_c if insulator: if insulator_thickness is None: insulator_thickness = contact_thickness * 3 if insulator_length is None: insulator_length = contact_length * 2 x_insulator = x_c + insulator_offset - insulator_length / 2 cyl1 = self.add_cylinder( x_insulator, y_c, z_c, insulator_length, self.Nerve_D / 2 + insulator_thickness, ) cyl2 = self.add_cylinder( x_insulator, y_c, z_c, insulator_length, self.Nerve_D / 2 ) self.model.occ.cut([(3, cyl1)], [(3, cyl2)]) for i in range(N): bar = self.add_box( x_active, y_active, z_c, contact_length, contact_width, contact_thickness + self.Nerve_D / 2, ) angles += [(2 * pi * i) / (N)] self.rotate(bar, angles[-1], x_c, y_c, z_c, ax=1) bar_fus += [(3, bar)] if ID is not None: self.electrodes[ID]["angles"] = angles cyl = self.add_cylinder(0, self.y_c, self.z_c, self.L, self.Nerve_D / 2) self.model.occ.cut(bar_fus, [(3, cyl)])
[docs] def add_CUFF_MEA( self, ID=None, N=4, x_c=0, y_c=0, z_c=0, size=None, thickness=100, inactive=True, inactive_th=None, inactive_L=None, ): """ """ rise_warning("add_CUFF_MEA not updated use add_CUFF_MP instead") if size is None: s = pi * self.Nerve_D * 4 / (5 * N) size = (s, s) elif not np.iterable(size): size = (size, size) elif len(size) != 2: size = (size[0], size[0]) self.add_CUFF_MP( ID=ID, N=N, x_c=x_c, contact_width=size[1], contact_length=size[0], is_volume=True, contact_thickness=thickness, insulator=inactive, insulator_thickness=inactive_th, insulator_length=inactive_L, ) if ID is not None: self.electrodes[ID]["type"] = "CUFF MP"