####################################################################################
#
# STEPS - STochastic Engine for Pathway Simulation
# Copyright (C) 2007-2023 Okinawa Institute of Science and Technology, Japan.
# Copyright (C) 2003-2006 University of Antwerp, Belgium.
#
# See the file AUTHORS for details.
# This file is part of STEPS.
#
# STEPS is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License version 3,
# as published by the Free Software Foundation.
#
# STEPS is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
#################################################################################
###
import contextlib
import copy
import functools
import numbers
import operator
import pickle
import types
import warnings
import numpy
from steps import stepslib
import steps.API_1.utilities.meshio as smeshio
import steps.API_1.utilities.geom_decompose as sgdecomp
import steps.API_1.utilities.metis_support as smetis
import steps.API_1.utilities.morph_support as smorph
from . import utils as nutils
from . import model as nmodel
__all__ = [
'Geometry',
'Compartment',
'Patch',
'Membrane',
'EndocyticZone',
'ROI',
'TetMesh',
'DistMesh',
'DiffBoundary',
'SDiffBoundary',
'MeshPartition',
'LinearMeshPartition',
'MetisPartition',
'GmshPartition',
'MorphPartition',
'Reference',
'TetReference',
'TriReference',
'BarReference',
'VertReference',
'RefList',
'TetList',
'TriList',
'BarList',
'VertList',
'Point',
'BoundingBox',
'Morph',
]
###################################################################################################
ELEM_VERTEX = stepslib._py_ElementType.ELEM_VERTEX
ELEM_TRI = stepslib._py_ElementType.ELEM_TRI
ELEM_TET = stepslib._py_ElementType.ELEM_TET
ELEM_UNDEFINED = stepslib._py_ElementType.ELEM_UNDEFINED
UNKNOWN_TET = stepslib.UNKNOWN_TET
UNKNOWN_TRI = stepslib.UNKNOWN_TRI
UNKNOWN_VERT = stepslib.UNKNOWN_VERT
INDEX_NUM_BYTES = stepslib.INDEX_NUM_BYTES
INDEX_DTYPE = numpy.uint32 if INDEX_NUM_BYTES == 4 else numpy.uint64
###################################################################################################
[docs]@nutils.FreezeAfterInit
class Geometry(nutils.UsableObject, nutils.StepsWrapperObject):
"""Top-level geometry container
A number of compartment objects and patches objects may be grouped.
Should be used as a context manager for the declaration of compartments etc::
geom = Geometry()
with geom:
comp1 = Compartment.Create(vsys, volume)
... # Declare other objects in geom
After having declared children objects in the ``with mesh:`` block, they can be accessed
as attributes of ``geom`` with their name (see :py:class:`steps.API_2.utils.NamedObject` and
:py:func:`steps.API_2.utils.NamedObject.Create`)::
geom.comp1
"""
def __init__(self, *args, _createObj=True, **kwargs):
super().__init__(*args, **kwargs)
self.stepsGeom = stepslib._py_Geom() if _createObj else None
self._local = False
self._owned = None
self._callKwargs = {}
self._lstArgs = {}
def _SetUpMdlDeps(self, mdl):
"""Set up structures that depend on objects declared in the model."""
# Start with patches because they can update compartments
for pl in self._getChildrenOfType(Patch):
pl._SetUpMdlDeps(mdl)
for pl in self._getChildrenOfType(Membrane):
pl._SetUpMdlDeps(mdl)
for pl in self._getChildrenOfType(Compartment):
pl._SetUpMdlDeps(mdl)
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsGeom]
@classmethod
def _FromStepsObject(cls, obj):
"""Create the interface object from a STEPS object."""
geom = cls(_createObj=False)
geom.stepsGeom = obj
with geom:
for comp in obj.getAllComps():
Compartment._FromStepsObject(comp, geom)
for patch in obj.getAllPatches():
Patch._FromStepsObject(patch, geom)
return geom
class _PhysicalLocation(nutils.UsingObjects(Geometry), nutils.StepsWrapperObject):
"""Base class for physical locations
Compartment, Patch, Membrane, etc. derive from this class
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.systems = []
self.sysNames = []
(self._geom,) = self._getUsedObjects()
def addSystem(self, sys, _loc=None):
"""Add a surface or a volume system to the location
:param sys: The volume or surface system to be added, or its name.
:type sys: :py:class:`steps.API_2.model.VolumeSystem`, :py:class:`steps.API_2.model.SurfaceSystem`\
or `str`
:returns: None
"""
if isinstance(sys, nmodel.SpaceSystem):
self.sysNames.append((sys.name, _loc))
elif isinstance(sys, str):
self.sysNames.append((sys, _loc))
else:
raise TypeError(f'Expected a VolumeSystem or a SurfaceSystem, got a {type(sys)} instead.')
def _SetUpMdlDeps(self, mdl):
"""Set up the structures that will allow species or reaction access from locations."""
for sname, loc in self.sysNames:
s = getattr(mdl, sname)
self.systems.append((s, loc))
# Add the reactions and the corresponding reactants to children
for name, c in s.children.items():
if _PhysicalLocation._canBeChild(self, c):
if self._canBeChild(c):
self._addChildren(c)
for re in c._getAllElems(loc):
if re.name not in self.children:
self._addChildren(re)
s._addLocation(self)
def _canBeChild(self, c):
"""Return wether c can be a child of self."""
return isinstance(c, (
nmodel.Reaction, nmodel.Diffusion, nmodel.Current, nmodel.Complex, nmodel.VesicleBind,
nmodel.VesicleUnbind, nmodel.Endocytosis, nmodel.RaftGen
))
def _createStepsObj(self, geom):
pass
def _simPathAutoMetaData(self):
"""Return a dictionary with string keys and string or numbers values."""
return {'loc_type': self.__class__._locStr, 'loc_id': self.name}
@property
def _callKwargs(self):
return self._geom._callKwargs
@property
def _local(self):
return self._geom._local
@property
def _owned(self):
return self._geom._owned
[docs]@nutils.FreezeAfterInit
class Compartment(_PhysicalLocation, nutils.ParameterizedObject, nutils.Facade):
"""Base class for compartment objects
The same class is used to declare compartments in well-mixed models and in 3D tetrahedral
meshes.
Compartments in well-mixed geometry can be declared in the following ways::
with geom:
wmComp1 = Compartment.Create()
wmComp2 = Compartment.Create(vsys)
wmComp3 = Compartment.Create(vol=volume)
wmComp4 = Compartment.Create(vsys, volume)
:param vsys: Optional, the volume system associated with this compartment.
:type vsys: :py:class:`steps.API_2.model.VolumeSystem` or :py:class:`str`
:param vol: Optional, volume of the compartment
:type vol: float
Compartments in a tetrahedral mesh can be declared in the following ways::
with mesh:
tetComp1 = Compartment.Create(tetLst)
tetComp2 = Compartment.Create(tetLst, vsys)
:param tetLst: List of tetrahedrons associated with the compartment.
:type tetLst: :py:class:`TetList` or any argument that can be used to build one
:param vsys: Optional, the volume system associated with this compartment.
:type vsys: :py:class:`steps.API_2.model.VolumeSystem` or :py:class:`str`
The volume of a compartment in a tetrahedral mesh is the total volume of the encapsulated
tetrahedrons.
Compartments in a distributed tetrahedral mesh can be declared in the following ways::
with mesh:
distComp1 = Compartment.Create()
distComp2 = Compartment.Create(vsys, conductivity=g)
distComp3 = Compartment.Create(physicalTag=10)
distComp4 = Compartment.Create(tetLst, vsys)
:param vsys: Optional, the volume system associated with this compartment.
:type vsys: :py:class:`steps.API_2.model.VolumeSystem` or :py:class:`str`
:param conductivity: Optional, the volume conductivity of the compartment.
:type conductivity: float
:param physicalTag: Optional, the index of the physical tag used to build the compartment
(see :py:class:`DistMesh`).
:type physicalTag: int
:param tetLst: Optional, List of tetrahedrons associated with the compartment. If it is ommited
or an empty list is given, the compartment name will be used to find tetrahedrons tagged with it.
:type tetLst: :py:class:`TetList` or any argument that can be used to build one
Some of the methods documented below are only available for tetrahedral compartments, they are
grouped accordingly.
"""
_locStr = 'Comp'
def __new__(cls, *args, **kwargs):
# Compartment is a facade to its subclasses, the user always uses the Compartment class
# but the appropriate subclass is used for object instantiation.
if cls is not Compartment:
return super(Compartment, cls).__new__(cls)
geom = nutils.UsableObject._getUsedObjectOfClass(Geometry)
if geom is None:
raise Exception(f'Compartments need to be declared in a geometry.')
if geom.__class__ == Geometry:
return _WmCompartment.__new__(_WmCompartment, *args, **kwargs)
elif geom.__class__ == TetMesh:
return _TetCompartment.__new__(_TetCompartment, *args, **kwargs)
elif geom.__class__ == DistMesh:
return _DistCompartment.__new__(_DistCompartment, *args, **kwargs)
else:
raise NotImplementedError(f'Compartments cannot be created in {type(geom)}')
@classmethod
def _getDisplayName(cls):
return Compartment.__name__
@classmethod
def _FromStepsObject(cls, obj, geom):
if stepslib._STEPS_USE_DIST_MESH and isinstance(obj, stepslib._py_DistComp):
return _DistCompartment._FromStepsObject(obj, geom)
elif isinstance(obj, stepslib._py_TmComp):
return _TetCompartment._FromStepsObject(obj, geom)
elif isinstance(obj, stepslib._py_Comp):
return _WmCompartment._FromStepsObject(obj, geom)
else:
raise NotImplementedError(f'Cannot create API2 version of {type(obj)}')
def __init__(self, vsys, _createObj=True, **kwargs):
super().__init__(**kwargs)
self.stepsComp = self._createStepsObj(self._geom) if _createObj else None
if vsys is not None:
self.addSystem(vsys)
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsComp]
def _createStepsObj(self, geom):
raise NotImplementedError()
def _SetUpMdlDeps(self, mdl):
"""Set up the structures that will allow species or reaction access from locations."""
super()._SetUpMdlDeps(mdl)
for ves in mdl._getChildrenOfType(nmodel.Vesicle):
self._addChildren(ves)
def _canBeChild(self, c):
"""Return wether c can be a child of self."""
if not super()._canBeChild(c):
return False
if ((isinstance(c, nmodel.Reaction) and c._isSurfaceReac()) or
(isinstance(c, nmodel.Diffusion) and c._isSurfaceDiff()) or
isinstance(c, (nmodel.Endocytosis, nmodel.RaftGen))
):
return False
return True
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return Compartment._locStr
[docs] def addSystem(self, vsys, _loc=None):
"""Add a volume system to the compartment
:param vsys: The volume system to be added, or its name.
:type vsys: :py:class:`steps.API_2.model.VolumeSystem` or `str`
:returns: None
"""
super().addSystem(vsys, _loc)
if _loc is None:
if isinstance(vsys, str):
self.stepsComp.addVolsys(vsys)
elif vsys.__class__ is nmodel.VolumeSystem:
self.stepsComp.addVolsys(vsys.name)
@property
@nutils.ParameterizedObject.RegisterGetter(units=nutils.Units('m^3'))
@nutils.IgnoresGhostElements
def Vol(self):
"""Volume of the compartment
:type: Union[float, :py:class:`steps.API_2.utils.Parameter`], read-only for tetrahedral compartments
"""
return self.stepsComp.getVol(**self._callKwargs)
@nutils.FreezeAfterInit
class _WmCompartment(Compartment):
def __init__(self, vsys=None, vol=None, **kwargs):
super().__init__(vsys=vsys, **kwargs)
if vol is not None:
self.Vol = vol
def _createStepsObj(self, geom):
return stepslib._py_Comp(self.name, geom.stepsGeom)
@classmethod
def _FromStepsObject(cls, obj, geom):
"""Create the interface object from a STEPS object."""
comp = cls(_createObj=False, name=obj.getID())
comp.stepsComp = obj
# setup volume systems
for vsysname in comp.stepsComp.getVolsys():
comp.addSystem(vsysname)
return comp
@Compartment.Vol.setter
@nutils.ParameterizedObject.RegisterSetter(units=nutils.Units('m^3'))
def Vol(self, v):
self.stepsComp.setVol(v)
@nutils.FreezeAfterInit
class _BaseTetCompartment(Compartment):
_FACADE_TITLE_STR = 'Only available for compartments defined in tetrahedral meshes'
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
@property
def surface(self):
"""All triangles in the surface of the compartment
:type: :py:class:`TriList`, read-only
"""
return self.tets.surface
@property
def tets(self):
"""All tetrahedrons in the compartment
:type: :py:class:`TetList`, read-only
"""
raise NotImplementedError()
@property
@nutils.IgnoresGhostElements
def bbox(self):
"""The bounding box of the compartment
:type: :py:class:`steps.API_2.geom.BoundingBox`, read-only
"""
return BoundingBox(
Point(*self.stepsComp.getBoundMin(**self._callKwargs)),
Point(*self.stepsComp.getBoundMax(**self._callKwargs))
)
@nutils.FreezeAfterInit
class _TetCompartment(_BaseTetCompartment):
_FACADE_TITLE_STR = None
def __init__(self, tetLst, vsys=None, **kwargs):
self._tetLst = tetLst
super().__init__(vsys=vsys, **kwargs)
def _createStepsObj(self, geom):
self._tetLst = TetList._toRefList(self._tetLst, geom)
return stepslib._py_TmComp(self.name, geom.stepsMesh, self._tetLst.indices)
@classmethod
def _FromStepsObject(cls, obj, geom):
"""Create the interface object from a STEPS object."""
tetLst = TetList(obj.getAllTetIndices(), geom)
comp = cls(tetLst, _createObj=False, name=obj.getID())
comp.stepsComp = obj
# setup volume systems
for vsysname in comp.stepsComp.getVolsys():
comp.addSystem(vsysname)
return comp
@_BaseTetCompartment.tets.getter
def tets(self):
return self._tetLst
@nutils.FreezeAfterInit
class _DistCompartment(_BaseTetCompartment):
_FACADE_TITLE_STR = 'Only available for compartments defined in distributed tetrahedral meshes'
def __init__(self, *args, conductivity=0, physicalTag=None, **kwargs):
self._physicalTag = physicalTag
# Allow any ordering of vsys and tetLst
vsys, tetLst, args, kwargs = nutils.extractArgs(args, kwargs, [
('vsys', (nmodel.VolumeSystem, str), None),
('tetLst', (TetList, list, tuple), None),
])
self._tetLst = tetLst
super().__init__(vsys=vsys, **kwargs)
self.Conductivity = conductivity
def _createStepsObj(self, geom):
if self._tetLst is None:
comp = stepslib._py_DistComp(self.name, geom.stepsMesh, None, self._physicalTag)
else:
tetLst = TetList._toRefList(self._tetLst, geom)
if tetLst.isLocal() and tetLst._owned != False:
raise Exception(
f'The compartment is being created with a local tetrahedron list but without '
f'including non-owned elements. It should instead be created inside a with '
f'mesh.asLocal(owned=False) block.'
)
comp = stepslib._py_DistComp(
self.name, geom.stepsMesh, tetLst.indices, self._physicalTag, local=tetLst.isLocal()
)
del self._tetLst
return comp
@_BaseTetCompartment.tets.getter
def tets(self):
return TetList(
self.stepsComp.getAllTetIndices(**self._geom._callKwargs),
mesh=self._geom, _immutable=True, **self._geom._lstArgs
)
@_BaseTetCompartment.surface.getter
@nutils.IgnoresGhostElements
def surface(self):
return TriList(
self.stepsComp.getSurfTris(**self._geom._callKwargs),
mesh=self._geom, _immutable=True, **self._geom._lstArgs
)
@property
@nutils.ParameterizedObject.RegisterGetter(units=nutils.Units('S m^-1'))
def Conductivity(self):
"""Conductivity of the compartment
:type: Union[float, :py:class:`steps.API_2.utils.Parameter`]
"""
return self.stepsComp.getConductivity()
@Conductivity.setter
@nutils.ParameterizedObject.RegisterSetter(units=nutils.Units('S m^-1'))
def Conductivity(self, v):
self.stepsComp.setConductivity(v)
[docs]@nutils.FreezeAfterInit
class Patch(_PhysicalLocation, nutils.UsableObject, nutils.ParameterizedObject, nutils.Facade):
"""Base class for patch objects
A patch is a piece of 2D surface surrounding (part of) a 3D compartment, which may be
connected to another compartment. The same class is used to declare patches in well-mixed
models and in 3D tetrahedral meshes.
Patches in well-mixed geometry can be declared in the following ways::
with geom:
patch1 = Patch.Create(inner)
patch2 = Patch.Create(inner, outer)
patch3 = Patch.Create(inner, outer, ssys)
patch4 = Patch.Create(inner, None , ssys)
patch5 = Patch.Create(inner, None , None, area)
patch6 = Patch.Create(inner, outer, None, area)
patch7 = Patch.Create(inner, outer, ssys, area)
:param inner: Inner compartment, see below
:type inner: :py:class:`steps.API_2.geom.Compartment`
:param outer: Optional, Outer compartment, see below
:type outer: :py:class:`steps.API_2.geom.Compartment`
:param ssys: Optional, the surface system associated with this patch.
:type ssys: :py:class:`steps.API_2.model.SurfaceSystem` or :py:class:`str`
:param area: Optional, surface area of the patch
:type area: float
Patches in a tetrahedral mesh can be declared in the following ways::
with mesh:
patch1 = Patch.Create(triLst, inner)
patch2 = Patch.Create(triLst, inner, outer)
patch3 = Patch.Create(triLst, inner, None , ssys)
patch4 = Patch.Create(triLst, inner, outer, ssys)
:param triLst: List of triangles associated with the patch.
:type triLst: :py:class:`TriList` or any argument that can be used to build one
:param inner: Inner compartment, see below
:type inner: :py:class:`steps.API_2.geom.Compartment`
:param outer: Optional, Outer compartment, see below
:type outer: :py:class:`steps.API_2.geom.Compartment`
:param ssys: Optional, the surface system associated with this patch.
:type ssys: :py:class:`steps.API_2.model.SurfaceSystem` or :py:class:`str`
Patches in a distributed tetrahedral mesh can be declared in the following ways::
with mesh:
comp1 = Patch.Create(inner)
comp2 = Patch.Create(inner, ssys=ssys)
comp3 = Patch.Create(inner, outer, ssys)
comp4 = Patch.Create(inner, physicalTag=10)
comp5 = Patch.Create(triList, inner, ssys=ssys)
:param inner: Inner compartment, see below
:type inner: :py:class:`steps.API_2.geom.Compartment`
:param outer: Optional, Outer compartment, see below
:type outer: :py:class:`steps.API_2.geom.Compartment`
:param ssys: Optional, the surface system associated with this patch.
:type ssys: :py:class:`steps.API_2.model.SurfaceSystem` or :py:class:`str`
:param physicalTag: Optional, the index of the physical tag used to build the patch
(see :py:class:`DistMesh`).
:type physicalTag: int
:param triLst: Optional, list of triangles associated with the patch. If an empy list is given
the patch will be created from triangles that are tagged with the patch name.
:type triLst: :py:class:`TriList` or any argument that can be used to build one
**Relationship between Compartments and Patches**
It is necessary to explain the inner/outer relationship between compartments
and patches. When a patch object is created. it is necessary to arbitrarily label the
compartment(s) "inner" and "outer" (if a patch is connected to only one compartment
then the compartment must be labelled "inner" by convention). This is necessary
in order to fully describe the surface reaction rules. Accordingly, compartments
also store a list of connections, "inner" patches and "outer" patches. So if a
patch1 is created with comp1 as it's "inner" compartment, comp1 knows patch1 as
an "outer" patch. The labelling is purely defined when creating the Patch objects,
bearing in mind the stoichiometry defined in the surface reaction objects. This may
seem a little confusing at first, but will become clearer when experience is gained
with these objects.
For patches in tetrahedral meshes, the area is the total area of the encapsulated triangles.
Some of the methods documented below are only available for tetrahedral patches, they are
grouped accordingly.
"""
_locStr = 'Patch'
def __new__(cls, *args, **kwargs):
# Patch is a facade to its subclasses, the user always uses the Patch class
# but the appropriate subclass is used for object instantiation.
if cls is not Patch:
return super(Patch, cls).__new__(cls)
geom = nutils.UsableObject._getUsedObjectOfClass(Geometry)
if geom is None:
raise Exception(f'Patches need to be declared in a geometry.')
if geom.__class__ == Geometry:
return _WmPatch.__new__(_WmPatch, *args, **kwargs)
elif geom.__class__ == TetMesh:
return _TetPatch.__new__(_TetPatch, *args, **kwargs)
elif geom.__class__ == DistMesh:
return _DistPatch.__new__(_DistPatch, *args, **kwargs)
else:
raise NotImplementedError(f'Patches cannot be created in {type(geom)}')
@classmethod
def _getDisplayName(cls):
return Patch.__name__
@classmethod
def _FromStepsObject(cls, obj, geom):
if stepslib._STEPS_USE_DIST_MESH and isinstance(obj, stepslib._py_DistPatch):
return _DistPatch._FromStepsObject(obj, geom)
elif isinstance(obj, stepslib._py_TmPatch):
return _TetPatch._FromStepsObject(obj, geom)
elif isinstance(obj, stepslib._py_Patch):
return _WmPatch._FromStepsObject(obj, geom)
else:
raise NotImplementedError(f'Cannot create API2 version of {type(obj)}')
def __init__(self, inner, outer=None, ssys=None, _createObj=True, **kwargs):
super().__init__(**kwargs)
if not isinstance(inner, Compartment):
raise TypeError(f'Expected an inner Compartment, got {inner} instead.')
if outer is not None and not isinstance(outer, Compartment):
raise TypeError(f'Expected a Compartment, got {outer} instead.')
self._innerComp = inner
self._outerComp = outer
icomp = self._innerComp.stepsComp if isinstance(self._innerComp, Compartment) else None
ocomp = self._outerComp.stepsComp if isinstance(self._outerComp, Compartment) else None
self.stepsPatch = self._createStepsObj(icomp, ocomp, self._geom) if _createObj else None
if ssys is not None:
self.addSystem(ssys)
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsPatch]
[docs] def addSystem(self, ssys):
"""Add a surface system to the patch
:param ssys: The surface system to be added, or its name.
:type ssys: :py:class:`steps.API_2.model.SurfaceSystem` or `str`
:returns: None
"""
super().addSystem(ssys, nmodel.Location.SURF)
if isinstance(ssys, str):
self.stepsPatch.addSurfsys(ssys)
elif ssys.__class__ is nmodel.SurfaceSystem:
self.stepsPatch.addSurfsys(ssys.name)
def _createStepsObj(self, icomp, ocomp, geom):
raise NotImplementedError()
def _canBeChild(self, c):
"""Return wether c can be a child of self."""
if not super()._canBeChild(c):
return False
if ((isinstance(c, nmodel.Reaction) and not c._isSurfaceReac()) or
(isinstance(c, nmodel.Diffusion) and not c._isSurfaceDiff()) or
isinstance(c, (nmodel.VesicleBind, nmodel.VesicleUnbind))
):
return False
return True
def _SetUpMdlDeps(self, mdl):
"""Set up the structures that will allow species or reaction access from locations."""
super()._SetUpMdlDeps(mdl)
for raft in mdl._getChildrenOfType(nmodel.Raft):
self._addChildren(raft)
for endoZone in self._getChildrenOfType(EndocyticZone):
endoZone._SetUpMdlDeps(mdl)
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return Patch._locStr
@property
@nutils.ParameterizedObject.RegisterGetter(units=nutils.Units('m^2'))
@nutils.IgnoresGhostElements
def Area(self):
"""Surface area of the patch
:type: Union[float, :py:class:`steps.API_2.utils.Parameter`], read-only for tetrahedral patches
"""
return self.stepsPatch.getArea(**self._callKwargs)
@property
def innerComp(self):
"""Inner compartment
:type: :py:class:`Compartment`
"""
return self._innerComp
@property
def outerComp(self):
"""Outer compartment
:type: :py:class:`Compartment`
"""
return self._outerComp
@nutils.FreezeAfterInit
class _BasePatch(Patch):
@classmethod
def _FromStepsObject(cls, obj, geom):
"""Create the interface object from a STEPS object."""
args = cls._ArgsFromStepsObject(obj, geom)
patch = cls(*args, _createObj=False, name=obj.getID())
patch.stepsPatch = obj
# setup surface systems
for ssysname in patch.stepsPatch.getSurfsys():
patch.addSystem(ssysname)
return patch
@classmethod
def _ArgsFromStepsObject(cls, obj, geom):
icomp = getattr(geom, obj.icomp.getID())
ocomp = getattr(geom, obj.ocomp.getID()) if obj.ocomp is not None else None
return (icomp, ocomp)
@nutils.FreezeAfterInit
class _WmPatch(_BasePatch):
def __init__(self, inner, outer=None, ssys=None, area=None, **kwargs):
super().__init__(inner, outer, ssys=ssys, **kwargs)
if area is not None:
self.Area = area
def _createStepsObj(self, icomp, ocomp, geom):
return stepslib._py_Patch(self.name, geom.stepsGeom, icomp, ocomp)
@Patch.Area.setter
@nutils.ParameterizedObject.RegisterSetter(units=nutils.Units('m^2'))
def Area(self, v):
self.stepsPatch.setArea(v)
@nutils.FreezeAfterInit
class _BaseTetPatch(_BasePatch):
_FACADE_TITLE_STR = 'Only available for patches defined in tetrahedral meshes'
def __init__(self, triLst, *args, **kwargs):
self._triLst = triLst
super().__init__(*args, **kwargs)
@property
@nutils.IgnoresGhostElements
def bbox(self):
"""The bounding box of the patch
:type: :py:class:`steps.API_2.geom.BoundingBox`, read-only
"""
return BoundingBox(
Point(*self.stepsPatch.getBoundMin(**self._callKwargs)),
Point(*self.stepsPatch.getBoundMax(**self._callKwargs))
)
@property
def tris(self):
"""All triangles in the patch
:type: :py:class:`TriList`, read-only
"""
raise NotImplementedError()
@property
def edges(self):
"""Perimeter of the triangles in the patch
I.e. all bars that are not shared by two triangles in the patch.
:type: :py:class:`BarList`, read-only
"""
return self._triLst.edges
@nutils.FreezeAfterInit
class _TetPatch(_BaseTetPatch):
_FACADE_TITLE_STR = None
def __init__(self, triLst, inner, outer=None, ssys=None, **kwargs):
super().__init__(triLst, inner, outer, ssys=ssys, **kwargs)
def _createStepsObj(self, icomp, ocomp, geom):
if icomp is None and ocomp is None:
raise Exception(
f'Cannot declare a Patch in a tetrahedral mesh with no neighboring compartment.'
)
self._triLst = TriList._toRefList(self._triLst, geom)
return stepslib._py_TmPatch(self.name, geom.stepsMesh, self._triLst.indices, icomp, ocomp)
@classmethod
def _FromStepsObject(cls, obj, geom):
"""Create the interface object from a STEPS object."""
patch = super()._FromStepsObject(obj, geom)
with patch:
for subObj in obj.getAllEndocyticZones():
EndocyticZone._FromStepsObject(subObj, geom)
return patch
@classmethod
def _ArgsFromStepsObject(cls, obj, geom):
return (TriList(obj.getAllTriIndices(), geom), ) + _BasePatch._ArgsFromStepsObject(obj, geom)
@_BaseTetPatch.tris.getter
def tris(self):
return self._triLst
@nutils.FreezeAfterInit
class _DistPatch(_BaseTetPatch):
_FACADE_TITLE_STR = None
def __init__(self, *args, physicalTag=None, **kwargs):
self._physicalTag = physicalTag
# Allow any ordering of inner, outer, ssys and triLst
inner, outer, ssys, triLst, args, kwargs = nutils.extractArgs(args, kwargs, [
('inner', Compartment, None),
('outer', Compartment, None),
('ssys', (nmodel.SurfaceSystem, str), None),
('triLst', (TriList, list, tuple), None),
])
self._triLst = triLst
super().__init__(triLst, inner, outer, ssys=ssys, **kwargs)
def _createStepsObj(self, icomp, ocomp, geom):
if self._triLst is None:
patch = stepslib._py_DistPatch(self.name, geom.stepsMesh, None, icomp, ocomp, self._physicalTag)
else:
triLst = TriList._toRefList(self._triLst, geom)
if triLst.isLocal() and triLst._owned != False:
raise Exception(
f'The patch is being created with a local triangle list but without '
f'including non-owned elements. It should instead be created inside a with '
f'mesh.asLocal(owned=False) block.'
)
patch = stepslib._py_DistPatch(
self.name, geom.stepsMesh, triLst.indices, icomp, ocomp, self._physicalTag, local=triLst.isLocal()
)
del self._triLst
return patch
@_BaseTetPatch.tris.getter
def tris(self):
return TriList(
self.stepsPatch.getAllTriIndices(**self._geom._callKwargs),
mesh=self._geom, _immutable=True, **self._geom._lstArgs
)
@_BaseTetPatch.edges.getter
def edges(self):
raise NotImplementedError()
[docs]@nutils.FreezeAfterInit
class EndocyticZone(nutils.UsingObjects(Patch)):
"""A set of triangles that model a zone in which endocytosis reactions can happen
Endocytosis reactions declared in surface systems do not happen at any point of the corresponding
patches, they only happen in endocytic zones. By default all endocytosis events are active in this
zone, but they can be deactivated during simulation.
:param tris: The list of triangle
:type tris: :py:class:`TriList`
Endocytic zones are declared by using a patch as a context manager:
with patch1:
zone1 = EndocyticZone.Create(triLst)
Note that the triangles in `triLst` all need to be part of `patch1`.
"""
def __init__(self, tris, _createObj=True, *args, **kwargs):
super().__init__(*args, **kwargs)
(patch,) = self._getUsedObjects()
if not isinstance(patch, _TetPatch):
raise TypeError('Endocytic zones can only be declared in patches from tetrahedral meshes.')
if not isinstance(tris, TriList):
mesh = patch._getParentOfType(Geometry)
if not isinstance(mesh, TetMesh):
raise TypeError(f'{patch} was not declared in a tetrahedral mesh.')
tris = TriList(tris, mesh=mesh)
self._patch = patch
self._tris = tris
if _createObj:
self.stepsEndoZone = stepslib._py_EndocyticZone(
self.name, self._patch._getStepsObjects()[0], self._tris.indices
)
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsEndoZone]
@classmethod
def _FromStepsObject(cls, obj, geom):
"""Create the interface object from a STEPS object."""
endoZone = cls(obj.getAllTriIndices(), _createObj=False, name=obj.getID())
endoZone.stepsEndoZone = obj
return endoZone
def _SetUpMdlDeps(self, mdl):
"""Set up the structures that will allow species or reaction access from locations."""
for endo in self._patch._getChildrenOfType(nmodel.Endocytosis):
self._addChildren(endo)
@property
def tris(self):
"""All triangles in the endocytic zone
:type: :py:class:`TriList`, read-only
"""
return self._tris
[docs]@nutils.FreezeAfterInit
class Membrane(_PhysicalLocation, nutils.Facade):
"""A set of patches on which membrane potential will be computed
This class provides annotation for a group of triangles that comprise
a surface describing a membrane in a Tetmesh. This may be the same
description of one or several TmPatches in order for voltage-dependent
transitions, currents and so on to be inserted in the membrane.
A Membrane object must be available if membrane potential calculation is to be
performed.
By default STEPS only adds the inner compartments of the patches to the conduction volume,
if other compartments should also be added to the conduction volume, they should be supplied
through the supplementary_comps keyword parameter.
:param patches: List of patches whose triangles will be added to the membrane
(only one patch for distributed meshes)
:type patches: :py:class:`list`
The following parameters are only available for non-distributed meshes:
:param verify: Perform some checks on suitability of membrane (see below)
:type verify: :py:class:`bool`
:param opt_method: Optimization method (see below)
:type opt_method: int
:param search_percent: See below
:type search_percent: float
:param opt_file_name: Full path to a membrane optimization file
:type opt_file_name: str
:param supplementary_comps: Supplementary compartments to be considered part of the conduction volume
:type supplementary_comps: List[:py:class:`Compartment`]
Perform some checks on suitability of membrane if *verify* is True: these checks will print
warnings if membrane forms an open surface or if any triangle is found to have
more than 3 neighbours.
Specify optimization method with *opt_method*:
1. principal axis ordering (quick to set up but usually results in slower simulation than
method 2).
2. breadth first search (can be time-consuming to set up, but usually faster simulation.
If breadth first search is chosen then argument *search_percent* can specify the number of
starting points to search for the lowest bandwidth.
If a filename (with full path) is given in optional argument *opt_file_name* the membrane
optimization will be loaded from file, which was saved previously for this membrane with
:py:func:`steps.API_2.sim.Simulation.saveMembOpt`
"""
_locStr = 'Memb'
def __new__(cls, *args, **kwargs):
# Membrane is a facade to its subclasses, the user always uses the Membrane class
# but the appropriate subclass is used for object instantiation.
if cls is not Membrane:
return super(Membrane, cls).__new__(cls)
geom = nutils.UsableObject._getUsedObjectOfClass(Geometry)
if geom is None:
raise Exception(f'Membranes need to be declared in a geometry.')
if geom.__class__ == TetMesh:
return _TetMembrane.__new__(_TetMembrane, *args, **kwargs)
elif geom.__class__ == DistMesh:
return _DistMembrane.__new__(_DistMembrane, *args, **kwargs)
else:
raise NotImplementedError(f'Membranes cannot be created in {type(geom)}')
@classmethod
def _getDisplayName(cls):
return Membrane.__name__
@classmethod
def _FromStepsObject(cls, obj, geom):
raise NotImplementedError(f'Cannot create API2 version of {type(obj)}')
def __init__(
self, patches, *args, **kwargs
):
super().__init__(*args, **kwargs)
(geom,) = self._getUsedObjects()
self.mesh = geom
self.patches = patches
self.stepsMemb = self._createStepsObj(geom)
def addSystem(self, _):
"""
Even though a Membrane is a PhysicalLocation, surface systems should be added to the
patches, not through the membrane.
:meta private:
"""
raise NotImplementedError(f'Cannot add space systems to a Membrane.')
def _SetUpMdlDeps(self, mdl):
"""Set up the structures that will allow species or reaction access from locations."""
for p in self.patches:
self.sysNames += p.sysNames
self.sysNames = list(set(self.sysNames))
super()._SetUpMdlDeps(mdl)
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsMemb]
def _createStepsObj(self, geom):
raise NotImplementedError()
@property
def open(self):
"""Whether the membrane is open (contains holes) or not
:type: `bool`, read-only
"""
return self.stepsMemb.open()
@property
def tris(self):
"""List of triangles associated with the membrane
:type: :py:class:`steps.API_2.geom.TriList`, read-only
"""
res = TriList([], mesh=self.mesh)
for p in self.patches:
res += p.tris
return res
@property
def Area(self):
"""Surface area of the membrane
:type: float, read-only
"""
return sum(p.Area for p in self.patches)
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return Membrane._locStr
@nutils.FreezeAfterInit
class _TetMembrane(Membrane):
def __init__(
self, patches, verify=False, opt_method=1, search_percent=100.0, opt_file_name='',
supplementary_comps=[], *args, **kwargs
):
self._tmpParams = dict(
verify=verify,
opt_method=opt_method,
search_percent=search_percent,
opt_file_name=opt_file_name,
supplementary_comps=[comp._getStepsObjects()[0] for comp in supplementary_comps],
)
super().__init__(patches, *args, **kwargs)
def _createStepsObj(self, geom):
(geomObj,) = geom._getStepsObjects()
patches = [p._getStepsObjects()[0] for p in self.patches]
return stepslib._py_Memb(
self.name,
geomObj,
patches,
**self._tmpParams,
)
@nutils.FreezeAfterInit
class _DistMembrane(Membrane, nutils.ParameterizedObject):
_FACADE_TITLE_STR = 'Only available for membranes defined in distributed tetrahedral meshes'
def __init__(self, patches, capacitance=0, **kwargs):
super().__init__(patches, **kwargs)
if len(patches) != 1:
raise ValueError(f'Membranes in distributed meshes can only have a single patch.')
self.Capacitance = capacitance
def _createStepsObj(self, geom):
patches = [p._getStepsObjects()[0] for p in self.patches]
return stepslib._py_DistMemb(self.name, geom.stepsMesh, patches)
@property
@nutils.ParameterizedObject.RegisterGetter(units=nutils.Units('F m^-2'))
def Capacitance(self):
"""Capacitance of the compartment
:type: Union[float, :py:class:`steps.API_2.utils.Parameter`]
"""
return self.stepsMemb.getCapacitance()
@Capacitance.setter
@nutils.ParameterizedObject.RegisterSetter(units=nutils.Units('F m^-2'))
def Capacitance(self, v):
self.stepsMemb.setCapacitance(v)
[docs]@nutils.FreezeAfterInit
class Point(numpy.ndarray):
"""Convenience class for representing 3D points
This class inherits from :py:class:`numpy.ndarray` and can thus be used in the same way
as a numpy array. The only difference is the possibility to access invidual coordinates
through the *x*, *y*, and *z* properties.
:param x: x coordinate
:type x: float
:param y: y coordinate
:type y: float
:param z: z coordinate
:type z: float
"""
def __new__(cls, x, y, z):
arr = numpy.zeros(3)
arr[0:3] = x, y, z
return arr.view(cls)
def __init__(self, *args, **kwargs):
pass
@property
def x(self):
"""x coordinate
:type: float, read-only
"""
return self[0]
@property
def y(self):
"""y coordinate
:type: float, read-only
"""
return self[1]
@property
def z(self):
"""z coordinate
:type: float, read-only
"""
return self[2]
def __hash__(self):
return hash(tuple(self))
[docs]@nutils.FreezeAfterInit
class BoundingBox:
"""Convenience class for holding bounding box information"""
def __init__(self, minP, maxP):
self._min = minP
self._max = maxP
@property
def min(self):
"""Minimum point of the bounding box
:type: :py:class:`Point`, read-only
"""
return self._min
@property
def max(self):
"""Maximum point of the bounding box
:type: :py:class:`Point`, read-only
"""
return self._max
@property
def center(self):
"""Center point of the bounding box
:type: :py:class:`Point`, read-only
"""
return (self._max - self._min) / 2 + self._min
class _BaseTetMesh(Geometry):
"""Base class for tetrahedral meshes"""
def __init__(self, *args, _createObj=True, **kwargs):
super().__init__(*args, _createObj=False, **kwargs)
self.stepsMesh = self._createStepsObj() if _createObj else None
def _createStepsObj(self):
raise NotImplementedError()
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsMesh]
@property
def tets(self):
"""All tetrahedrons in the mesh
:type: :py:class:`TetList`, read-only
"""
return TetList(
range(self.stepsMesh.countTets(**self._callKwargs)), self, _immutable=True, **self._lstArgs
)
@property
def tris(self):
"""All triangles in the mesh
:type: :py:class:`TriList`, read-only
"""
return TriList(
range(self.stepsMesh.countTris(**self._callKwargs)), self, _immutable=True, **self._lstArgs
)
@property
def surface(self):
"""All triangles in the surface of the mesh
:type: :py:class:`TriList`, read-only
"""
return TriList(
self.stepsMesh.getSurfTris(**self._callKwargs), self, _immutable=True, **self._lstArgs
)
@property
def verts(self):
"""All vertices in the mesh
:type: :py:class:`VertList`, read-only
"""
return VertList(
range(self.stepsMesh.countVertices(**self._callKwargs)), self, _immutable=True, **self._lstArgs
)
@property
def bars(self):
"""All bars in the mesh
:type: :py:class:`BarList`, read-only
"""
return BarList(
range(self.stepsMesh.countBars(**self._callKwargs)), self, _immutable=True, **self._lstArgs
)
@property
@nutils.IgnoresGhostElements
def bbox(self):
"""The bounding box of the mesh
:type: :py:class:`steps.API_2.geom.BoundingBox`, read-only
"""
return BoundingBox(
Point(*self.stepsMesh.getBoundMin(**self._callKwargs)),
Point(*self.stepsMesh.getBoundMax(**self._callKwargs))
)
@property
@nutils.IgnoresGhostElements
def Vol(self):
"""Total volume of the mesh
:type: float, read-only
"""
return self.stepsMesh.getMeshVolume(**self._callKwargs)
def intersect(self, points, sampling=-1):
"""Computes the intersection of the current mesh and line segment(s), given their vertices
:param points: A 2-D NumPy array (/memview) of 3D points where each element contains the 3 point
coordinates
:type points: numpy.ndarray
:param sampling: if > 0, use montecarlo method with sampling points, otherwise use deterministic
method.
:type sampling: int
:returns: A list of lists of tuples representing the intersected tetrahedrons, one element per line
segment. Each tuple is made of 2 elements, a tetrahedron global identifier, and its respective
intersection ratio.
:rtype: List[List[Tuple[:py:class:`TetReference`, float]]]
"""
res = self.stepsMesh.intersect(points, sampling)
return [[(TetReference(idx, mesh=self, **self._lstArgs), rat) for idx, rat in seg] for seg in res]
def intersectIndependentSegments(self, points, sampling=-1):
"""Similar to the intersect method but here we deal with independent segments, i.e.
every two points we have a segment not related to previous or following ones.
E.g. seg0 = (points[0], points[1]), seg1 = (points[2], points[3]), etc.
:param points: A 2-D NumPy array (/memview) of 3D points where each element contains the 3 point
coordinates
:type points: numpy.ndarray
:param sampling: if > 0, use montecarlo method with sampling points, otherwise use deterministic
method.
:type sampling: int
:returns: A list of lists of tuples representing the intersected tetrahedrons, one element per line
segment. Each tuple is made of 2 elements, a tetrahedron global identifier, and its respective
intersection ratio.
:rtype: List[List[Tuple[:py:class:`TetReference`, float]]]
"""
res = self.stepsMesh.intersectIndependentSegments(points, sampling)
return [[(TetReference(idx, mesh=self, **self._lstArgs), rat) for idx, rat in seg] for seg in res]
[docs]@nutils.FreezeAfterInit
class TetMesh(_BaseTetMesh):
"""Container class for static tetrahedral meshes
This class stores the vertices points, 3D tetrahedral and 2D triangular elements that comprise
the mesh. The indices of the elements will be stored as unsigned integers (a positive integer
or zero) beginning at zero and incremented by 1. For example, if there are ntets number of
tetrahedrons in the mesh, the indices of the tetrahedrons will be [0,1,2,..., (ntets-1)].
A TetMesh object should not be created from scratch but should instead be created with one of
the dedicated class methods:
- :py:func:`TetMesh.Load`
- :py:func:`TetMesh.LoadAbaqus`
- :py:func:`TetMesh.LoadGmsh`
- :py:func:`TetMesh.LoadVTK`
- :py:func:`TetMesh.LoadTetGen`
TetMesh objects can be used as a context manager in the same way as
:py:class:`Geometry`::
mesh = TetMesh.Load('/path/to/file')
with mesh:
comp1 = Comp.Create(vsys)
... # Declare other objects in mesh
"""
def __init__(self, *args, _createObj=True, **kwargs):
self._vertGroups = {}
self._triGroups = {}
self._tetGroups = {}
self._vertProxy = None
self._triProxy = None
self._tetProxy = None
if _createObj:
raise Exception('Cannot create a bare TetMesh, use one of the class methods.')
super().__init__(*args, _createObj=False, **kwargs)
@classmethod
def _FromStepsObject(cls, obj, comps=None, patches=None, name=None):
"""Create the interface object from a STEPS object."""
mesh = cls(_createObj=False, name=name)
mesh.stepsMesh = obj
mesh.stepsGeom = obj
comps = obj.getAllComps() if comps is None else comps
patches = obj.getAllPatches() if patches is None else patches
with mesh:
for scomp in comps:
Compartment._FromStepsObject(scomp, mesh)
for spatch in patches:
Patch._FromStepsObject(spatch, mesh)
for ROIname in obj.getAllROINames():
ROI._FromStepsObject(obj.getROI(ROIname), mesh, ROIname)
return mesh
[docs] @classmethod
def FromData(cls, verts, tets, tris=[], name=None):
"""Create a mesh from geometrical data
Construct a Tetmesh container with the folowing method: Supply a list of all
vertices verts (by Cartesian coordinates), supply a list of all tetrahedrons
tets (by indices of the 4 vertices) and supply a full or partial list of
triangles tris (by indices of the 3 vertices). Indexing in STEPS begins at
zero, so the first 3 coordinates in verts will describe the zeroth vertex, t
he next 3 coordinates will describe the 1st vertex and so on. Labelling of
the vertices in tets and tris should follow this indexing. Lists must be
one-dimensional. Length of verts = nverts*3 where nverts is the total
number of vertices; length of tets = ntets*4 where ntets is the total
number of tetrahedrons; maximum length of tris ntris*3 where ntris is
the total number of triangles. For example, if we have just three tetrahedrons;
tet0=[0,1,2,3], tet1=[0,1,3,4] and tet2=[1,3,4,5] then the required
one-dimensional list tets=[0,1,2,3,0,1,3,4,1,3,4,5].
:param verts: Vertex data
:type verts: `List[int]`
:param tets: Tetrahedron data
:type tets: `List[int]`
:param tris: Optional triangle data
:type tris: `List[int]`
:param name: Optional name for the mesh
:type name: str
:returns: The constructed TetMesh
:rtype: :py:class:`TetMesh`
"""
mesh = cls(_createObj=False, name=name)
obj = stepslib._py_Tetmesh(verts, tets, tris)
mesh.stepsMesh = obj
mesh.stepsGeom = obj
return mesh
[docs] @classmethod
def Load(cls, path, scale=1, strict=False, name=None):
"""Load a mesh in STEPS
The mesh is loaded from an XML file that was previously generated by the
:py:func:`TetMesh.Save` method.
:param path: The root of the path where the file(s) are stored. e.g. with 'meshes/spine1'
this function will look for the file /meshes/spine1.xml
:type path: str
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param strict: Apply strict(-er) checking to the input XML
:type strict: bool
:param name: Optional name for the mesh
:type name: str
:returns: The loaded TetMesh
:rtype: :py:class:`TetMesh`
"""
smesh, scomps, spatches = smeshio.loadMesh(path, scale, strict)
return TetMesh._FromStepsObject(smesh, comps=scomps, patches=spatches, name=name)
def _loadElementProxys(self, ndprx, tetprx, triprx):
"""Load the blocks and groups from ElementProxy objects."""
self._vertProxy, self._triProxy, self._tetProxy = ndprx, triprx, tetprx
tmp = []
for prx, lstCls in zip([ndprx, triprx, tetprx], [VertList, TriList, TetList]):
# Merge blocks and groups
dct = {n: lstCls(range(b[0], b[1] + 1), mesh=self) for n, b in prx.getBlocks().items()}
for n, l in prx.getGroups().items():
if n in dct:
warnings.warn(
f'Key {n} is used by a {lstCls._refCls._locStr} block and a '
f'group, taking the group.'
)
dct[n] = lstCls(l, mesh=self)
tmp.append(dct)
self._vertGroups, self._triGroups, self._tetGroups = tmp
# The shadow_mesh keyword parameter is left here for compatibility but is not documented since the
# STEPS-CUBIT toolkit is deprecated.
[docs] @classmethod
def LoadAbaqus(cls, filename, scale=1, ebs=None, shadow_mesh=None, name=None):
"""Load a mesh from an ABAQUS-formated mesh file
If blocks or groups of elements are present in the file, they will also be loaded.
:param filename: The Abaqus filename (or path) including any suffix. Can also take a
2-tuple containing separate paths for tetrahedron and triangle data (in this order).
:type filename: str or `Tuple[str, str]`
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param ebs: Names of selected element blocks which are included in the mesh (does not apply
if a 2-tuple is given as filename)
:type ebs: `List[str]`
:param name: Optional name for the mesh
:type name: str
:returns: The loaded TetMesh
:rtype: :py:class:`TetMesh`
"""
if isinstance(filename, tuple) and len(filename) == 2:
tetfile, trifile = filename
stepsMesh, ndprx, tetprx, triprx = smeshio.importAbaqus2(
tetfile, trifile, scale, shadow_mesh=shadow_mesh
)
elif isinstance(filename, str):
stepsMesh, ndprx, tetprx, triprx = smeshio.importAbaqus(
filename, scale, ebs=ebs, shadow_mesh=shadow_mesh
)
else:
raise TypeError(f'Expected a string or a 2-tuple of strings, got {filename} instead.')
mesh = TetMesh._FromStepsObject(stepsMesh, name=name)
mesh._loadElementProxys(ndprx, tetprx, triprx)
return mesh
[docs] @classmethod
def LoadGmsh(cls, filename, scale=1, name=None):
"""Load a mesh from a Gmsh (2.2 ASCII)-formated mesh file
If blocks or groups of elements are present in the file, they will also be loaded.
:param filename: The Gmsh filename (or path) including any suffix.
:type filename: str
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param name: Optional name for the mesh
:type name: str
:returns: The loaded TetMesh
:rtype: :py:class:`TetMesh`
"""
stepsMesh, ndprx, tetprx, triprx = smeshio.importGmsh(filename, scale)
mesh = TetMesh._FromStepsObject(stepsMesh, name=name)
mesh._loadElementProxys(ndprx, tetprx, triprx)
return mesh
[docs] @classmethod
def LoadVTK(cls, filename, scale=1, name=None):
"""Load a mesh from a VTK (legacy ASCII)-formated mesh file
If blocks or groups of elements are present in the file, they will also be loaded.
:param filename: The VTK filename (or path) including any suffix.
:type filename: str
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param name: Optional name for the mesh
:type name: str
:returns: The loaded TetMesh
:rtype: :py:class:`TetMesh`
"""
stepsMesh, ndprx, tetprx, triprx = smeshio.importVTK(filename, scale)
mesh = TetMesh._FromStepsObject(stepsMesh, name=name)
mesh._loadElementProxys(ndprx, tetprx, triprx)
return mesh
[docs] @classmethod
def LoadTetGen(cls, pathroot, scale, name=None):
"""Load a mesh from a TetGen-formated set of files
If blocks or groups of elements are present, they will also be loaded.
:param pathroot: The root of the path to the mesh files. E.g. mesh/torus should be given
to read files mesh/torus.node, mesh/torus.ele, and optionally mesh/torus.face
:type pathroot: str
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param name: Optional name for the mesh
:type name: str
:returns: The loaded TetMesh
:rtype: :py:class:`TetMesh`
The TetGen file format for describing a mesh actually consists of
multiple files with varying suffixes. Currently, this class method only
reads meshes consisting of three files:
* <input>.node: describing the tetrahedral mesh node points.
* <input>.ele: describing tetrahedral elements, each of which
consists of 4 pointers into the <input>.node file. (TetGen
also supports 10-node elements; these 6 extra nodes are obviously
not read by STEPS.)
* <input>.face: describing triangular faces, each of which
consists of 3 pointers into the <input>.node file. This file is optional.
Other files are .vol (list of maximum volumes), .var (variant constraints)
.neigh (list of neighbours), .smesh (simple PLC descriptions) and .edge
(list of boundary edges) files. None of these seem relevant for our
use cases, so we don't load them even when they are there. In particular,
the .neigh file is computed by STEPS itself.
Please refer to the TetGen manual (pages 31-40 in the last edition)
for more information on these file formats
tetgen.berlios.de/files/tetgen-manual.pdf
"""
stepsMesh, ndprx, tetprx, triprx = smeshio.importTetGen(pathroot, scale)
mesh = TetMesh._FromStepsObject(stepsMesh, name=name)
mesh._loadElementProxys(ndprx, tetprx, triprx)
return mesh
[docs] def Save(self, path):
"""Save the mesh to an XML file
:param path: The root of the path to store the files. e.g. 'meshes/spine1' will save
data in /meshes/spine1.xml
:type path: str
This file stores the basic information about the mesh which tends to be
common information for any software that supports tetrahedral meshes.
* NODES are stored by cartesian coordinates.
* TRIANGLES are stored by the indices of their 3 nodes.
* TETRAHEDRONS are stored by the indices of their 4 nodes.
The XML file also stores information about any compartments or
patches created in STEPS (class :py:class:`steps.API_2.geom.Comp`
:py:class:`steps.API_2.geom.Patch` respectively).
* COMPARTMENT(S) are stored by:
* their string identification.
* a list of any volume systems added to the compartment at time of saving.
* a list of tetrahedrons belonging to the compartment
* PATCH(ES) are stored by:
* their string identification.
* a list of any surface systems added to the patch at time of saving.
* the inner compartment id.
* the outer compartment id (if it exists).
* a list of triangles belonging to this patch.
"""
smeshio.saveMesh(path, self.stepsMesh)
[docs] def ConvertToMetis(self, path):
"""Convert the mesh data to metis connectivity data
:param path: The path to the file in which to save the Metis data.
:type path: str
"""
smetis.tetmesh2metis(self.stepsMesh, path)
@property
def vertGroups(self):
"""Groups of vertices defined in the loaded mesh
:type: Dict[str, :py:class:`VertList`], read-only
.. note::
Both blocks and groups are accessed through the same property.
"""
return self._vertGroups
@property
def triGroups(self):
"""Groups of triangles defined in the loaded mesh
:type: Dict[str, :py:class:`TriList`], read-only
.. note::
Both blocks and groups are accessed through the same property.
"""
return self._triGroups
@property
def tetGroups(self):
"""Groups of tetrahedrons defined in the loaded mesh
:type: Dict[str, :py:class:`TetList`], read-only
.. note::
Both blocks and groups are accessed through the same property.
"""
return self._tetGroups
@property
def vertProxy(self):
"""Element proxy object for vertices (see :py:class:`steps.API_1.utilities.meshio.ElementProxy`)
:type: :py:class:`steps.API_1.utilities.meshio.ElementProxy`
.. note::
Only available with some import methods.
"""
return self._vertProxy
@property
def triProxy(self):
"""Element proxy object for triangles (see :py:class:`steps.API_1.utilities.meshio.ElementProxy`)
:type: :py:class:`steps.API_1.utilities.meshio.ElementProxy`
.. note::
Only available with some import methods.
"""
return self._triProxy
@property
def tetProxy(self):
"""Element proxy object for tetrahedrons (see :py:class:`steps.API_1.utilities.meshio.ElementProxy`)
:type: :py:class:`steps.API_1.utilities.meshio.ElementProxy`
.. note::
Only available with some import methods.
"""
return self._tetProxy
class _DistElemGroupsProxy:
"""Utility class for calling methods like DistMesh::getTaggedTetrahedrons"""
def __init__(self, mesh, lstCls, method):
self._mesh = mesh
self._lstCls = lstCls
self._method = method
def keys(self):
return self._mesh.stepsMesh.getTags(self._lstCls._dim)
def __getitem__(self, tag):
return self._lstCls(
self._method(tag, **self._mesh._callKwargs),
mesh=self._mesh, _immutable=True, **self._mesh._lstArgs
)
def __setitem__(self, n, v):
raise NotImplementedError(f'Cannot add group of elements to distributed meshes.')
def __delitem__(self, n):
raise NotImplementedError(f'Cannot remove group of elements from distributed meshes.')
[docs]@nutils.FreezeAfterInit
class DistMesh(_BaseTetMesh):
"""Container class for distributed tetrahedral meshes
In contrast to :py:class:`TetMesh`, only gmsh mesh files are accepted.
:param filename: The Gmsh filename (or path) including any suffix.
:type filename: str
:param scale: Optionally rescale the mesh on loading by given factor
:type scale: float
:param mpiComm: mpi4py MPI communicator, defaults to COMM_WORLD
:type mpiComm: :py:class:`mpi4py.MPI.Intracomm`
``filename`` should either be a full path to a single gmsh mesh or a path prefix for
pre-partitioned meshes. If a single gmsh file is provided, it will be partitioned automatically.
Pre-partitioned meshes can be loaded by providing a path prefix e.g. ``'path/to/meshName'`` and
individual partition files should be named ``'path/to/meshName_1.msh'``, ``'path/to/meshName_2.msh'``,
etc. The numbering of files starts at 1 and goes up to :py:attr:`steps.API_2.sim.MPI.nhosts`.
Note that, for pre-partitioned meshes, elements are not necessarily numbered between ``0`` and ``n``
(with ``n`` the number of elements).
"""
def __init__(self, filename, scale=1, mpiComm=None, _createObj=True, **kwargs):
self._filename = filename
if _createObj:
if mpiComm is None:
# Only import mpi4py when necessary
import mpi4py.MPI
mpiComm = mpi4py.MPI.COMM_WORLD
self._comm = mpiComm
# In dist steps, scale == 0 means no scaling
self._scale = 0 if scale == 1 else scale
super().__init__(_createObj=_createObj, **kwargs)
self._callKwargs['local'] = self._local
self._lstArgs['local'] = self._local
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsMesh]
def _createStepsObj(self):
library = stepslib._py_Library(self._comm)
return stepslib._py_DistMesh(library, self._filename, self._scale)
[docs] @contextlib.contextmanager
def asGlobal(self):
"""Context manager method to use the mesh globally
This should be used with the ``with`` statement in the following way::
with mesh.asGlobal():
lst1 = mesh.tets # Returns a global list of tetrahedron indices
lst2 = mesh.surface # Returns a global list of boundary triangles
It is equivalent to using only ``with mesh:`` since this defaults to using global indices.
``with mesh.asGlobal():`` is however preferable since it is more explicit.
:rtype: :py:class:`DistMesh`
"""
try:
self.__enter__()
yield self
finally:
self.__exit__(None, None, None)
[docs] @contextlib.contextmanager
def asLocal(self, owned=True):
"""Context manager method to use the mesh locally
This should be used with the ``with`` statement in the following way::
with mesh.asLocal():
lst1 = mesh.tets # Returns a local list of tetrahedron indices
lst2 = mesh.surface # Returns a local list of boundary triangles
All properties and methods that would usually return lists of global indices will return local
lists instead when wrapped with this context manager.
In the above example, ``lst1`` would end up having the same value as ``mesh.tets.toLocal()``
but the latter is less efficient since it first queries the full list of all tetrahedrons
and then restricts it to only the local ones.
.. note::
This does not affect existing lists or newly created lists. For example, a tetrahedron list
created in a ``with mesh.asLocal()`` block will still default to containing global elements
(unless the ``local=True`` keyword argument is given at creation).
This context manager only affects the behavior of the mesh object, and its subparts like
compartments and patches. For example ``comp.tets`` will return a list of local indices
in this context manager.
:param owned: Whether only owned elements should be considered (defaults to True)
:type owned: bool
:rtype: :py:class:`DistMesh`
"""
self._local = True
self._owned = owned
self._callKwargs['local'] = self._local
self._lstArgs['local'] = self._local
if not owned:
self._callKwargs['owned'] = self._owned
self._lstArgs['_owned'] = self._owned
try:
self.__enter__()
yield self
finally:
self.__exit__(None, None, None)
self._local = False
self._owned = None
self._callKwargs['local'] = self._local
self._lstArgs['local'] = self._local
if not owned:
del self._callKwargs['owned']
del self._lstArgs['_owned']
[docs] def intersect(self, points, sampling=-1, raw=False, local=True):
"""Computes the intersection of the current mesh and line segment(s), given their vertices
:param points: A 2-D NumPy array (/memview) of 3D points where each element contains the 3 point
coordinates
:type points: numpy.ndarray
:param sampling: if > 0, use montecarlo method with sampling points, otherwise use deterministic
method.
:type sampling: int
:param raw: If True, return raw integer tetrahedron indices.
:type raw: bool
:param local: if False, return global tetrahedron identifiers instead of local ones.
:type local: bool
This method can only be called when `mesh.asLocal()` is used, and it will return local tetrahedron
references by default.
:returns: A list of lists of tuples representing the intersected tetrahedrons, one element per line
segment. Each tuple is made of 2 elements, a tetrahedron identifier (local or global depending
on `local`), and its respective intersection ratio.
:rtype: List[List[Tuple[Union[:py:class:`TetReference`, int], float]]]
"""
if not self._local:
raise Exception('Cannot use intersect method without using "with mesh.asLocal():".')
res = self.stepsMesh.intersect(points, sampling)
if not local:
res = [[(self.stepsMesh.getTetGlobalIndex(idx), rat) for idx, rat in seg] for seg in res]
if raw:
return res
else:
return [[(TetReference(idx, mesh=self, local=local), rat) for idx, rat in seg] for seg in res]
[docs] def intersectIndependentSegments(self, points, sampling=-1, raw=False, local=True):
"""Similar to the intersect method but here we deal with independent segments, i.e.
every two points we have a segment not related to previous or following ones.
E.g. seg0 = (points[0], points[1]), seg1 = (points[2], points[3]), etc.
:param points: A 2-D NumPy array (/memview) of 3D points where each element contains the 3 point
coordinates
:type points: numpy.ndarray
:param sampling: if > 0, use montecarlo method with sampling points, otherwise use deterministic
method.
:type sampling: int
:param raw: If True, return raw integer tetrahedron indices.
:type raw: bool
:param local: if False, return global tetrahedron identifiers instead of local ones.
:type local: bool
This method can only be called when `mesh.asLocal()` is used, and it will return local tetrahedron
references by default.
:returns: A list of lists of tuples representing the intersected tetrahedrons, one element per line
segment. Each tuple is made of 2 elements, a tetrahedron identifier (local or global depending
on `local`), and its respective intersection ratio.
:rtype: List[List[Tuple[Union[:py:class:`TetReference`, int], float]]]
"""
if not self._local:
raise Exception('Cannot use intersectIndependentSegments method without using "with mesh.asLocal():".')
res = self.stepsMesh.intersectIndependentSegments(points, sampling)
if not local:
res = [[(self.stepsMesh.getTetGlobalIndex(idx), rat) for idx, rat in seg] for seg in res]
if raw:
return res
else:
return [[(TetReference(idx, mesh=self, local=local), rat) for idx, rat in seg] for seg in res]
@_BaseTetMesh.tets.getter
def tets(self):
"""All tetrahedrons in the mesh
:type: :py:class:`TetList`, read-only
"""
if self._local:
return TetList(
self.stepsMesh.getAllTetIndices(**self._callKwargs), self, _immutable=True, **self._lstArgs
)
else:
return _BaseTetMesh.tets.fget(self)
@_BaseTetMesh.tris.getter
def tris(self):
"""All triangles in the mesh
:type: :py:class:`TriList`, read-only
"""
if self._local:
return TriList(
self.stepsMesh.getAllTriIndices(**self._callKwargs), self, _immutable=True, **self._lstArgs
)
else:
return _BaseTetMesh.tris.fget(self)
@_BaseTetMesh.verts.getter
def verts(self):
"""All vertices in the mesh
:type: :py:class:`VertList`, read-only
"""
if self._local:
return VertList(
self.stepsMesh.getAllVertIndices(**self._callKwargs), self, _immutable=True, **self._lstArgs
)
else:
return _BaseTetMesh.verts.fget(self)
@property
def tetGroups(self):
return _DistElemGroupsProxy(self, TetList, self.stepsMesh.getTaggedTetrahedrons)
@property
def triGroups(self):
return _DistElemGroupsProxy(self, TriList, self.stepsMesh.getTaggedTriangles)
@property
def vertGroups(self):
return _DistElemGroupsProxy(self, VertList, self.stepsMesh.getTaggedVertices)
@staticmethod
def _use_gmsh():
return stepslib._py_DistMesh._use_gmsh()
DistMesh.intersect.__doc__ = _BaseTetMesh.intersect.__doc__
[docs]class Reference(nutils.UsingObjects(nutils.Optional(_BaseTetMesh)), nutils.SolverPathObject, nutils.Facade):
"""Base class for all element references.
:param idx: Index of the element
:type idx: int
:param mesh: Mesh object that contains the element
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
The actual index of the element can be accessed with the idx attribute.
"""
def __new__(cls, idx, mesh=None, _newCalled=False, **kwargs):
# TetReference, TriReference, etc. are facades to their subclasses, the user always uses
# the main class but the appropriate subclass is used for object instantiation.
if cls not in Reference.__subclasses__() + _DistReference.__subclasses__() or _newCalled:
return super(Reference, cls).__new__(cls)
idx, mesh = cls._getIdxAndMeshFromParams(idx, mesh)
if mesh.__class__ == DistMesh:
return cls._distCls.__new__(cls._distCls, idx, mesh=mesh, _newCalled=True, **kwargs)
else:
return cls.__new__(cls, idx, mesh=mesh, _newCalled=True, **kwargs)
def __init__(self, idx, mesh=None, anonymous=False, *args, **kwargs):
if not anonymous:
# Only call parent __init__ if we need the reference to be named
super().__init__(*args, addAsElement=False, **kwargs)
idx, mesh = self._getIdxAndMeshFromParams(idx, mesh)
self.mesh = mesh
self._idx = idx
self._callKwargs = {}
self._cloneArgs = dict(mesh=self.mesh)
[docs] def toList(self):
"""Get a list holding only this element
:returns: A list with only this element
:rtype: :py:class:`RefList`
"""
return self.__class__._lstCls([self], **self._cloneArgs)
@property
def idx(self):
"""Index of the element
:Type: int, read-only
"""
return self._idx
@classmethod
def _getIdxAndMeshFromParams(cls, idx, mesh):
if mesh is None:
# Try to infer the mesh from context-managers
mesh = nutils.UsableObject._getUsedObjectOfClass(_BaseTetMesh)
if mesh is None:
# Try to infer the mesh from the idx argument
if isinstance(idx, Reference):
mesh = idx.mesh
else:
raise Exception('Cannot infer which mesh is associated with the reference.')
if not isinstance(mesh, _BaseTetMesh):
raise TypeError(f'Expected a TetMesh or DistMesh object, got {mesh} instead.')
idx = cls._getActualIdx(idx, mesh)
return idx, mesh
@classmethod
def _getActualIdx(cls, idx, mesh):
"""Return the integer idx."""
if isinstance(idx, cls):
if mesh != idx.mesh:
raise Exception(
f'Cannot create a reference from a reference related to a different mesh.'
)
idx = idx._idx
if not isinstance(idx, numbers.Integral):
raise TypeError(f'Expected an integer index, got {idx} instead.')
return idx
def _getPhysicalLocation(self):
"""Return the location associated with the reference (comp for tet, patch for tri)."""
pass
def __hash__(self):
return hash((id(self.mesh), self._idx))
def __eq__(self, other):
return isinstance(other, self.__class__) and self.mesh == other.mesh and self._idx == other._idx
def __repr__(self):
return f'{self.__class__._locStr}({self._idx})'
def _solverId(self):
"""Return the id that is used to identify an object in the solver."""
return (self._idx,)
def _simPathAutoMetaData(self):
"""Return a dictionary with string keys and string or numbers values."""
mtdt = {'loc_type': self.__class__._locStr, 'loc_id': self._idx}
# Add information about parent comp or patch
parent = self._getPhysicalLocation()
if parent is not None:
for key, val in parent._simPathAutoMetaData().items():
mtdt['parent_' + key] = val
return mtdt
class _DistReference(Reference):
_FACADE_TITLE_STR = 'Only available for references defined in distributed tetrahedral meshes'
def __init__(self, idx, mesh=None, anonymous=False, *args, local=False, _owned=None, **kwargs):
super().__init__(idx, mesh, anonymous, *args, **kwargs)
if isinstance(idx, _DistReference):
local = idx.isLocal()
self._local = local
self._callKwargs['local'] = self._local
self._cloneArgs['local'] = self._local
# The owned flag is just indicative of how the reference was created
self._owned = _owned
@classmethod
def _getToLocalFunc(cls, mesh):
"""This is mainly used to avoid having to instantiate a VertReference when manipulating
vert indices.
"""
return getattr(mesh.stepsMesh, f'get{cls._locStr}LocalIndex')
@classmethod
def _getToGlobalFunc(cls, mesh):
return getattr(mesh.stepsMesh, f'get{cls._locStr}GlobalIndex')
def isLocal(self):
"""Return whether the reference uses a local index
:rtype: bool
"""
return self._local
def toLocal(self, owned=True):
"""Get the local version of the reference
If the reference is local, it returns itself. If the reference is global, and if the element
exists in the current MPI process, a local reference to it is returned. Otherwise,
`None` is returned. The returned reference is thus different for each MPI process.
:param owned: Whether only owned elements should be considered (defaults to True)
:type owned: bool
:rtype: :py:class:`Reference`
"""
if self._local:
if self._owned == owned:
return self
else:
return self.toGlobal().toLocal(owned)
method = self._getToLocalFunc(self.mesh)
ind = method(self.idx, owned=owned)
if ind is None:
return None
else:
return self.__class__(ind, mesh=self.mesh, local=True, _owned=owned)
def toGlobal(self, **_kwargs):
"""Get the global version of the reference
If the reference is not local, it returns itself. Otherwise, the local index of the original
reference is converted to global mesh index and the corresponding global reference is returned.
:rtype: :py:class:`Reference`
"""
if not self._local:
return self
method = self._getToGlobalFunc(self.mesh)
return self.__class__(method(self.idx), mesh=self.mesh, local=False, **_kwargs)
def _solverKeywordParams(self):
"""Return the additional keyword parameters that should be passed to the solver"""
return dict(local=self._local)
def _simPathAutoMetaData(self):
"""Return a dictionary with string keys and string or numbers values."""
mtdt = super()._simPathAutoMetaData()
if self._local:
# The location id in the metadata should always be global
# Get global id without creating a new reference because it causes an issue
# for VertReferences that querry their position upon creation
method = getattr(self.mesh.stepsMesh, f'get{self._locStr}GlobalIndex')
mtdt['loc_id'] = method(self.idx)
return mtdt
def __hash__(self):
return hash((id(self.mesh), self._idx, self._local))
def __eq__(self, other):
return isinstance(other, self.__class__) and\
(self.mesh, self._idx, self._local) == (other.mesh, other._idx, other._local)
def __repr__(self):
if self._local:
return f'Local{self.__class__._locStr}({self._idx})'
else:
return super().__repr__()
[docs]@nutils.FreezeAfterInit
class TetReference(Reference):
"""Convenience class for accessing properties of a tetrahedron
:param idx: Index of the tetrahedron
:type idx: int
:param mesh: Mesh object that contains the tetrahedron
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
"""
_locStr = 'Tet'
def __init__(self, idx, mesh=None, *args, **kwargs):
super().__init__(idx, mesh=mesh, *args, **kwargs)
@classmethod
def _getActualIdx(cls, idx, mesh):
"""Return the integer idx."""
try:
idx = super()._getActualIdx(idx, mesh)
except TypeError as ex:
if hasattr(idx, '__iter__') and len(idx) == 3:
idx = mesh.stepsMesh.findTetByPoint(list(idx))
if idx == UNKNOWN_TET:
raise Exception(f'No tetrahedron at position {idx}.')
else:
raise ex
return idx
def __getattr__(self, name):
"""Handle the references to species, reactions, or diffusion rules."""
if self.comp is not None:
return getattr(self.comp, name)
else:
raise AttributeError
@property
def children(self):
if self.comp is not None:
return self.comp.children
else:
return {}
[docs] def containsPoint(self, point):
"""Check whether a given 3D point is inside the tetrahedron
:param point: The 3D point
:type point: List[float]
:rtype: bool
"""
return self.mesh.stepsMesh.isPointInTet(point, self._idx, **self._callKwargs)
@property
def neighbs(self):
"""Neighboring tetrahedrons (between 0 and 4)
:type: :py:class:`TetList`, read-only
"""
inds = self.mesh.stepsMesh.getTetTetNeighb(self._idx, **self._callKwargs)
return TetList([ind for ind in inds if ind != UNKNOWN_TET], **self._cloneArgs)
@property
def faces(self):
"""Faces of the tetrahedron
:type: :py:class:`TriList`, read-only
"""
return TriList(self.mesh.stepsMesh.getTetTriNeighb(self._idx, **self._callKwargs), **self._cloneArgs)
@property
def center(self):
"""Barycenter of the tetrahedron
:type: :py:class:`Point`, read-only
"""
return Point(*self.mesh.stepsMesh.getTetBarycenter(self._idx, **self._callKwargs))
@property
def verts(self):
"""Vertices of the tetrahedron
:type: :py:class:`VertList`, read-only
"""
return VertList(self.mesh.stepsMesh.getTet(self._idx, **self._callKwargs), **self._cloneArgs)
@property
def comp(self):
"""Compartment associated to the tetrahedron (if any)
:type: :py:class:`Comp` or None, read-only
"""
try:
c = self.mesh.stepsMesh.getTetComp(self._idx, **self._callKwargs)
if c is None:
return None
else:
return getattr(self.mesh, c.getID())
except AttributeError as ex:
# Since comp is called in __getattr__, raising an AttributeError here can lead to
# infinite recursion.
raise Exception(ex)
@property
def Vol(self):
"""Volume of the tetrahedron
:type: float, read-only
"""
return self.mesh.stepsMesh.getTetVol(self._idx, **self._callKwargs)
@property
def qualityRER(self):
"""Radius-edge-ratio (a quality measurement) of the tetrahedron
:type: float, read-only
"""
return self.mesh.stepsMesh.getTetQualityRER(self._idx, **self._callKwargs)
def _getPhysicalLocation(self):
"""Return the location associated with the reference (comp for tet, patch for tri)."""
return self.comp
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return TetReference._locStr
class _DistTetReference(TetReference, _DistReference):
@TetReference.neighbs.getter
def neighbs(self):
if self.mesh._local:
inds = self.mesh.stepsMesh.getTetTetNeighb(self._idx, **self._callKwargs, owned=self.mesh._owned)
return TetList([ind for ind in inds if ind != UNKNOWN_TET], **self._cloneArgs)
else:
return TetReference.neighbs.fget(self)
TetReference._distCls = _DistTetReference
[docs]@nutils.FreezeAfterInit
class TriReference(Reference):
"""Convenience class for accessing properties of a triangle
:param idx: Index of the triangle
:type idx: int
:param mesh: Mesh object that contains the triangle
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
"""
_locStr = 'Tri'
def __init__(self, idx, mesh=None, *args, **kwargs):
super().__init__(idx, mesh=mesh, *args, **kwargs)
def __getattr__(self, name):
"""Handle the references to species or reactions."""
if self.patch is not None:
return getattr(self.patch, name)
else:
raise AttributeError
@property
def children(self):
if self.patch is not None:
return self.patch.children
else:
return {}
@property
def verts(self):
"""Vertices of the triangle
:type: :py:class:`VertList`, read-only
"""
return VertList(self.mesh.stepsMesh.getTri(self._idx, **self._callKwargs), **self._cloneArgs)
@property
def center(self):
"""Barycenter of the triangle
:type: :py:class:`Point`, read-only
"""
return Point(*self.mesh.stepsMesh.getTriBarycenter(self._idx, **self._callKwargs))
@property
def Area(self):
"""Area of the triangle
:type: float, read-only
"""
return self.mesh.stepsMesh.getTriArea(self._idx, **self._callKwargs)
@property
def tetNeighbs(self):
"""Neighboring tetrahedrons
:type: :py:class:`TetList`, read-only
"""
inds = self.mesh.stepsMesh.getTriTetNeighb(self._idx, **self._callKwargs)
return TetList([ind for ind in inds if ind != UNKNOWN_TET], **self._cloneArgs)
@property
def triNeighbs(self):
"""Neighboring triangles
:type: :py:class:`TriList`, read-only
"""
inds = self.mesh.stepsMesh.getTriTriNeighbs(self._idx, **self._callKwargs)
return TriList(inds, **self._cloneArgs)
@property
def patchTriNeighbs(self):
"""Neighboring triangles that are also part of the same patch
If the triangle is not part of a :py:class:`Patch`, returns an empty :py:class:`TriList`
:type: :py:class:`TriList`, read-only
"""
patch = self.patch
if patch is not None:
inds = self.mesh.stepsMesh.getTriTriNeighb(self._idx, patch.stepsPatch, **self._callKwargs)
inds = [i for i in inds if i != UNKNOWN_TRI]
else:
inds = []
return TriList(inds, **self._cloneArgs)
@property
def patch(self):
"""Patch associated to the triangle (if any)
:type: :py:class:`Patch` or None, read-only
"""
p = self.mesh.stepsMesh.getTriPatch(self._idx, **self._callKwargs)
if p is None:
return None
else:
return getattr(self.mesh, p.getID())
@property
def bars(self):
"""bars of the triangle
:type: :py:class:`VertList`, read-only
"""
return BarList(self.mesh.stepsMesh.getTriBars(self._idx, **self._callKwargs), **self._cloneArgs)
@property
def norm(self):
"""Normal vector of the triangle
:type: :py:class:`Point`, read-only
"""
return Point(*self.mesh.stepsMesh.getTriNorm(self._idx, **self._callKwargs))
def _getPhysicalLocation(self):
"""Return the location associated with the reference (comp for tet, patch for tri)."""
return self.patch
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return TriReference._locStr
class _DistTriReference(TriReference, _DistReference):
@TriReference.tetNeighbs.getter
def tetNeighbs(self):
if self.mesh._local:
inds = self.mesh.stepsMesh.getTriTetNeighb(self._idx, **self._callKwargs, owned=self.mesh._owned)
return TetList([ind for ind in inds if ind != UNKNOWN_TET], **self._cloneArgs)
else:
return TriReference.tetNeighbs.fget(self)
TriReference._distCls = _DistTriReference
[docs]@nutils.FreezeAfterInit
class BarReference(Reference):
"""Convenience class for accessing properties of a bar
:param idx: Index of the bar
:type idx: int
:param mesh: Mesh object that contains the bar
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
"""
_locStr = 'Bar'
def __init__(self, idx, mesh=None, *args, **kwargs):
super().__init__(idx, mesh=mesh, *args, **kwargs)
@property
def verts(self):
"""Vertices of the bar
:type: :py:class:`VertList`, read-only
"""
return VertList(self.mesh.stepsMesh.getBar(self._idx, **self._callKwargs), **self._cloneArgs)
@property
def center(self):
v1, v2 = self.verts
return Point(*((v1 + v2) / 2))
class _DistBarReference(BarReference, _DistReference):
pass
BarReference._distCls = _DistBarReference
[docs]@nutils.FreezeAfterInit
class VertReference(Reference, Point):
"""Convenience class for accessing properties of a vertex
Can be used in the same way as a :py:class:`Point`.
:param idx: Index of the vertex
:type idx: int
:param mesh: Mesh object that contains the vertex
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
.. warning::
When used with MPI, keep in mind that the creation of a :py:class:`VertReference`
retrieves the coordinates of the vertex and thus potentially requires synchronization
across MPI ranks. When using :py:class:`DistMesh`, there is no need for synchronization
if the `local=True` keyword argument is used.
.. note::
In contrast to other references, VertReference can only be created by explicitely
providing a mesh object, it cannot use the context manager syntax described in
:py:class:`TetMesh`.
"""
_locStr = 'Vert'
def __new__(cls, idx, mesh, anonymous=False, *args, _newCalled=False, **kwargs):
if not _newCalled:
return Reference.__new__(cls, idx, mesh, anonymous=anonymous, *args, **kwargs)
if not isinstance(mesh, _BaseTetMesh):
raise TypeError(f'Expected a TetMesh or DistMesh object, got {mesh} instead.')
idx = cls._getActualIdx(idx, mesh)
return Point.__new__(cls, numpy.nan, numpy.nan, numpy.nan)
def __init__(self, idx, mesh=None, *args, _getData=True, **kwargs):
super().__init__(idx, mesh=mesh, *args, **kwargs)
if _getData:
self[0:3] = self.mesh.stepsMesh.getVertex(self._idx, **self._callKwargs)
# Need to redefine __hash__ and __eq__ to be sure the Reference version of these methods
# is called.
def __hash__(self):
return super().__hash__()
def __eq__(self, other):
return super().__eq__(other)
def __ne__(self, other):
return not super().__eq__(other)
def __repr__(self):
return super().__repr__()
def __str__(self):
return super().__repr__()
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return VertReference._locStr
class _DistVertReference(VertReference, _DistReference):
def __init__(self, idx, mesh=None, *args, local=False, _fromLocal=False, **kwargs):
if not local and _fromLocal:
warnings.warn(
f'Global vertex references created from local references or lists are loaded without '
f'coordinates data (which would require global MPI synchronization).'
)
super().__init__(idx, mesh=mesh, *args, local=local, _getData=not _fromLocal, **kwargs)
def toGlobal(self, **_kwargs):
"""Get the global version of the reference
If the reference is not local, it returns itself.
Otherwise, the local index of the original reference is converted to global mesh index and
the corresponding global reference is returned. The returned reference will not contain coordinate
data, which would require an MPI synchronization. In order to access coordinate data, the global
index first needs to be synchronized across ranks.
:rtype: :py:class:`Reference`
"""
return super().toGlobal(_fromLocal=self._local)
VertReference._distCls = _DistVertReference
[docs]class RefList(nutils.UsingObjects(nutils.Optional(_BaseTetMesh)), nutils.SolverPathObject, nutils.Facade):
"""Base convenience class for geometrical element lists
:param lst: The list of elements
:type lst: Iterable[int] or range or Iterable[:py:class:`Reference`] or :py:class:`RefList`
:param mesh: Mesh object that contains the elements
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
Behaves like a python list but with additional functionalities.
"""
class _OptimizationCM:
"""Context manager for handling the updating and invalidation of optimization data
All code that modifies the list needs to be wrapped using ``with self._optimCM:``.
Code that does not modify the list can directly use self._optimCM.
Beware of not wrapping non-modifying code, it can lead to strange behavior.
"""
# TODO Not urgent: Refactor to have a clear distinction (better granularity) between data
# setting linked to actual modification of the list and data setting for local caching
def __init__(self, parent):
self.parent = parent
self._optimData = {}
self._modifData = set()
def __enter__(self):
if self.parent._immutable:
raise Exception(f'Cannot modify an immutable list.')
self._modifData = set()
return self
def __exit__(self, exc_type, exc_val, exc_tb):
self._optimData = {key: val for key, val in self._optimData.items() if key in self._modifData}
def __getitem__(self, key):
if key in self._optimData:
return self._optimData[key]
return None
def __setitem__(self, key, val):
self._optimData[key] = val
self._modifData.add(key)
def __contains__(self, key):
return key in self._optimData
def __new__(cls, lst=None, mesh=None, _immutable=False, _newCalled=False, **kwargs):
# TetList, TriList, etc. are facades to their subclasses, the user always uses the main class
# but the appropriate subclass is used for object instantiation.
if cls not in RefList.__subclasses__() or _newCalled:
return super(RefList, cls).__new__(cls)
if mesh is None:
# Try to infer the mesh from context-managers
mesh = nutils.UsableObject._getUsedObjectOfClass(_BaseTetMesh)
# Try to infer the mesh from the lst argument
if mesh is None:
# If lst is a generator, do not consume it by trying to infer the mesh
lst = None if isinstance(lst, types.GeneratorType) else lst
lst, mesh = cls._getLstAndMeshFromParams(lst, mesh)
if mesh.__class__ == DistMesh:
return cls._distCls.__new__(
cls._distCls, lst, mesh=mesh, _immutable=_immutable, _newCalled=True, **kwargs
)
else:
return cls.__new__(cls, lst, mesh=mesh, _immutable=_immutable, _newCalled=True, **kwargs)
def __init__(self, lst, mesh=None, _immutable=False, *args, **kwargs):
super().__init__(*args, addAsElement=False, **kwargs)
if mesh is None:
(mesh,) = self._getUsedObjects(addAsElement=False)
elif not isinstance(mesh, _BaseTetMesh):
raise TypeError(f'Expected a TetMesh or DistMesh object, got a {type(mesh)} instead.')
self.lst, self.mesh = self._getLstAndMeshFromParams(lst, mesh)
self._immutable = _immutable
# Optimization data wrapper, used as a context manager when modifying the list
self._optimCM = RefList._OptimizationCM(self)
# Keyword dict to initialize a list with the same parameters as this one
self._cloneArgs = dict(mesh=self.mesh)
@classmethod
def _getIdxLst(cls, lst):
lst = list(lst)
if all(isinstance(e, numbers.Integral) for e in lst):
return lst
elif all(isinstance(e, cls._refCls) for e in lst):
return [e._idx for e in lst]
else:
raise TypeError(f'Cannot use {lst} to initialize a list of {cls._refCls}.')
@classmethod
def _toRefList(cls, lst, mesh):
if isinstance(lst, RefList):
if not isinstance(lst, cls):
raise TypeError(f'Expected a {cls}, got a {lst.__class__} instead.')
return lst
else:
return cls(lst, mesh)
@classmethod
def _getLstAndMeshFromParams(cls, lst, mesh):
actualLst = []
actualMesh = mesh
if cls in lst.__class__.__mro__:
actualMesh = lst.mesh
actualLst = copy.copy(lst.lst)
elif isinstance(lst, range):
actualLst = lst
elif hasattr(lst, '__iter__'):
lst = list(lst)
actualLst = cls._getIdxLst(lst)
if len(lst) > 0 and isinstance(lst[0], cls._refCls):
actualMesh = lst[0].mesh
elif lst is not None:
raise TypeError(f'Cannot create an element list from a {type(lst)}')
if actualMesh is None:
raise Exception(f'Unable to identify which mesh should the list be bound to.')
return actualLst, actualMesh
def _splitByLocation(self):
"""Split the list in several lists that contain elements in the same comp / patch
The order of elements is conserved. The name of the list is conserved but can be postfixed.
"""
res = []
loc = (None,)
for elem in self:
newLoc = elem._getPhysicalLocation()
if newLoc != loc:
loc = newLoc
if not self._autoNamed:
name = f'{self.name}({loc})' if len(res) > 0 else self.name
else:
name=None
res.append(self.__class__([], **self._cloneArgs, name=name))
res[-1].append(elem)
return res
[docs] def splitToROIs(self, ROIPrefix, method='grid', gridSize=None):
"""Split an element list into a square grid and define each square as a Region Of Interest
:param ROIPrefix: A prefix for the names of the ROIs, it will be suffixed with 'x_y_z' with x, y, and
z the integer grid coordinate of a given square ROI.
:type ROIPrefix: str
:param method: The method used to define the ROIs, only 'grid' is available for now.
:type method: str
:param gridSize: The side length of a single square from the grid, in meters.
:type gridSize: float
"""
rois = []
if method == 'grid':
gridCells = {}
for elem in self:
key = tuple(elem.center // gridSize)
gridCells.setdefault(key, self.__class__([], **self._cloneArgs)).append(elem)
for key, elems in gridCells.items():
name = ROIPrefix + '_' + '_'.join(map(str, key))
rois.append(ROI(elems, name=name))
return rois
def _modify(self):
"""Return a context manager and signal that the wrapped statements will modify the list"""
# Simply return the _OptimizationCM for now
return self._optimCM
[docs] def __getitem__(self, key):
"""Access one or several element(s) from the list
:param key: Index of the element or slice
:type key: int or slice
:returns: The element(s)
:rtype: :py:class:`Reference` or :py:class:`RefList`
Usage::
>>> lst = RefList(range(10, 20), mesh=mesh)
>>> lst[0] # Returns the 1st element from the list
Tet(10)
>>> lst[0].idx
10
>>> lst2 = lst[1:5] # Returns elements between the 2nd and the 5th (included).
>>> lst2[0]
Tet(11)
.. note::
The index of an element in the list is in general different from its index in the mesh.
:meta public:
"""
idxres = self.lst[key]
if isinstance(idxres, (list, range)):
return self.__class__(idxres, **self._cloneArgs)
else:
return self.__class__._refCls(idxres, **self._cloneArgs)
[docs] def __iter__(self):
"""Get an iterator over the list elements
Usage::
>>> lst = RefList(range(3), mesh=mesh)
>>> for ref in lst:
... print(ref.idx)
...
0
1
2
>>> [ref.idx for ref in lst]
[0, 1, 2]
.. note:
The references that it returns does not have a name, in contrast to a
:py:class:`Reference` created by the user directly.
:meta public:
"""
for idx in self.lst:
yield self.__class__._refCls(idx, **self._cloneArgs, anonymous=True)
[docs] def append(self, e):
"""Append an element to the list
:param e: The element
:type e: :py:class:`Reference`
.. note::
The type of the element depends on the type of the list, one cannot add a
:py:class:`TriReference` to a :py:class:`TetList`.
"""
if not isinstance(e, self.__class__._refCls):
raise TypeError(f'Cannot append {e}, expected a {self.__class__._refCls.__name__} object')
self.lst = self._getLst()
# Modifying the list
with self._modify():
self.lst.append(e._idx)
[docs] def remove(self, e):
"""Remove an element from the list
:param e: The element
:type e: :py:class:`Reference`
"""
if not isinstance(e, self.__class__._refCls):
raise TypeError(f'Cannot remove {e}, expected a {self.__class__._refCls.__name__} object')
self.lst = self._getLst()
# Modifying the list
with self._modify():
self.lst.remove(e._idx)
_ElemsAsSet_K = 'ElemsAsSet'
[docs] def __contains__(self, elem):
"""Check whether the list contains an element
:param elem: The element
:type elem: :py:class:`Reference`
:rtype: bool
Usage::
>>> ref1 = Reference(15, mesh=mesh)
>>> ref2 = Reference(25, mesh=mesh)
>>> lst = RefList(range(20), mesh=mesh)
>>> ref1 in lst
True
>>> ref2 in lst
False
:meta public:
"""
if isinstance(elem, self.__class__._refCls):
# Optimization
if RefList._ElemsAsSet_K not in self._optimCM:
self._optimCM[RefList._ElemsAsSet_K] = set(self.lst)
return elem._idx in self._optimCM[RefList._ElemsAsSet_K]
else:
return False
[docs] def __len__(self):
"""Get the length of the list
:rtype: int
Usage::
>>> lst = RefList(range(5), mesh=mesh)
>>> len(lst)
5
:meta public:
"""
return len(self.lst)
[docs] def index(self, elem):
"""Return the index of the element elem in this list
:param elem: The element
:type elem: :py:class:`Reference`
:returns: The index of the element in the class (not its index in the mesh)
:rtype: int
Behaves like :py:func:`list.index`.
"""
if isinstance(elem, self.__class__._refCls):
if self.mesh is not elem.mesh:
raise Exception(
f'Cannot retrieve the index of an element that is associated to a different mesh.'
)
return self.lst.index(elem._idx)
else:
raise TypeError(f'Expected a {self.__class__._refCls}, got {elem} instead.')
def _checkSameType(self, other):
if other.__class__ != self.__class__:
raise TypeError('Cannot combine lists of different types.')
if self.mesh is not other.mesh:
raise Exception('Cannot combine lists associated to different meshes.')
def _getLst(self):
return list(self.lst) if isinstance(self.lst, range) else self.lst
[docs] def __and__(self, other):
"""Compute the intersection between two lists
:returns: All the elements that are in both lists. Keeps the order of the left operand.
:rtype: :py:class:`RefList`
Usage::
>>> lst1 = RefList(range(0, 10), mesh=mesh)
>>> lst2 = RefList(range(5, 15), mesh=mesh)
>>> lst3 = lst1 & lst2
>>> [ref.idx for ref in lst3]
[5, 6, 7, 8, 9]
:meta public:
"""
self._checkSameType(other)
s2 = set(other.lst)
res = [i for i in self.lst if i in s2]
return self.__class__(res, **self._cloneArgs)
[docs] def __or__(self, other):
"""Compute the union of two lists
:returns: All the elements of the left operand followed by all the elements of the right
operand that are not in the left one.
:rtype: :py:class:`RefList`
Usage::
>>> lst1 = RefList(range(0, 10), mesh=mesh)
>>> lst2 = RefList(range(5, 15), mesh=mesh)
>>> lst3 = lst1 | lst2
>>> [ref.idx for ref in lst3]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ,14]
:meta public:
"""
self._checkSameType(other)
s1 = set(self.lst)
res = self._getLst() + [j for j in other.lst if j not in s1]
return self.__class__(res, **self._cloneArgs)
[docs] def __sub__(self, other):
"""Compute the substraction of one list from another
:returns: All elements that are in the left operand but are absent from the right operand.
:rtype: :py:class:`RefList`
Usage::
>>> lst1 = RefList(range(0, 10), mesh=mesh)
>>> lst2 = RefList(range(5, 15), mesh=mesh)
>>> lst3 = lst1 - lst2
>>> [ref.idx for ref in lst3]
[0, 1, 2, 3, 4]
>>> lst4 = lst2 - lst1
>>> [ref.idx for ref in lst4]
[10, 11, 12, 13, 14]
:meta public:
"""
self._checkSameType(other)
s2 = set(other.lst)
res = [i for i in self.lst if i not in s2]
return self.__class__(res, **self._cloneArgs)
[docs] def __xor__(self, other):
"""Compute the symetric difference of two lists
:returns: All elements in the left operand that are not in the right operand followed
by all elements of the right operand that are not in the left operand.
:rtype: :py:class:`RefList`
Usage::
>>> lst1 = RefList(range(0, 10), mesh=mesh)
>>> lst2 = RefList(range(5, 15), mesh=mesh)
>>> lst3 = lst1 ^ lst2
>>> [ref.idx for ref in lst3]
[0, 1, 2, 3, 4, 10, 11, 12, 13, 14]
>>> lst4 = lst2 ^ lst1
>>> [ref.idx for ref in lst4]
[10, 11, 12, 13, 14, 0, 1, 2, 3, 4]
:meta public:
"""
self._checkSameType(other)
s1, s2 = set(self.lst), set(other.lst)
res = [i for i in self.lst if i not in s2] + [j for j in other.lst if j not in s1]
return self.__class__(res, **self._cloneArgs)
[docs] def __add__(self, other):
"""Concatenate two lists
:returns: All elements of the left operand followed by all elements of the right operand.
:rtype: :py:class:`RefList`
Usage::
>>> lst1 = RefList(range(0, 10), mesh=mesh)
>>> lst2 = RefList(range(5, 15), mesh=mesh)
>>> lst3 = lst1 + lst2
>>> [ref.idx for ref in lst3]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
>>> lst4 = lst2 + lst1
>>> [ref.idx for ref in lst4]
[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
:meta public:
"""
self._checkSameType(other)
return self.__class__(self._getLst() + other._getLst(), **self._cloneArgs)
[docs] def __eq__(self, other):
"""Test for the equality of two lists
:returns: True if the lists have the same elements in the same order, False otherwise
:rtype: bool
Usage::
>>> lst1 = RefList(range(0, 5), mesh=mesh)
>>> lst2 = RefList(range(0, 10), mesh=mesh)
>>> lst1 == lst2
False
>>> lst1 == lst2[0:5]
True
:meta public:
"""
try:
self._checkSameType(other)
return len(self) == len(other) and tuple(self.lst) == tuple(other.lst)
except Exception:
return False
@property
def indices(self):
"""The indices of elements in the list
:type: List[int], read-only
"""
return copy.copy(self.lst)
def __hash__(self):
return hash((len(self), tuple(self.lst)))
def __repr__(self):
if self._autoNamed:
# If the list does not have a given name, represent the list with a hash of its contents
return f'{self.__class__.__name__}#{hash(self)}'
else:
return self.name
def _solverId(self):
"""Return the id that is used to identify an object in the solver."""
return (self._getLst(),)
def _simPathWalkExpand(self):
"""Return an iterable of the elements that should be part of ResultSelector paths."""
return iter(self)
class _DistRefList(RefList):
"""Base class for distributed element lists
Distributed lists can either be constituted of local or global indices
"""
_FACADE_TITLE_STR = 'Only available for lists defined in distributed tetrahedral meshes'
_LocalList = 'LocalList'
_GlobalList = 'GlobalList'
def __init__(self, lst=None, mesh=None, _immutable=False, *args, local=False, _owned=None, **kwargs):
super().__init__(lst, mesh, _immutable, *args, **kwargs)
if isinstance(lst, _DistRefList):
local = lst.isLocal()
self._local = local
self._cloneArgs['local'] = self._local
# The owned flag is just indicative of how the list was created
self._owned = _owned
self._cloneArgs['_owned'] = self._owned
def isLocal(self):
"""Return whether the list is composed of local element indices
:returns: ``True`` if the list is local
:rtype: bool
"""
return self._local
def toLocal(self, owned=True, _returnInds=False):
"""Get the local version of a distributed list
If the list is local, it returns itself. Otherwise, the orginal list is filtered and only
the elements that are hosted by the current MPI process are added to the local list, with
their local index. The returned list is thus different for each MPI process.
:param owned: Whether only owned elements should be considered (defaults to True)
:type owned: bool
:returns: A local version of the list
:rtype: :py:class:`RefList`
"""
# if _returnInds is True, also return the positions of the local elements in the original list
if self._local:
if owned:
# Check that the returned list does not contain any non-owned elements
return self.toGlobal().toLocal(owned=True, _returnInds=_returnInds)
else:
return (self, list(range(len(self)))) if _returnInds else self
cmkey = (_DistRefList._LocalList, owned)
if cmkey not in self._optimCM:
method = getattr(self.mesh.stepsMesh, f'get{self._refCls._locStr}LocalIndex')
localInds = []
listInds = []
for i, gInd in enumerate(self.indices):
localInd = method(gInd, owned=owned)
if localInd is not None:
localInds.append(localInd)
listInds.append(i)
self._optimCM[cmkey] = (
self.__class__(localInds, mesh=self.mesh, local=True, _owned=owned),
listInds,
)
if _returnInds:
return self._optimCM[cmkey]
else:
return self._optimCM[cmkey][0]
def toGlobal(self, **_kwargs):
"""Get the global version of a distributed list
If the list is not local, it returns itself. Otherwise, the local indices of the original
list are converted to global mesh indices and the list of these global indices is returned.
:returns: A global version of the list
:rtype: :py:class:`RefList`
"""
if not self._local:
return self
if _DistRefList._GlobalList not in self._optimCM:
method = getattr(self.mesh.stepsMesh, f'get{self._refCls._locStr}GlobalIndex')
self._optimCM[_DistRefList._GlobalList] = self.__class__(
[ind for ind in map(method, self.indices) if ind is not None],
mesh=self.mesh,
local=False,
**_kwargs
)
return self._optimCM[_DistRefList._GlobalList]
def _checkSameType(self, other):
if self._local != other._local:
raise TypeError('Cannot combine local and global index lists.')
return super()._checkSameType(other)
def combineWithOperator(self, binOp):
"""Combine the local versions of this list across ranks with the given operator
:param binOp: The binary operator used to combine the global lists
:type binOp: Callable[[:py:class:`RefList`, :py:class:`RefList`], :py:class:`RefList`]
:returns: The global element list resulting from the combination of local lists
:rtype: :py:class:`RefList`
"""
if not self._local:
raise Exception(f'The element list is already global, cannot combine it across ranks.')
lists = self.mesh._comm.gather(self.toGlobal().indices, root=0)
if self.mesh._comm.Get_rank() == 0:
allIdx = self.__class__(lists[0], mesh=self.mesh, local=False)
for lst in lists[1:]:
allIdx = binOp(allIdx, self.__class__(lst, mesh=self.mesh, local=False))
else:
allIdx = self.__class__([], mesh=self.mesh, local=False)
lst = self.mesh._comm.bcast(allIdx.indices, root=0)
return self.__class__(lst, mesh=self.mesh, local=False)
def _solverKeywordParams(self):
"""Return the additional keyword parameters that should be passed to the solver"""
return dict(local=self._local)
[docs]@nutils.FreezeAfterInit
class TetList(RefList):
"""Convenience class for tetrahedron lists
:param lst: The list of tetrahedrons
:type lst: Iterable[int] or range or Iterable[:py:class:`TetReference`] or :py:class:`TetList`
:param mesh: Mesh object that contains the elements
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If `lst` is omitted, the list will be initialized to an empty list.
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
Behaves like a :py:class:`RefList` but with additional properties specific to tetrahedrons.
Note that :py:func:`__getitem__` and :py:func:`__contains__` are modified to accept a 3D
point as parameter.
"""
_refCls = TetReference
_dim = 3
def __init__(self, lst=None, mesh=None, _immutable=False, *args, **kwargs):
super().__init__(lst, mesh, _immutable, *args, **kwargs)
[docs] def __getitem__(self, key):
"""Access one or several tetrahedron(s) from the list
:param key: Index of the tetrahedron, or slice, or 3D point
:type key: int or slice or Tuple[float, float, float] or :py:class:`Point`
:returns: The tetrahedron(s)
:rtype: :py:class:`TetReference` or :py:class:`TetList`
Usage::
>>> lst = TetList(range(10, 20), mesh=mesh)
>>> lst[0] # Returns the 1st element from the list
Tet(10)
>>> lst[0].idx
10
>>> lst2 = lst[1:5] # Returns elements between the 2nd and the 5th (included).
>>> lst2[0]
Tet(11)
>>> mesh.tets[0, 0, 0] # Returns the tetrahedron which contains the point x=0, y=0, z=0
Tet(123)
.. note::
The index of an element in the list is in general different from its index in the mesh.
:py:class:`TetList` is the only subclass of :py:class:`RefList` that can take a 3D point
as key of its __getitem__ special method.
:meta public:
"""
try:
return super().__getitem__(key)
except TypeError:
if hasattr(key, '__iter__') and len(key) == 3:
idx = self._findTetIdxByPoint(list(key))
if idx != UNKNOWN_TET and idx in self.lst:
return TetReference(idx, **self._cloneArgs)
else:
raise KeyError(f'No Tetrahedron exists at position {key} in this list.')
else:
raise
[docs] def __contains__(self, key):
"""Check whether the list contains a tetrahedron
:param key: The element or a 3D point
:type key: :py:class:`TetReference` or Tuple[float, float, float] or :py:class:`Point`
:rtype: bool
In addition to the arguments accepted in :py:func:`RefList.__contains__`, it is possible
to check whether a 3D point is inside one tetrahedron in the list.
Usage::
>>> ref1 = TetReference(15, mesh=mesh)
>>> ref2 = TetReference(25, mesh=mesh)
>>> lst = TetList(range(20), mesh=mesh)
>>> ref1 in lst
True
>>> ref2 in lst
False
>>> Point(0, 0, 0) in lst
False
>>> Point(1e-6, 0, 5e-6) in lst
True
:meta public:
"""
if isinstance(key, Reference):
return super().__contains__(key)
elif hasattr(key, '__iter__') and len(key) == 3:
try:
self[key]
return True
except Exception:
return False
return False
_CurrShell_K = 'CurrShell'
[docs] def dilate(self, d=1):
"""Dilates the list to tetrahedron neighbors
One cycle of dilation corresponds to adding all tetrahedrons that are neighbors of the
current tetrahedron list but are not part of the list. This operation is repeated ``d``
times.
:param d: Topological distance to grow, defaults to 1
:type d: int
"""
prevShell = set(self.indices)
if TetList._CurrShell_K not in self._optimCM:
self._optimCM[TetList._CurrShell_K] = TetList(
[tet for tet in self if any(neighb not in self for neighb in tet.neighbs)],
**self._cloneArgs
)
shell = self._optimCM[TetList._CurrShell_K]
for i in range(d):
newShell = set()
for tet in shell:
for neighb in tet.neighbs:
if neighb._idx not in prevShell:
newShell.add(neighb._idx)
prevShell.update(newShell)
newShell = sorted(newShell)
self.lst = self._getLst() + newShell
shell = TetList(newShell, **self._cloneArgs)
with self._modify():
self._optimCM[TetList._CurrShell_K] = shell
self.lst = sorted(set(self.lst))
[docs] def erode(self, d=1):
"""Erodes the list, removing surface tetrahedrons
One cycle of erosion corresponds to removing all tetrahedrons that are on the surface of
the list (i.e. at least one of their four neighbors is not in the list). This operation is
repeated ``d`` times.
:param d: Topological distance to erode, defaults to 1
:type d: int
"""
shell = TetList(
[
tet
for tet in self
if any(neighb not in self for neighb in tet.neighbs) or len(tet.neighbs) < 4
],
**self._cloneArgs,
)
shell.dilate(d - 1)
with self._modify():
self.lst = [tet._idx for tet in self if tet not in shell]
_CurrSurface_K = 'CurrSurface'
_LastSurface_K = 'LastSurface'
@property
def surface(self):
"""Surface of the tetrahedrons in the list
I.e. all triangles that are not shared by two tetrahedrons in the list.
:type: :py:class:`TriList`, read-only
"""
if TetList._CurrSurface_K not in self._optimCM:
if TetList._LastSurface_K in self._optimCM:
# Only compute the changes due to the addition or removal of tetrahedrons
oldLst, oldSurf = self._optimCM[TetList._LastSurface_K]
surf = set(oldSurf)
for tet in TetList(oldLst, **self._cloneArgs) ^ self:
surf ^= set(tet.faces)
else:
# Compute the surface from scratch
surf = set()
for tet in sorted(set(self), key=lambda x:x._idx):
surf ^= set(tet.faces)
# Keep the surface until next list change
self._optimCM[TetList._CurrSurface_K] = TriList(sorted(surf, key=lambda x: x._idx), **self._cloneArgs)
# Keep a version of the current list and surface, for diff computation later
self._optimCM[TetList._LastSurface_K] = (copy.copy(self.lst), self._optimCM[TetList._CurrSurface_K])
return self._optimCM[TetList._CurrSurface_K]
@property
def Vol(self):
"""The summed volume of all tetrahedrons in the list
.. note::
If the list contains duplicate elements, they will be counted only once in the
total volume.
:type: float, read-only
"""
return sum(tet.Vol for tet in sorted(set(self), key=lambda t: t._idx))
@property
def tris(self):
"""All faces of all tetrahedrons in the list
:type: :py:class:`TriList`, read-only
"""
tris = TriList([], **self._cloneArgs)
for tet in self:
tris |= tet.faces
return tris
@property
def bars(self):
"""All bars of all tetrahedrons in the list
:type: :py:class:`BarList`, read-only
"""
return self.tris.bars
@property
def verts(self):
"""All vertices of all tetrahedrons in the list
:type: :py:class:`VertList`, read-only
"""
return self.tris.verts
@nutils.NamedObject.children.getter
def children(self):
try:
allComps = set()
for elem in self:
c = elem.comp
if c is not None:
allComps.add(c)
allChildren = {}
for c in allComps:
allChildren.update(c.children)
return allChildren
except AttributeError as ex:
# Since children is called in __getattr__, raising an AttributeError here can lead to
# infinite recursion.
raise Exception(ex)
def _findTetIdxByPoint(self, pos):
return self.mesh.stepsMesh.findTetByPoint(pos)
class _DistTetList(TetList, _DistRefList):
"""Convenience class for distributed tetrahedron list
"""
@TetList.surface.getter
def surface(self):
if self._local:
return TetList.surface.fget(self)
else:
return self.toLocal().surface.combineWithOperator(operator.xor)
@TetList.tris.getter
def tris(self):
if self._local:
return TetList.tris.fget(self)
else:
return self.toLocal().tris.combineWithOperator(operator.or_)
@TetList.verts.getter
def verts(self):
if self._local:
verts = VertList([], **self._cloneArgs)
for tet in self:
verts |= tet.verts
return verts
else:
return self.toLocal().verts.combineWithOperator(operator.or_)
def _findTetIdxByPoint(self, pos):
return self.mesh.stepsMesh.findTetByPoint(pos, local=self._local)
TetList._distCls = _DistTetList
[docs]@nutils.FreezeAfterInit
class TriList(RefList):
"""Convenience class for triangle lists
:param lst: The list of triangles
:type lst: Iterable[int] or range or Iterable[:py:class:`TriReference`] or :py:class:`TriList`
:param mesh: Mesh object that contains the elements
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If `lst` is omitted, the list will be initialized to an empty list.
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
Behaves like a :py:class:`RefList` but with additional properties specific to triangles.
"""
_refCls = TriReference
_dim = 2
def __init__(self, lst=None, mesh=None, _immutable=False, *args, **kwargs):
super().__init__(lst if lst is not None else [], mesh, _immutable, *args, **kwargs)
@property
def edges(self):
"""Perimeter of the triangles in the list
I.e. all bars that are not shared by two triangles in the list.
:type: :py:class:`BarList`, read-only
"""
edg = set()
for tri in self:
edg ^= set(tri.bars)
return BarList(sorted(edg, key=lambda x: x._idx), **self._cloneArgs)
@property
def bars(self):
"""All bars of all triangles in the list
:type: :py:class:`BarList`, read-only
"""
bars = BarList([], **self._cloneArgs)
for tri in self:
bars |= tri.bars
return bars
@property
def verts(self):
"""All vertices of all triangles in the list
:type: :py:class:`VertList`, read-only
"""
return self.bars.verts
@property
def Area(self):
"""The summed area of all triangles in the list
.. note::
If the list contains duplicate elements, they will be counted only once in the
total area.
:type: float, read-only
"""
return sum(tri.Area for tri in sorted(set(self), key=lambda t: t._idx))
@nutils.NamedObject.children.getter
def children(self):
try:
allPatches = set()
for elem in self:
p = elem.patch
if p is not None:
allPatches.add(p)
allChildren = {}
for p in allPatches:
allChildren.update(p.children)
return allChildren
except AttributeError as ex:
# Since children is called in __getattr__, raising an AttributeError here can lead to
# infinite recursion.
raise Exception(ex)
class _DistTriList(TriList, _DistRefList):
@TetList.verts.getter
def verts(self):
if self._local:
verts = VertList([], **self._cloneArgs)
for tri in self:
verts |= tri.verts
return verts
else:
return self.toLocal().verts.combineWithOperator(operator.or_)
TriList._distCls = _DistTriList
[docs]@nutils.FreezeAfterInit
class BarList(RefList):
"""Convenience class for bar lists
:param lst: The list of bars
:type lst: Iterable[int] or range or Iterable[:py:class:`BarReference`] or :py:class:BarList`
:param mesh: Mesh object that contains the elements
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If `lst` is omitted, the list will be initialized to an empty list.
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
Behaves like a :py:class:`RefList` but with additional properties specific to bars.
"""
_refCls = BarReference
_dim = 1
def __init__(self, lst=None, mesh=None, _immutable=False, *args, **kwargs):
super().__init__(lst if lst is not None else [], mesh, _immutable, *args, **kwargs)
@property
def verts(self):
"""All vertices of all bars in the list
:type: :py:class:`VertList`, read-only
"""
verts = VertList([], **self._cloneArgs)
for bar in self:
verts |= bar.verts
return verts
class _DistBarList(BarList, _DistRefList):
pass
BarList._distCls = _DistBarList
[docs]@nutils.FreezeAfterInit
class VertList(RefList):
"""Convenience class for vertex lists
:param lst: The list of vertices
:type lst: Iterable[int] or range or Iterable[:py:class:`VertReference`] or :py:class:VertList`
:param mesh: Mesh object that contains the elements
:type mesh: Union[:py:class:`TetMesh`, :py:class:`DistMesh`, None]
If `lst` is omitted, the list will be initialized to an empty list.
If mesh is not provided, it will be inferred, if possible, from context managers
(see :py:class:`TetMesh`).
Behaves like a :py:class:`RefList`.
"""
_refCls = VertReference
_dim = 0
def __init__(self, lst=None, mesh=None, _immutable=False, *args, **kwargs):
super().__init__(lst if lst is not None else [], mesh, _immutable, *args, **kwargs)
class _DistVertList(VertList, _DistRefList):
def __init__(self, *args, _fromLocal=False, **kwargs):
super().__init__(*args, **kwargs)
self._cloneArgs['_fromLocal'] = _fromLocal
VertList._distCls = _DistVertList
# Add a link from the reference classes to the list class
for cls in RefList.__subclasses__():
if cls != _DistRefList:
cls._refCls._lstCls = cls
# Do not leave cls as a global variable
del cls
[docs]@nutils.FreezeAfterInit
class ROI(nutils.UsingObjects(_BaseTetMesh), nutils.StepsWrapperObject):
"""Region of Interest class
:param lst: Element list
:type lst: :py:class:`TetList`, or :py:class:`TriList`, or :py:class:`VertList`
Should be declared using the context manager syntax with :py:func:`TetMesh`::
with mesh:
tetROI1 = ROI.Create(mesh.tets[10:30])
The ROI is then treated as a physical location, like patches or compartment,
and can be accessed from the mesh with its name::
mesh.tetROI1
"""
_locStr = 'ROI'
def __init__(self, lst=None, *args, _createObj=True, **kwargs):
super().__init__(*args, **kwargs)
(mesh,) = self._getUsedObjects()
self.mesh = mesh
if _createObj:
self.stepsROI, self._elemType = self._createStepsObj(mesh, lst)
else:
self.stepsROI, self._elemType = None, None
if self._elemType is not None and lst is not None:
self.lst = self._elemType._lstCls._toRefList(lst, mesh)
else:
self.lst = None
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsROI]
@classmethod
def _FromStepsObject(cls, obj, mesh, name):
"""Create the interface object from a STEPS object."""
revTypeMap = {
ELEM_TET: TetReference,
ELEM_TRI: TriReference,
ELEM_VERTEX: VertReference,
}
region = cls(_createObj=False, name=name)
region.stepsROI = obj
region._elemType = revTypeMap[obj.type]
region.lst = region._elemType._lstCls._toRefList(obj.indices, mesh)
return region
def _createStepsObj(self, mesh, lst):
typeMap = {
TetReference: ELEM_TET,
TriReference: ELEM_TRI,
VertReference: ELEM_VERTEX,
}
if len(lst) == 0:
raise Exception('An empty list was given.')
elem = next(iter(lst))
tpe = None
for cls, _tpe in typeMap.items():
if isinstance(elem, cls):
tpe = _tpe
if tpe is None:
raise TypeError('Unsupported element type for creating an ROI.')
mesh.stepsMesh.addROI(self.name, tpe, lst.indices)
return mesh.stepsMesh.getROI(self.name), elem.__class__
def __getitem__(self, key):
"""
Indexing can only be used with ':' or '...' to signify that we want to consider all
elements separately instead of the full ROI.
"""
if key == slice(None) or key is Ellipsis:
return self.lst
else:
raise KeyError(
f'ROIs can only be indexed with an empty slice ":" or the ellipsis operator "...".'
)
@property
def Vol(self):
"""Total volume of the ROI, if applicable
:type: float, read-only
"""
return self.mesh.stepsMesh.getROIVol(self.name)
@property
def Area(self):
"""Total area of the ROI, if applicable
:type: float, read-only
"""
return self.mesh.stepsMesh.getROIArea(self.name)
def _checkType(self, tpe):
if self._elemType != tpe:
raise AttributeError(f'ROI {self} is not composed of {tpe}.')
@property
def tets(self):
"""All tetrahedrons in the ROI
Raises a :py:class:`TypeError` if the ROI was not created from a :py:class:`TetList`.
:type: :py:class:`TetList`, read-only
"""
self._checkType(TetReference)
return self.lst
@property
def tris(self):
"""All triangles in the ROI
Raises a :py:class:`TypeError` if the ROI was not created from a :py:class:`TriList`.
:type: :py:class:`TriList`, read-only
"""
self._checkType(TriReference)
return self.lst
@property
def verts(self):
"""All vertices in the ROI
Raises a :py:class:`TypeError` if the ROI was not created from a :py:class:`VertList`.
:type: :py:class:`VertList`, read-only
"""
self._checkType(VertReference)
return self.lst
@property
def children(self):
# Handle the references to species etc.
for elem in self.lst:
try:
return elem._getPhysicalLocation().children
except:
continue
return {}
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return ROI._locStr
def _simPathAutoMetaData(self):
"""Return a dictionary with string keys and string or numbers values."""
return {'loc_type': self.__class__._locStr, 'loc_id': self.name}
[docs]@nutils.FreezeAfterInit
class DiffBoundary(nutils.UsingObjects(_BaseTetMesh), nutils.StepsWrapperObject):
"""Diffusion Boundary class
Annotation of a group of triangles in a Tetmesh. The triangles form a boundary between two
compartments, that may allow diffusion of some specified species.
:param lst: List of triangles forming the boundary
:type lst: :py:class:`TriList`
"""
_locStr = 'DiffBoundary'
_lstType = TriList
_direcStr = 'direction_comp'
def __init__(self, lst, *args, **kwargs):
super().__init__(*args, **kwargs)
(mesh,) = self._getUsedObjects()
self.lst = self.__class__._lstType(lst, mesh)
self.stepsDiffBound = self._createStepsObj(mesh)
self._direc = None
def _getStepsObjects(self):
"""Return a list of the steps objects that this named object holds."""
return [self.stepsDiffBound]
def _createStepsObj(self, mesh):
if isinstance(mesh, TetMesh):
return stepslib._py_DiffBoundary(self.name, mesh.stepsMesh, [tri._idx for tri in self.lst])
elif isinstance(mesh, DistMesh):
tet1, tet2 = self.lst[0].tetNeighbs
comp1 = tet1.comp.name
comp2 = tet2.comp.name
mesh._getStepsObjects()[0].addDiffusionBoundary(self.name, comp1, comp2, self.lst.indices)
return None
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return DiffBoundary._locStr
def _solverKeywordParams(self):
"""Return the additional keyword parameters that should be passed to the solver"""
if self._direc is not None:
return {DiffBoundary._direcStr: self._direc._solverId()[0]}
else:
return {}
[docs] def __call__(self, direc=None):
"""Get a version of the diffusion boundary with added information
Parentheses notation (function call) can be used on a diffusion boundary to represent a
diffusion boundary further specified with additional information. This should never be needed
during model declaration but can become handy for simulation control and data saving
(see :py:class:`steps.API_2.sim.SimPath`).
:param direc: A direction for the diffusion boundary (i.e. a target compartment for
diffusion boundaries and a target patch for surface diffusion boundary)
:type key: :py:class:`Comp`
:returns: A version of the diffusion boundary that is specific to direction ``direc``.
:meta public:
"""
db = copy.copy(self)
db._direc = direc
return db
@property
def tris(self):
"""All triangles in the boundary, if applicable
:type: :py:class:`TriList`, read-only
"""
if not isinstance(self.lst, self._lstType):
raise NotImplementedError()
return self.lst
@property
def bars(self):
""":meta private:"""
raise NotImplementedError()
@property
def children(self):
# Handle the references to species.
for tet in self.lst[0].tetNeighbs:
try:
return tet.comp.children
except:
continue
return {}
[docs]@nutils.FreezeAfterInit
class SDiffBoundary(DiffBoundary):
"""Surface Diffusion Boundary class
Annotation of a group of bars in a Tetmesh. The bars form a boundary between two patches,
that may allow diffusion of some specified species. The patches to be connected by this
surface diffusion boundary need to be specified to avoid potential ambiguity.
:param lst: List of bars forming the boundary
:type lst: :py:class:`BarList`
:param patch1: One of the patches connected to the boundary
:type patch1: :py:class:`Patch`
:param patch2: The other patch connected to the boundary
:type patch2: :py:class:`Patch`
"""
_locStr = 'SDiffBoundary'
_lstType = BarList
_direcStr = 'direction_patch'
def __init__(self, lst, patch1, patch2, *args, **kwargs):
self.patches = [patch1, patch2]
super().__init__(lst, *args, **kwargs)
def _createStepsObj(self, mesh):
return stepslib._py_SDiffBoundary(
self.name, mesh.stepsMesh, [bar._idx for bar in self.lst], [p.stepsPatch for p in self.patches]
)
def _solverStr(self):
"""Return the string that is used as part of method names for this specific object."""
return SDiffBoundary._locStr
def _solverKeywordParams(self):
"""Return the additional keyword parameters that should be passed to the solver"""
if self._direc is not None:
return {SDiffBoundary._direcStr: self._direc._solverId()[0]}
else:
return {}
@property
def tris(self):
""":meta private:"""
raise NotImplementedError()
@property
def bars(self):
"""All bars in the boundary, if applicable
:type: :py:class:`BarList`, read-only
"""
if not isinstance(self.lst, self._lstType):
raise NotImplementedError()
return self.lst
@property
def children(self):
# Handle the references to species.
for patch in self.patches:
try:
return patch.children
except:
continue
return {}
###################################################################################################
# Mesh partitioning
[docs]@nutils.FreezeAfterInit
class MeshPartition(nutils.Params):
"""Mesh element partition for parallel simulations
Wrapper class that groups tetrahedrons, triangles and well-mixed partitions.
:param mesh: Mesh being partitioned
:type mesh: :py:class:`TetMesh`
:param tet_hosts: List of host process indices for each tetrahedron
:type tet_hosts: List[int]
:param tri_hosts: Mapping between triangle index and host process index
:type tri_hosts: Mapping[int, int]
:param wm_hosts: List of host process indices for each well mixed compartment
:type wm_hosts: List[int]
"""
def __init__(self, mesh, tet_hosts=[], tri_hosts={}, wm_hosts=[]):
super().__init__(tet_hosts=tet_hosts, tri_hosts=tri_hosts, wm_hosts=wm_hosts)
self._mesh = mesh
self._tet_hosts = tet_hosts
self._tri_hosts = tri_hosts
self._wm_hosts = wm_hosts
@property
def tetPart(self):
"""The partition data for tetrahedrons
:type: List[int], read-only
"""
return self._tet_hosts
@property
def triPart(self):
"""The partition data for triangles
:type: Dict[int, int], read-only
"""
return self._tri_hosts
@property
def wmPart(self):
"""The partition data for well mixed locations
:type: List[int], read-only
"""
return self._wm_hosts
@property
def tetTable(self):
"""A dictionary in the format of {host0: TetList([tet0, tet1, ...]), ...}
:type: Dict[int, :py:class:`TetList`], read-only
"""
return {
h: TetList(lst, mesh=self._mesh)
for h, lst in sgdecomp.getTetPartitionTable(self._tet_hosts).items()
}
@property
def triTable(self):
"""A dictionary in the format of {host0: TriList([tri0, tri1, ...]), ...}
:type: Dict[int, :py:class:`TriList`], read-only
"""
return {
h: TriList(lst, mesh=self._mesh)
for h, lst in sgdecomp.getTriPartitionTable(self._tri_hosts).items()
}
[docs] def validate(self):
"""Validate the partitioning of the mesh"""
sgdecomp.validatePartition(self._mesh.stepsMesh, self._tet_hosts, self._tri_hosts)
[docs] def printStats(self):
"""Print out partitioning stastics
:returns: ``(tet_stats, tri_stats, wm_stats, num_hosts, min_degree, max_degree, mean_degree)``
[tet/tri/wm]_stats contains the number of tetrahedrons/triangles/well-mixed volumes in
each hosting process, num_hosts provide the number of hosting processes,
[min/max/mean]_degree provides the minimum/maximum/average connectivity degree of the
partitioning.
"""
return sgdecomp.printPartitionStat(
self._tet_hosts, self._tri_hosts, self._wm_hosts, self._mesh.stepsMesh
)
def _getTriPartitionFromTet(mesh, tet_hosts, default_tris=None):
"""Return the tri partition corresponding to the tet partition given as an argument."""
triInds = set(default_tris.indices) if default_tris is not None else set()
for patch in mesh._getChildrenOfType(Patch):
triInds |= set(patch.tris.indices)
return sgdecomp.partitionTris(mesh.stepsMesh, tet_hosts, sorted(triInds))
[docs]def LinearMeshPartition(mesh, xbin, ybin, zbin, default_tris=None):
"""Partition the mesh by linearly binning along the 3 axes
First partition the tetrahedrons and then compute a triangle partition that
matches the tetrahedron one.
:param mesh: The mesh to be partitioned
:type mesh: :py:class:`TetMesh`
:param xbin: Number of bins on the x axis
:type xbin: int
:param ybin: Number of bins on the y axis
:type ybin: int
:param zbin: Number of bins on the z axis
:type zbin: int
:param default_tris: Optional list of triangles that should be partitioned even if they are
not part of any patch
:type default_tris: :py:class:`TriList`
:returns: The partition object
:rtype: :py:class:`MeshPartition`
.. note::
The triangle partition only takes into account triangles that are part of a
:py:class:`Patch`.
"""
tet_hosts = sgdecomp.linearPartition(mesh.stepsMesh, (xbin, ybin, zbin))
tri_hosts = _getTriPartitionFromTet(mesh, tet_hosts, default_tris)
return MeshPartition(mesh, tet_hosts=tet_hosts, tri_hosts=tri_hosts)
[docs]def MetisPartition(mesh, path, default_tris=None):
"""Partition the mesh using a Metis .epart file
First partition the tetrahedrons according to the Metis file and then compute a triangle
partition that matches the tetrahedron one.
:param mesh: The mesh to be partitioned
:type mesh: :py:class:`TetMesh`
:param path: Path to the .epart file generated by Metis
:type path: str
:param default_tris: Optional list of triangles that should be partitioned even if they are
not part of any patch
:type default_tris: :py:class:`TriList`
:returns: The partition object
:rtype: :py:class:`MeshPartition`
.. note::
The triangle partition only takes into account triangles that are part of a
:py:class:`Patch`.
"""
tet_hosts = smetis.readPartition(path)
tri_hosts = _getTriPartitionFromTet(mesh, tet_hosts, default_tris)
return MeshPartition(mesh, tet_hosts=tet_hosts, tri_hosts=tri_hosts)
[docs]def GmshPartition(mesh, default_tris=None):
"""Retrive partition information stored in a pre-partitioned gmsh.
The mesh needs to be pre-partitioned and stored in a single
Gmsh file with format version 2.2 ascii.
:param mesh: The mesh to be partitioned
:type mesh: :py:class:`TetMesh`
:param default_tris: Optional list of triangles that should be partitioned even if they are
not part of any patch
:type default_tris: :py:class:`TriList`
:returns: The partition object
:rtype: :py:class:`MeshPartition`
.. note::
The triangle partition only takes into account triangles that are part of a
:py:class:`Patch`.
"""
tet_hosts = [None] * len(mesh.tets)
for k, v in mesh.tetGroups.items():
# partitioning tag is the forth one in gmsh v2.2 format
# starting from 0
if k[0] == 3:
# gmsh partition indices start from 1
# mpi ranks start from 0
partition_rank = k[1] - 1
partition_tets = v
for tet_ref in partition_tets:
tet_hosts[tet_ref._idx] = partition_rank
if None in tet_hosts:
raise Exception("Partition information not available for one or more tets.")
tri_hosts = _getTriPartitionFromTet(mesh, tet_hosts, default_tris)
return MeshPartition(mesh, tet_hosts=tet_hosts, tri_hosts=tri_hosts)
[docs]def MorphPartition(mesh, morph, scale=1e-6, default_tris=None):
"""Partition the mesh using morphological sectioning data
First partition the tetrahedrons according to the morph data and then compute a triangle
partition that matches the tetrahedron one.
:param mesh: The mesh to be partitioned
:type mesh: :py:class:`TetMesh`
:param morph: Morphological sectioning data
:type morph: :py:class:`Morph`
:param scale: Scaling factor from morphological sectioning data to Tetmesh data,
the default value is 1e-6, that is 1 unit (usually micron) in the morph file equals
1e-6 unit (meter) in Tetmesh
:type scale: float
:param default_tris: Optional list of triangles that should be partitioned even if they are
not part of any patch
:type default_tris: :py:class:`TriList`
:returns: The partition object and a mapping between partition index and morphological
section name
:rtype: Tuple[:py:class:`MeshPartition`, Dict[int, str]]
"""
tet_map = smorph.mapMorphTetmesh(morph._sections, mesh.stepsMesh, scale)
tet_hosts = []
ind2sec = {}
sec2ind = {}
i = 0
for name in tet_map:
if name not in sec2ind:
sec2ind[name] = i
ind2sec[i] = name
i += 1
tet_hosts.append(sec2ind[name])
tri_hosts = _getTriPartitionFromTet(mesh, tet_hosts, default_tris)
return MeshPartition(mesh, tet_hosts=tet_hosts, tri_hosts=tri_hosts), ind2sec
###################################################################################################
# Morph sectioning
[docs]@nutils.FreezeAfterInit
class Morph:
"""Morphological sectioning data
Contains data describing the sections of a NEURON morphology.
This class should not be instantiated directly, it should be created using one of the dedicated
class methods.
"""
def __init__(self, sections):
self._sections = sections
[docs] @classmethod
def Load(self, path):
"""Load morphological data from a morph file
:param path: Path to the file
:type path: str
:returns: The morphological sectioning data
:rtype: :py:class:`Morph`
"""
with open(path, 'rb') as f:
return Morph(pickle.load(f))
[docs] @classmethod
def LoadHOC(self, path):
"""Load morphological data from a NEURON hoc file
:param path: Path to the file
:type path: str
:returns: The morphological sectioning data
:rtype: :py:class:`Morph`
"""
return Morph(smorph.hoc2morph(path))
[docs] @classmethod
def LoadSWC(self, path):
"""Load morphological data from an swc file
:param path: Path to the file
:type path: str
:returns: The morphological sectioning data
:rtype: :py:class:`Morph`
"""
return Morph(smorph.swc2morph(path))
[docs] def Save(self, path):
"""Save to a morph file
:param path: Path to the file
:type path: str
"""
with open(path, 'wb') as f:
pickle.dump(self._sections, f)