diff --git a/docs/_images/lung4_human1_coarse.jpg b/docs/_images/lung4_human1_coarse.jpg
new file mode 100644
index 00000000..862b8800
Binary files /dev/null and b/docs/_images/lung4_human1_coarse.jpg differ
diff --git a/docs/_images/lung4_left_right_fitted.jpg b/docs/_images/lung4_left_right_fitted.jpg
new file mode 100644
index 00000000..d47f53aa
Binary files /dev/null and b/docs/_images/lung4_left_right_fitted.jpg differ
diff --git a/docs/scaffolds/lung.rst b/docs/scaffolds/lung.rst
index 89df71b6..6d56ee9f 100644
--- a/docs/scaffolds/lung.rst
+++ b/docs/scaffolds/lung.rst
@@ -1,148 +1,12 @@
Lung Scaffold
=============
-The current recommended lung scaffold is ``3D Lung 2`` built from ``class MeshType_3d_Lung2``;
-the human variant is shown in :numref:`fig-scaffoldmaker-human-lung`.
+For human lung studies the most recent :doc:`lung4` with open fissures is recommended.
+It supports higher, variable resolution meshes.
-.. _fig-scaffoldmaker-human-lung:
+For other species including mouse, rat and pig with accessory lobes and either open or closed fissures, :doc:`lung2` can be used.
-.. figure:: ../_images/scaffoldmaker_human_lung.jpg
- :align: center
+.. toctree::
- Human lung scaffold.
-
-The lung scaffold is a 3-D volumetric model of the lungs representing left lung, right lung and its lobes.
-Depending on the species the left lung may have lower and upper lobes, while the right lung has lower, middle and
-upper lobes, and non-human right lungs usually have a 4th accessory lobe (also known as the diaphragmatic lobe).
-
-The lung scaffold is only a representation of the volumetric spaces of the lobes and does not include representations
-of the pulmonary airway, blood vessels or alveoli.
-
-.. note::
-
- The lung scaffold contains two (or more) independent meshes for the left lung, right lung, and accessory lobes, and
- optionally has open fissures so the lower, middle and upper lobes of each side can be independent meshes.
-
-Variants
---------
-
-The lung scaffold is provided with parameter sets for the following four species, which differ in shape, and in
-particular have different numbers of lobes:
-
-* Human (2 lobes in the left, 3 lobes in the right lung)
-* Mouse (1 lobe in the left, 4 lobes in the right lung)
-* Pig (2 lobes in the left, 4 lobes in the right lung)
-* Rat (1 lobe in the left, 4 lobes in the right lung)
-
-These variants' geometry and annotations are best viewed in the **Scaffold Creator** tool in the ABI Mapping Tools.
-On the web, the latest published generic lung scaffold variants can be viewed on the
-`SPARC Portal `_ by searching for ``lung``, filtering for anatomical models, selecting a variant
-and viewing the scaffold in its Gallery tab or via the `Organ Scaffolds
-`_ help article.
-
-The lung scaffold script generates the scaffold mesh and geometry from an idealization of their shapes. The left and
-right lung (excluding accessory lobe) are generated as half ellipsoids which are then reshaped by smooth functions for
-which parameters are provided on the scaffold, to give approximately realistic geometry for the species.
-The accessory lobe is similarly created as a triangular prism and reshaped.
-
-The generic lung scaffolds are parameterized and fitted to segmentation data from CT and MRI images from the following
-sources: human (`Osanlouy et al. `_), mouse
-(`Beichel et al. `_), pig
-(`Lee et al. `_), rat
-(`NeuroRat V4.0 `_).
-
-[A special ``Material`` parameter set is provided to allow new species' parameters to be developed from the material
-coordinates definition (see below).]
-These parameters were carefully tuned for each species, and it is not recommended that these be edited.
-
-An advanced optional feature is to check *Open fissures* (set parameter to ``true``) which separates the lobes into
-independent meshes allowing elements on opposite sides of each fissure to move independently.
-
-Coordinates
------------
-
-The lung scaffold defines both geometric and material coordinates.
-
-The geometric ``coordinates`` field gives an approximate, idealized unit-scale representation of the lung shape for the
-species, which is intended to be fitted to actual data for a specimen.
-
-The material coordinates field ``lung coordinates`` defines a highly idealized coordinate system to give
-permanent locations for embedding structures in the lungs, defined as 2 half ellipsoids for the left and right lung
-(excluding accessory lobe) and a triangular wedge for the accesory lobe, if present. These can be viewed by visualising
-this field in the *Display* tab of **Scaffold Creator** or by switching to the special ``Material`` parameter set.
-
-The lung scaffold supports limited refinement/resampling by checking *Refine* (set parameter to ``true``) with chosen
-*Refine number of elements* parameter. Be aware that only the ``coordinates`` field is currently defined on the refined
-mesh (but annotations are transferred).
-
-Annotations
------------
-
-Important anatomical regions of the lungs are defined by groups of elements (or faces, edges and nodes/points) and
-annotated with standard term names and identifiers from a controlled vocabulary.
-
-Annotated 3-dimensional volume regions are defined by groups of 3-D elements including (using only one of the items
-separated by slash /):
-
-* left/right lung
-* lower/middle/upper lobe of left/right lung
-* lung
-* right lung accessory lobe
-
-**Terms for volume regions such as the above are not to be used for digitized contours!** They are used for applying
-different material properties in models and the strain/curvature penalty (stiffness) parameters in fitting.
-
-Annotated 2-dimensional surface regions are defined for matching annotated contours digitized from medical images
-including (where ``surface`` is the outside boundary on the meshes and using only one of the items separated by slash
-/):
-
-* base of left/right lung surface
-* base of lower lobe of left/right lung surface
-* base of middle lobe of right lung surface
-* base of upper lobe of left lung surface
-* horizontal fissure of left/right lung
-* horizontal fissure of lower/middle/upper lobe of left/right lung
-* lateral/medial surface of left/right lung
-* lateral/medial surface of lower/middle/upper lobe of left/right lung
-* left/right lung surface
-* lower/middle/upper lobe of left/right lung surface
-* oblique fissure of left/right lung
-* oblique fissure of lower/middle/upper lobe of left/right lung
-* right lung accessory lobe surface
-* base of right lung accessory lobe surface
-* left/right/dorsal/ventral surface of right lung accessory lobe
-
-Annotated 1-dimensional line regions are defined for matching annotated contours digitized from medical images including
-(using only one of the items separated by slash /):
-
-* anterior border of left/right lung
-
-Several fiducial marker points are defined on the lung scaffold, of which the followings are potentially usable when
-digitizing:
-
-* apex of left/right lung
-* laterodorsal tip of middle lobe of right lung
-* medial/ventral base of left/right lung
-* dorsal/ventral apex of right lung accessory lobe
-* left/right dorsal/ventral base of right lung accessory lobe
-
-**Digitization tips to assist fitting:**
-
-1. A proper lung model requires accurate location of all surfaces of the lobes. This requires digitizing the fissures
-and exterior surfaces of the lungs. It's not a requirement to use the most specific annotation group for a surface
-(e.g. for a particular lobe instead of the whole left/right lung) but it may make fitting more efficient.
-
-2. Digitize any fiducial markers you can identify as these are gold standard locations which can be highly weighted in
-the fit.
-
-3. The sharp anterior edges of the left/right lungs are difficult to fit. At a minimum it's important to have data
-marker points for the ``ventral base of left/right lung`` fiducial markers. To properly fit the rest of the edges it's
-best to have data points/contours along these edges annotated with the 1-D ``anterior border of left/right lung``
-terms, which can be weighted highly in the fit. Annotating with lateral/medial surfaces may work in some cases, but just
-using lung/lobe surface groups can be problematic; in both cases there may not be enough data to definitively pull the
-edge into position during the fit.
-
-4. For fitting a lung scaffold with open fissures, annotate fissure data points/contours which clearly belong to only
-one lobe with the fissure term specific to that lobe. Where fissures are too close to distinguish the lobe they are on,
-annotate digitized data the generic fissure term (not for a specific lobe). Doing this allows the same data to fit all
-lobes correctly.
+ lung2
+ lung4
diff --git a/docs/scaffolds/lung2.rst b/docs/scaffolds/lung2.rst
new file mode 100644
index 00000000..f6f9c68e
--- /dev/null
+++ b/docs/scaffolds/lung2.rst
@@ -0,0 +1,150 @@
+Lung 2 Scaffold
+===============
+
+The ``3D Lung 2`` Scaffold built from ``class MeshType_3d_lung2`` has variants for human, mouse, rat and pig, the non-human models having an accessory lobe.
+Open fissures are optionally supported.
+
+The human variant which is shown in :numref:`fig-scaffoldmaker-human-lung2`, but note :doc:`lung4` is now preferred for human lung studies.
+
+.. _fig-scaffoldmaker-human-lung2:
+
+.. figure:: ../_images/scaffoldmaker_human_lung.jpg
+ :align: center
+
+ Human lung scaffold.
+
+The lung scaffold is a 3-D volumetric model of the lungs representing left lung, right lung and its lobes.
+Depending on the species the left lung may have lower and upper lobes, while the right lung has lower, middle and
+upper lobes, and non-human right lungs usually have a 4th accessory lobe (also known as the diaphragmatic lobe).
+
+The lung scaffold is only a representation of the volumetric spaces of the lobes and does not include representations
+of the pulmonary airway, blood vessels or alveoli.
+
+.. note::
+
+ The lung scaffold contains two (or more) independent meshes for the left lung, right lung, and accessory lobes, and
+ optionally has open fissures so the lower, middle and upper lobes of each side can be independent meshes.
+
+Variants
+--------
+
+The lung scaffold is provided with parameter sets for the following four species, which differ in shape, and in
+particular have different numbers of lobes:
+
+* Human (2 lobes in the left, 3 lobes in the right lung)
+* Mouse (1 lobe in the left, 4 lobes in the right lung)
+* Pig (2 lobes in the left, 4 lobes in the right lung)
+* Rat (1 lobe in the left, 4 lobes in the right lung)
+
+These variants' geometry and annotations are best viewed in the **Scaffold Creator** tool in the ABI Mapping Tools.
+On the web, the latest published generic lung scaffold variants can be viewed on the
+`SPARC Portal `_ by searching for ``lung``, filtering for anatomical models, selecting a variant
+and viewing the scaffold in its Gallery tab or via the `Organ Scaffolds
+`_ help article.
+
+The lung scaffold script generates the scaffold mesh and geometry from an idealization of their shapes. The left and
+right lung (excluding accessory lobe) are generated as half ellipsoids which are then reshaped by smooth functions for
+which parameters are provided on the scaffold, to give approximately realistic geometry for the species.
+The accessory lobe is similarly created as a triangular prism and reshaped.
+
+The generic lung scaffolds are parameterized and fitted to segmentation data from CT and MRI images from the following
+sources: human (`Osanlouy et al. `_), mouse
+(`Beichel et al. `_), pig
+(`Lee et al. `_), rat
+(`NeuroRat V4.0 `_).
+
+[A special ``Material`` parameter set is provided to allow new species' parameters to be developed from the material
+coordinates definition (see below).]
+These parameters were carefully tuned for each species, and it is not recommended that these be edited.
+
+An advanced optional feature is to check *Open fissures* (set parameter to ``true``) which separates the lobes into
+independent meshes allowing elements on opposite sides of each fissure to move independently.
+
+Coordinates
+-----------
+
+The lung scaffold defines both geometric and material coordinates.
+
+The geometric ``coordinates`` field gives an approximate, idealized unit-scale representation of the lung shape for the
+species, which is intended to be fitted to actual data for a specimen.
+
+The material coordinates field ``lung coordinates`` defines a highly idealized coordinate system to give
+permanent locations for embedding structures in the lungs, defined as 2 half ellipsoids for the left and right lung
+(excluding accessory lobe) and a triangular wedge for the accesory lobe, if present. These can be viewed by visualising
+this field in the *Display* tab of **Scaffold Creator** or by switching to the special ``Material`` parameter set.
+
+The lung scaffold supports limited refinement/resampling by checking *Refine* (set parameter to ``true``) with chosen
+*Refine number of elements* parameter. Be aware that only the ``coordinates`` field is currently defined on the refined
+mesh (but annotations are transferred).
+
+Annotations
+-----------
+
+Important anatomical regions of the lungs are defined by groups of elements (or faces, edges and nodes/points) and
+annotated with standard term names and identifiers from a controlled vocabulary.
+
+Annotated 3-dimensional volume regions are defined by groups of 3-D elements including (using only one of the items
+separated by slash /):
+
+* left/right lung
+* lower/middle/upper lobe of left/right lung
+* lung
+* right lung accessory lobe
+
+**Terms for volume regions such as the above are not to be used for digitized contours!** They are used for applying
+different material properties in models and the strain/curvature penalty (stiffness) parameters in fitting.
+
+Annotated 2-dimensional surface regions are defined for matching annotated contours digitized from medical images
+including (where ``surface`` is the outside boundary on the meshes and using only one of the items separated by slash
+/):
+
+* base of left/right lung surface
+* base of lower lobe of left/right lung surface
+* base of middle lobe of right lung surface
+* base of upper lobe of left lung surface
+* horizontal fissure of left/right lung
+* horizontal fissure of lower/middle/upper lobe of left/right lung
+* lateral/medial surface of left/right lung
+* lateral/medial surface of lower/middle/upper lobe of left/right lung
+* left/right lung surface
+* lower/middle/upper lobe of left/right lung surface
+* oblique fissure of left/right lung
+* oblique fissure of lower/middle/upper lobe of left/right lung
+* right lung accessory lobe surface
+* base of right lung accessory lobe surface
+* left/right/dorsal/ventral surface of right lung accessory lobe
+
+Annotated 1-dimensional line regions are defined for matching annotated contours digitized from medical images including
+(using only one of the items separated by slash /):
+
+* anterior border of left/right lung
+
+Several fiducial marker points are defined on the lung scaffold, of which the followings are potentially usable when
+digitizing:
+
+* apex of left/right lung
+* laterodorsal tip of middle lobe of right lung
+* medial/ventral base of left/right lung
+* dorsal/ventral apex of right lung accessory lobe
+* left/right dorsal/ventral base of right lung accessory lobe
+
+**Digitization tips to assist fitting:**
+
+1. A proper lung model requires accurate location of all surfaces of the lobes. This requires digitizing the fissures
+and exterior surfaces of the lungs. It's not a requirement to use the most specific annotation group for a surface
+(e.g. for a particular lobe instead of the whole left/right lung) but it may make fitting more efficient.
+
+2. Digitize any fiducial markers you can identify as these are gold standard locations which can be highly weighted in
+the fit.
+
+3. The sharp anterior edges of the left/right lungs are difficult to fit. At a minimum it's important to have data
+marker points for the ``ventral base of left/right lung`` fiducial markers. To properly fit the rest of the edges it's
+best to have data points/contours along these edges annotated with the 1-D ``anterior border of left/right lung``
+terms, which can be weighted highly in the fit. Annotating with lateral/medial surfaces may work in some cases, but just
+using lung/lobe surface groups can be problematic; in both cases there may not be enough data to definitively pull the
+edge into position during the fit.
+
+4. For fitting a lung scaffold with open fissures, annotate fissure data points/contours which clearly belong to only
+one lobe with the fissure term specific to that lobe. Where fissures are too close to distinguish the lobe they are on,
+annotate digitized data the generic fissure term (not for a specific lobe). Doing this allows the same data to fit all
+lobes correctly.
diff --git a/docs/scaffolds/lung4.rst b/docs/scaffolds/lung4.rst
new file mode 100644
index 00000000..f9a5b11d
--- /dev/null
+++ b/docs/scaffolds/lung4.rst
@@ -0,0 +1,98 @@
+Lung 4 Scaffold
+===============
+
+The `3D Lung 4` Scaffold built from ``class MeshType_3d_lung4`` builds a volumetric model of the left and right human lungs with open fissures between lobes, but with a line of common nodes joining lobes of each lung on the hilum axis.
+
+It is built in an idealised form on an ellipsoid and shaped slightly to give it a sharp anterior edge, concave base and slight curvature around the heart.
+The scaffold mesh consists of all hexahedral (cube) elements without any collapsed faces/edges, with geometry interpolated by smooth tricubic Hermite Serendipity basis functions. The resolution / element density of the scaffold's mesh is fully controllable without changing its shape. It is suitable for fitting to a wide range of human lung shapes.
+
+.. _fig-scaffoldmaker-lung4-human:
+
+.. figure:: ../_images/lung4_human1_coarse.jpg
+ :align: center
+
+ 3D Lung 4 Scaffold with `Human 1 Coarse` parameter set. The oblique fissures of the left and right lungs and the horizontal fissure of the right lung are shaded.
+
+.. _fig-scaffoldmaker-lung4-fitted-left-right:
+
+.. figure:: ../_images/lung4_left_right_fitted.jpg
+ :align: center
+
+ Lung 4 Scaffold fitted to left and right human lung data. Lower lobes are blue, right middle lobe is pink and upper lobes are green.
+
+The lung scaffold is only a representation of the volumetric spaces of the lobes and does not include representations of the pulmonary airways, blood vessels or alveoli.
+
+Variants
+--------
+
+This scaffold currently only provides parameters for human lungs, but coarse, medium and fine mesh resolutions are provided in the idealised human shape, plus an unmodified ellipsoid.
+
+Scaffold parameters are presented to allow other numbers of elements and shapes, but for a given population study it's expected that all fitted subjects start from the same scaffold. Lungs have quite a lot of variability in shape between subjects, so in the early stages of fitting it may be preferred to perform a gross shape change, resolve most non-linearities over 3 or more fit steps/iterations, and select `Update reference state` on the last of these iterations.
+
+Coordinates
+-----------
+
+The lung scaffold defines only a geometric ``coordinates`` field which gives an approximate, idealized unit-scale representation of the lung shapes, which is intended to be fitted to actual data for a subject.
+
+It is planned in future to use the coordinates generated for the ``Ellipsoid`` parameter set as common material coordinates (i.e. ``lung coordinates``) across all subjects, to use for embedding data within the lobes.
+
+The lung scaffold supports limited refinement/resampling to trilinear elements by checking *Refine* (set parameter to ``True``) with chosen
+*Refine number of elements* parameter. Be aware that only the ``coordinates`` field is currently defined on the refined
+mesh (but annotations are transferred).
+
+Annotations
+-----------
+
+Important anatomical regions of the lungs are defined by groups of elements (or faces, edges and nodes/points) and
+annotated with standard term names and identifiers from a controlled vocabulary.
+
+Annotated 3-dimensional volume regions are defined by groups of 3-D elements including (using only one of the items
+separated by slash /, and noting the middle lobe is only on the right so the upper lobe of the right lung never touches the base):
+
+* lung (everything)
+* left/right lung
+* lower/middle/upper lobe of left/right lung
+
+.. note::
+
+ Terms for volume regions such as the above are not to be used for digitized contours. They are used for applying different material properties in models and the strain/curvature penalty (stiffness) parameters in fitting.
+
+
+Annotated 2-dimensional surface regions are defined for matching annotated contours digitized from medical images
+including:
+
+* base of lower/middle/upper lobe of left/right lung surface
+* horizontal fissure of middle/upper lobe of right lung
+* lateral/medial surface of left/right lung
+* lateral/medial surface of lower/middle/upper lobe of left/right lung
+* lower/middle/upper lobe of left/right lung surface
+* oblique fissure of lower/middle/upper lobe of left/right lung
+
+.. note::
+
+ 1. Horizontal and oblique fissure surfaces (and edges below) are labelled specifically for the lobe touching them. For fitting, data points on the fissures must be duplicated for each side to keep the respective surfaces together on the fissure data.
+ 2. The surfaces of the oblique fissure of the left-upper lobe below the hilum, and the right-middle lobe are the same as the base of the left-upper and right-middle lobes, respectively. This is to allow material to slide in or out of the fissure as there are huge variations of relative fissure and base sizes across lung populations. This is seen in :numref:`fig-scaffoldmaker-lung4-fitted-left-right` where the base of the middle lobe is very large on the right.
+
+Annotated 1-dimensional line regions are defined for matching annotated contours digitized from medical images including:
+
+* anterior edge of middle lobe of right lung
+* antero-posterior edge of upper lobe of left/right lung
+* base edge of oblique fissure of lower lobe of left/right lung
+* lateral/medial edge of most surface groups listed earlier
+* posterior edge of lower lobe of left/right lung
+
+.. note::
+
+ 1. For fitting, it is a good idea to give edge data higher weights e.g. 10x default weights to pull material into some of the sharp extremeties.
+ 2. For fitting, edges of the fissures may need to be weighted higher to help constrain them to stay together, particularly the medial edges of the oblique fissure of the lower/upper lobe of the left lung, as the medial surface of the left lung can be quite distorted.
+
+Several fiducial marker points are defined on the lung scaffold, as shown in :numref:`fig-scaffoldmaker-lung4-human`:
+
+* apex of left/right lung
+* dorsal/medial/ventral base of left/right lung
+* laterodorsal tip of middle lobe of right lung
+
+.. note::
+
+ 1. It is not recommended that these fiducial markers be used in the fitting process as their positions are somewhat amorphous and without them the fit will smoothly spread out elements between fissures and prescribed edges. Fitting with markers and high weights may cause severe distortions.
+ 2. These points are very useful for subsequent processing as they can help align fitted scaffolds for PCA and other tasks.
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_ellipsoid1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_ellipsoid1.py
index 78d30cb0..fe407301 100644
--- a/src/scaffoldmaker/meshtypes/meshtype_3d_ellipsoid1.py
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_ellipsoid1.py
@@ -129,8 +129,7 @@ def generateBaseMesh(cls, region, options):
fieldmodule = region.getFieldmodule()
coordinates = find_or_create_field_coordinates(fieldmodule)
- ellipsoid = EllipsoidMesh(a, b, c, element_counts, transition_element_count,
- axis2_x_rotation_radians, axis3_x_rotation_radians, surface_only)
+ ellipsoid = EllipsoidMesh(a, b, c, element_counts, transition_element_count, surface_only)
left_group = AnnotationGroup(region, ("left", ""))
right_group = AnnotationGroup(region, ("right", ""))
@@ -157,7 +156,7 @@ def generateBaseMesh(cls, region, options):
ellipsoid.set_nway_derivative_factor(nway_derivative_factor)
ellipsoid.set_surface_d3_mode(surface_d3_mode)
- ellipsoid.build()
+ ellipsoid.build(axis2_x_rotation_radians, axis3_x_rotation_radians)
ellipsoid.generate_mesh(fieldmodule, coordinates)
return annotation_groups, None
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py
index 1227765b..255b8235 100644
--- a/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py
@@ -505,10 +505,6 @@ def createEsophagusMesh3d(region, options, networkLayout, nextNodeIdentifier, ne
markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location")
nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
- markerPoints = markerGroup.getOrCreateNodesetGroup(nodes)
- markerTemplateInternal = nodes.createNodetemplate()
- markerTemplateInternal.defineField(markerName)
- markerTemplateInternal.defineField(markerLocation)
markerNames = ["proximodorsal midpoint on serosa of upper esophageal sphincter",
"proximoventral midpoint on serosa of upper esophageal sphincter",
@@ -532,19 +528,15 @@ def createEsophagusMesh3d(region, options, networkLayout, nextNodeIdentifier, ne
[0.0, 1.0, 1.0],
[xi1Pi, 1.0, 1.0]]
+ esophagusNodesetGroup = esophagusGroup.getNodesetGroup(nodes)
for n in range(len(markerNames)):
markerGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region,
- get_esophagus_term(markerNames[n]))
+ get_esophagus_term(markerNames[n]), isMarker=True)
markerElement = mesh.findElementByIdentifier(markerElementIdentifiers[n])
markerXi = markerXis[n]
- cache.setMeshLocation(markerElement, markerXi)
- markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal)
+ markerNode = markerGroup.createMarkerNode(nodeIdentifier, element=markerElement, xi=markerXi)
+ esophagusNodesetGroup.addNode(markerNode)
nodeIdentifier += 1
- cache.setNode(markerPoint)
- markerName.assignString(cache, markerGroup.getName())
- markerLocation.assignMeshLocation(cache, markerElement, markerXi)
- for group in [esophagusGroup, markerGroup]:
- group.getNodesetGroup(nodes).addNode(markerPoint)
fm.endChange()
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py
index 97bbdd55..d00be9c6 100644
--- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py
@@ -129,7 +129,7 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid
# AnnotationGroup(region, get_heart_term("posterior cusp of aortic valve")),
# AnnotationGroup(region, get_heart_term("right cusp of aortic valve")),
# AnnotationGroup(region, get_heart_term("left cusp of aortic valve")) ]
- cuspGroup = AnnotationGroup(region, get_heart_term("aortic valve leaflet")),
+ cuspGroup = AnnotationGroup(region, get_heart_term("aortic valve leaflet"))
cuspGroups = [cuspGroup] * 3
else:
arterialRootGroup = AnnotationGroup(region, get_heart_term("root of pulmonary trunk"))
@@ -138,7 +138,7 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid
# AnnotationGroup(region, get_heart_term("right cusp of pulmonary valve")),
# AnnotationGroup(region, get_heart_term("anterior cusp of pulmonary valve")),
# AnnotationGroup(region, get_heart_term("left cusp of pulmonary valve")) ]
- cuspGroup = AnnotationGroup(region, get_heart_term("pulmonary valve leaflet")),
+ cuspGroup = AnnotationGroup(region, get_heart_term("pulmonary valve leaflet"))
cuspGroups = [cuspGroup] * 3
allGroups = [ arterialRootGroup ] # groups that all elements in scaffold will go in
@@ -528,11 +528,9 @@ def refineMesh(cls, meshrefinement, options):
numberInXi2 = refineElementsCountSurface
numberInXi3 = refineElementsCountThroughWall
for cusp in range(3):
- lastShareNodeIds = lastShareNodeCoordinates = None
for e in range(2):
element = meshrefinement._sourceElementiterator.next()
- lastShareNodeIds, lastShareNodeCoordinates = meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3,
- addNewNodesToOctree=False, shareNodeIds=lastShareNodeIds, shareNodeCoordinates=lastShareNodeCoordinates)
+ meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3)
def getSemilunarValveSinusPoints(centre, axisSide1, axisSide2, radius, sinusRadialDisplacement, startMidCusp, elementsCountAround = 6):
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung3.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung3.py
index c4107a1e..9b9d1f7c 100644
--- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung3.py
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung3.py
@@ -1,5 +1,5 @@
"""
-Generates a lung scaffold by deforming a hemisphere.
+Generates a lung scaffold by deforming an ellipsoid.
"""
from cmlibs.utils.zinc.field import find_or_create_field_coordinates
from cmlibs.zinc.field import Field
@@ -17,7 +17,7 @@
class MeshType_3d_lung3(Scaffold_base):
"""
- Generates a lung scaffold by deforming a hemisphere.
+ Generates a lung scaffold by deforming an ellipsoid.
"""
@classmethod
@@ -263,7 +263,7 @@ def generateBaseMesh(cls, region, options):
axis3_x_rotation_radians = math.radians(90) - oblique_slope_radians
ellipsoid = EllipsoidMesh(halfDepth, halfBreadth, halfHeight, elementCounts, elementsCountTransition,
- axis2_x_rotation_radians, axis3_x_rotation_radians, surface_only)
+ surface_only)
if lung == leftLung:
octant_group_lists = []
@@ -292,7 +292,7 @@ def generateBaseMesh(cls, region, options):
ellipsoid.set_octant_group_lists(octant_group_lists)
- ellipsoid.build()
+ ellipsoid.build(axis2_x_rotation_radians, axis3_x_rotation_radians)
nodeIdentifier, elementIdentifier = ellipsoid.generate_mesh(fieldmodule, coordinates, nodeIdentifier, elementIdentifier)
for lung in lungs:
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung4.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung4.py
new file mode 100644
index 00000000..f316fe34
--- /dev/null
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung4.py
@@ -0,0 +1,904 @@
+"""
+Generates a lung scaffold with common hilum but open fissures with lobes built from 1/6 ellipsoid segments.
+"""
+from cmlibs.maths.vectorops import magnitude
+from cmlibs.utils.zinc.field import find_or_create_field_coordinates
+from cmlibs.zinc.element import Element
+from cmlibs.zinc.field import Field
+
+from scaffoldmaker.annotation.annotationgroup import (
+ AnnotationGroup, findAnnotationGroupByName, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm)
+from scaffoldmaker.annotation.lung_terms import get_lung_term
+from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
+from scaffoldmaker.utils.geometry import getEllipsePointAtTrueAngle
+from scaffoldmaker.utils.interpolation import getNearestLocationOnCurve
+from scaffoldmaker.utils.meshrefinement import MeshRefinement
+from scaffoldmaker.utils.ellipsoidmesh import EllipsoidMesh, EllipsoidSurfaceD3Mode
+from scaffoldmaker.utils.zinc_utils import get_mesh_first_element_with_node, translate_nodeset_coordinates
+
+import copy
+import math
+
+
+class MeshType_3d_lung4(Scaffold_base):
+ """
+ Generates a lung scaffold with common hilum but open fissures with lobes built from 1/6 ellipsoid segments.
+ """
+
+ @classmethod
+ def getName(cls):
+ return "3D Lung 4"
+
+ @classmethod
+ def getParameterSetNames(cls):
+ return [
+ "Default",
+ "Human 1 Coarse",
+ "Human 1 Medium",
+ "Human 1 Fine",
+ "Ellipsoid"
+ ]
+
+ @classmethod
+ def getDefaultOptions(cls, parameterSetName="Default"):
+ options = {}
+ useParameterSetName = "Human 1 Coarse" if (parameterSetName == "Default") else parameterSetName
+ options["Left lung"] = True
+ options["Right lung"] = True
+ options["Ellipsoid height"] = 1.0
+ options["Ellipsoid dorsal-ventral size"] = 0.8
+ options["Ellipsoid medial-lateral size"] = 0.4
+ options["Left-right lung spacing"] = 0.5
+ options["Lower lobe extension"] = 0.3
+ options["Refine"] = False
+ options["Refine number of elements"] = 4
+
+ if "Medium" in useParameterSetName:
+ options["Number of elements lateral"] = 4
+ options["Number of elements lower lobe extension"] = 3
+ options["Number of elements oblique"] = 8
+ options["Number of transition elements"] = 1
+ elif "Fine" in useParameterSetName:
+ options["Number of elements lateral"] = 6
+ options["Number of elements lower lobe extension"] = 4
+ options["Number of elements oblique"] = 12
+ options["Number of transition elements"] = 1
+ else:
+ options["Number of elements lateral"] = 4
+ options["Number of elements lower lobe extension"] = 2
+ options["Number of elements oblique"] = 6
+ options["Number of transition elements"] = 1
+
+ if "Human" in useParameterSetName:
+ options["Lower lobe base concavity"] = 0.1
+ options["Ventral edge sharpness factor"] = 0.9
+ options["Cardiac curvature"] = 2.0
+ else:
+ options["Lower lobe base concavity"] = 0.0
+ options["Ventral edge sharpness factor"] = 0.0
+ options["Cardiac curvature"] = 0.0
+
+ return options
+
+ @classmethod
+ def getOrderedOptionNames(cls):
+ return [
+ "Left lung",
+ "Right lung",
+ "Number of elements lateral",
+ "Number of elements lower lobe extension",
+ "Number of elements oblique",
+ "Number of transition elements",
+ "Ellipsoid height",
+ "Ellipsoid dorsal-ventral size",
+ "Ellipsoid medial-lateral size",
+ "Left-right lung spacing",
+ "Lower lobe base concavity",
+ "Lower lobe extension",
+ "Ventral edge sharpness factor",
+ "Cardiac curvature",
+ "Refine",
+ "Refine number of elements"
+ ]
+
+ @classmethod
+ def checkOptions(cls, options):
+ dependent_changes = False
+
+ max_transition_count = None
+ for key in [
+ "Number of elements lateral",
+ "Number of elements oblique"
+ ]:
+ if options[key] < 4:
+ options[key] = 4
+ elif options[key] % 2:
+ options[key] += 1
+ transition_count = (options[key] // 2) - 1
+ if (max_transition_count is None) or (transition_count < max_transition_count):
+ max_transition_count = transition_count
+
+ if options["Number of transition elements"] < 1:
+ options["Number of transition elements"] = 1
+ elif options["Number of transition elements"] > max_transition_count:
+ options["Number of transition elements"] = max_transition_count
+ dependent_changes = True
+
+ for key, default in [
+ ("Ellipsoid height", 1.0),
+ ("Ellipsoid dorsal-ventral size", 0.8),
+ ("Ellipsoid medial-lateral size", 0.4)
+ ]:
+ if options[key] <= 0.0:
+ options[key] = default
+ for key in [
+ "Lower lobe base concavity",
+ "Lower lobe extension",
+ "Cardiac curvature"
+ ]:
+ if options[key] < 0.0:
+ options[key] = 0.0
+ if options["Lower lobe extension"] == 0.0:
+ options["Number of elements lower lobe extension"] = 0
+ dependent_changes = True
+ elif options["Number of elements lower lobe extension"] < 1:
+ options["Number of elements lower lobe extension"] = 1
+ dependent_changes = True
+ depth = options["Ellipsoid dorsal-ventral size"]
+ height = options["Ellipsoid height"]
+ max_extension = 0.99 * magnitude(getEllipsePointAtTrueAngle(depth / 2.0, height / 2.0, math.pi / 3.0))
+ if options["Lower lobe extension"] > max_extension:
+ options["Lower lobe extension"] = max_extension
+
+ if options["Left-right lung spacing"] < 0.0:
+ options["Left-right lung spacing"] = 0.0
+
+ for dimension in [
+ "Ventral edge sharpness factor"
+ ]:
+ if options[dimension] < 0.0:
+ options[dimension] = 0.0
+ elif options[dimension] > 1.0:
+ options[dimension] = 1.0
+
+ if options['Refine number of elements'] < 1:
+ options['Refine number of elements'] = 1
+
+ return dependent_changes
+
+ @classmethod
+ def generateBaseMesh(cls, region, options):
+ """
+ Generate the base tricubic Hermite mesh. See also generateMesh().
+ :param region: Zinc region to define model in. Must be empty.
+ :param options: Dict containing options. See getDefaultOptions().
+ :return: list of AnnotationGroup, None
+ """
+ has_left_lung = options["Left lung"]
+ has_right_lung = options["Right lung"]
+
+ elements_count_lateral = options["Number of elements lateral"]
+ elements_count_lower_extension = options["Number of elements lower lobe extension"]
+ elements_count_oblique = options["Number of elements oblique"]
+ elements_count_transition = options["Number of transition elements"]
+ lung_spacing = options["Left-right lung spacing"] * 0.5
+ lower_lobe_extension = options["Lower lobe extension"]
+ lower_lobe_base_concavity = options["Lower lobe base concavity"]
+ ventral_sharpness_factor = options["Ventral edge sharpness factor"]
+ cardiac_curvature = options["Cardiac curvature"]
+ ellipsoid_ml_size = options["Ellipsoid medial-lateral size"]
+ ellipsoid_dv_size = options["Ellipsoid dorsal-ventral size"]
+ ellipsoid_height = options["Ellipsoid height"]
+
+ fieldmodule = region.getFieldmodule()
+ fieldcache = fieldmodule.createFieldcache()
+ coordinates = find_or_create_field_coordinates(fieldmodule)
+ nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+ mesh = fieldmodule.findMeshByDimension(3)
+
+ # annotation groups & nodeset groups
+ lung_group = AnnotationGroup(region, get_lung_term("lung"))
+
+ if has_left_lung:
+ left_lung_group = AnnotationGroup(region, get_lung_term("left lung"))
+ lower_left_lung_group = AnnotationGroup(region, get_lung_term("lower lobe of left lung"))
+ upper_left_lung_group = AnnotationGroup(region, get_lung_term("upper lobe of left lung"))
+ left_anterior_lung_group = AnnotationGroup(region, ("anterior left lung", ""))
+ left_lateral_lung_group = AnnotationGroup(region, ("lateral left lung", ""))
+ left_medial_lung_group = AnnotationGroup(region, ("medial left lung", ""))
+ left_annotation_groups = [left_lung_group, lower_left_lung_group, upper_left_lung_group,
+ left_anterior_lung_group, left_lateral_lung_group, left_medial_lung_group]
+ else:
+ left_lung_group = lower_left_lung_group = upper_left_lung_group = None
+ left_lateral_lung_group = left_medial_lung_group = left_anterior_lung_group = None
+ left_annotation_groups = []
+
+ if has_right_lung:
+ right_lung_group = AnnotationGroup(region, get_lung_term("right lung"))
+ lower_right_lung_group = AnnotationGroup(region, get_lung_term("lower lobe of right lung"))
+ middle_right_lung_group = AnnotationGroup(region, get_lung_term("middle lobe of right lung"))
+ upper_right_lung_group = AnnotationGroup(region, get_lung_term("upper lobe of right lung"))
+ right_anterior_lung_group = AnnotationGroup(region, ("anterior right lung", ""))
+ right_lateral_lung_group = AnnotationGroup(region, ("lateral right lung", ""))
+ right_medial_lung_group = AnnotationGroup(region, ("medial right lung", ""))
+ right_annotation_groups = [
+ right_lung_group, lower_right_lung_group, middle_right_lung_group, upper_right_lung_group,
+ right_anterior_lung_group, right_lateral_lung_group, right_medial_lung_group]
+ else:
+ right_lung_group = lower_right_lung_group = middle_right_lung_group = upper_right_lung_group = None
+ right_anterior_lung_group = right_lateral_lung_group = right_medial_lung_group = None
+ right_annotation_groups = []
+
+ box_group = AnnotationGroup(region, ("box", ""))
+ transition_group = AnnotationGroup(region, ("transition", ""))
+
+ annotation_groups = [lung_group] + left_annotation_groups + right_annotation_groups + \
+ [box_group, transition_group]
+
+ half_ml_size = ellipsoid_ml_size * 0.5
+ half_dv_size = ellipsoid_dv_size * 0.5
+ half_height = ellipsoid_height * 0.5
+
+ pi__3 = math.pi / 3.0
+ sin_pi__3 = math.sin(pi__3)
+
+ left_lung, right_lung = 0, 1
+ lungs = [lung for show, lung in [(has_left_lung, left_lung), (has_right_lung, right_lung)] if show]
+ node_identifier, element_identifier = 1, 1
+
+ marker_name_element_xi = []
+
+ # need to build both left and right lungs then delete any not used so element, face, line, node numbers
+ # are the same if only left or right
+ for lung in [left_lung, right_lung]:
+
+ if lung == left_lung:
+ if has_left_lung:
+ lower_octant_group_lists = []
+ upper_octant_group_lists = []
+ for octant in range(8):
+ octant_group_list = [group.getGroup() for group in
+ [lung_group, left_lung_group, lower_left_lung_group] +
+ [left_medial_lung_group if (octant & 1) else left_lateral_lung_group]]
+ lower_octant_group_lists.append(octant_group_list)
+ octant_group_list = [group.getGroup() for group in
+ [lung_group, left_lung_group, upper_left_lung_group] +
+ [left_medial_lung_group if (octant & 1) else left_lateral_lung_group]]
+ if octant & 2:
+ octant_group_list.append(left_anterior_lung_group.getGroup())
+ upper_octant_group_lists.append(octant_group_list)
+ else:
+ lower_octant_group_lists = upper_octant_group_lists = None
+ middle_octant_group_lists = None
+ else:
+ if has_right_lung:
+ lower_octant_group_lists = []
+ middle_octant_group_lists = []
+ upper_octant_group_lists = []
+ for octant in range(8):
+ octant_group_list = [group.getGroup() for group in
+ [lung_group, right_lung_group, lower_right_lung_group] +
+ [right_lateral_lung_group if (octant & 1) else right_medial_lung_group]]
+ lower_octant_group_lists.append(octant_group_list)
+ octant_group_list = [group.getGroup() for group in
+ [lung_group, right_lung_group, middle_right_lung_group] +
+ [right_lateral_lung_group if (octant & 1) else right_medial_lung_group]]
+ if octant & 2:
+ octant_group_list.append(right_anterior_lung_group.getGroup())
+ middle_octant_group_lists.append(octant_group_list)
+ octant_group_list = [group.getGroup() for group in
+ [lung_group, right_lung_group, upper_right_lung_group] +
+ [right_lateral_lung_group if (octant & 1) else right_medial_lung_group]]
+ if octant & 2:
+ octant_group_list.append(right_anterior_lung_group.getGroup())
+ upper_octant_group_lists.append(octant_group_list)
+ else:
+ lower_octant_group_lists = middle_octant_group_lists = upper_octant_group_lists = None
+
+ element_counts = [elements_count_lateral, elements_count_oblique, elements_count_oblique]
+ lower_ellipsoid = EllipsoidMesh(
+ half_ml_size, half_dv_size, half_height, element_counts, elements_count_transition)
+ lower_ellipsoid.set_surface_d3_mode(EllipsoidSurfaceD3Mode.SURFACE_NORMAL_PLANE_PROJECTION)
+ lower_ellipsoid.set_box_transition_groups(box_group.getGroup(), transition_group.getGroup())
+ upper_ellipsoid = EllipsoidMesh(
+ half_ml_size, half_dv_size, half_height, element_counts, elements_count_transition)
+ upper_ellipsoid.set_surface_d3_mode(EllipsoidSurfaceD3Mode.SURFACE_NORMAL_PLANE_PROJECTION)
+ upper_ellipsoid.set_box_transition_groups(box_group.getGroup(), transition_group.getGroup())
+ if lung == right_lung:
+ middle_ellipsoid = EllipsoidMesh(
+ half_ml_size, half_dv_size, half_height, element_counts, elements_count_transition)
+ middle_ellipsoid.set_surface_d3_mode(EllipsoidSurfaceD3Mode.SURFACE_NORMAL_PLANE_PROJECTION)
+ middle_ellipsoid.set_box_transition_groups(box_group.getGroup(), transition_group.getGroup())
+ else:
+ middle_ellipsoid = upper_ellipsoid
+ half_counts = [count // 2 for count in element_counts]
+ octant1 = middle_ellipsoid.build_octant(half_counts, -pi__3, 0.0)
+ middle_ellipsoid.merge_octant(octant1, quadrant=3)
+ if lung == right_lung:
+ middle_ellipsoid.copy_to_negative_axis1()
+
+ # save hilum coordinates for all other lobes
+ hilum_x = []
+ ox = octant1.get_parameters()
+ box_count1 = octant1.get_box_counts()[0]
+ for n1 in range(element_counts[0] + 1):
+ mirror_x = n1 < half_counts[0]
+ o1 = abs(n1 - half_counts[0])
+ parameters = ox[0][0][o1]
+ obox = o1 <= box_count1
+ parameters = [
+ copy.copy(parameters[0]),
+ copy.copy(parameters[3 if obox else 2]),
+ [-d for d in parameters[2 if obox else 1]],
+ copy.copy(parameters[1 if obox else 3])
+ ]
+ if mirror_x:
+ for i in range(3):
+ parameters[i][0] = -parameters[i][0]
+ hilum_x.append(parameters)
+
+ octant2 = upper_ellipsoid.build_octant(half_counts, 0.0, pi__3)
+ upper_ellipsoid.merge_octant(octant2, quadrant=0)
+ octant3 = upper_ellipsoid.build_octant(half_counts, pi__3, 2.0 * pi__3)
+ upper_ellipsoid.merge_octant(octant3, quadrant=1)
+ upper_ellipsoid.copy_to_negative_axis1()
+
+ octant4 = lower_ellipsoid.build_octant(half_counts, 2.0 * pi__3, math.pi,
+ lower_lobe_extension, elements_count_lower_extension)
+ # merge into separate lower ellipsoid to have space for extension elements
+ lower_ellipsoid_mesh = EllipsoidMesh(
+ half_ml_size, half_dv_size, half_height,
+ [element_counts[0], element_counts[1], element_counts[2] + 2 * elements_count_lower_extension],
+ elements_count_transition)
+ lower_ellipsoid_mesh.set_surface_d3_mode(EllipsoidSurfaceD3Mode.SURFACE_NORMAL_PLANE_PROJECTION)
+ lower_ellipsoid_mesh.set_box_transition_groups(box_group.getGroup(), transition_group.getGroup())
+ lower_ellipsoid_mesh.merge_octant(octant4, quadrant=1)
+ lower_ellipsoid_mesh.copy_to_negative_axis1()
+
+ node_layout_manager = lower_ellipsoid.get_node_layout_manager()
+ node_layout_permuted = node_layout_manager.getNodeLayoutRegularPermuted(d3Defined=True)
+ for n1 in range(element_counts[0] + 1):
+ lower_ellipsoid_mesh.set_node_parameters(
+ n1, half_counts[1], element_counts[2] + 2 * elements_count_lower_extension - half_counts[2],
+ hilum_x[n1], node_layout=node_layout_permuted)
+ lower_ellipsoid_mesh.set_octant_group_lists(lower_octant_group_lists)
+ node_identifier, element_identifier = lower_ellipsoid_mesh.generate_mesh(
+ fieldmodule, coordinates, node_identifier, element_identifier)
+
+ for ellipsoid in [middle_ellipsoid, upper_ellipsoid] if (lung == right_lung) else [upper_ellipsoid]:
+ node_layout_manager = ellipsoid.get_node_layout_manager()
+ node_layout_6way = node_layout_manager.getNodeLayout6Way12(d3Defined=True)
+ for n1 in range(element_counts[0] + 1):
+ nid = (lower_ellipsoid_mesh.get_node_identifier(
+ n1, half_counts[1], element_counts[2] + 2 * elements_count_lower_extension - half_counts[2])
+ if (((lung == left_lung) and (n1 >= half_counts[0])) or
+ ((lung == right_lung) and (n1 <= half_counts[0]))) else None)
+ ellipsoid.set_node_parameters(n1, half_counts[1], half_counts[2], hilum_x[n1],
+ nid, node_layout=node_layout_6way)
+ ellipsoid.set_octant_group_lists(
+ middle_octant_group_lists if ((ellipsoid == middle_ellipsoid) and
+ (ellipsoid != upper_ellipsoid)) else upper_octant_group_lists)
+ node_identifier, element_identifier = ellipsoid.generate_mesh(
+ fieldmodule, coordinates, node_identifier, element_identifier)
+
+ # find elements for marker points
+ if ((lung == left_lung) and has_left_lung) or ((lung == right_lung) and has_right_lung):
+ nid = upper_ellipsoid.get_node_identifier(
+ elements_count_lateral // 2, 0, element_counts[2])
+ group = upper_left_lung_group if (lung == left_lung) else upper_right_lung_group
+ if nid and group:
+ marker_name_element_xi.append((
+ "apex of left lung" if (lung == left_lung) else "apex of right lung",
+ get_mesh_first_element_with_node(
+ group.getMeshGroup(mesh), coordinates, nodes.findNodeByIdentifier(nid)),
+ [1.0, 1.0, 1.0]))
+ nid = lower_ellipsoid_mesh.get_node_identifier(
+ elements_count_lateral // 2, 0, elements_count_oblique // 2 + elements_count_lower_extension)
+ group = lower_left_lung_group if (lung == left_lung) else lower_right_lung_group
+ if nid and group:
+ marker_name_element_xi.append((
+ "dorsal base of left lung" if (lung == left_lung) else "dorsal base of right lung",
+ get_mesh_first_element_with_node(
+ group.getMeshGroup(mesh), coordinates, nodes.findNodeByIdentifier(nid)),
+ [1.0, 0.0, 1.0]))
+ if lung == right_lung:
+ nid = middle_ellipsoid.get_node_identifier(
+ elements_count_lateral, elements_count_oblique // 2, elements_count_oblique // 2)
+ group = middle_right_lung_group
+ if nid and group:
+ marker_name_element_xi.append((
+ "laterodorsal tip of middle lobe of right lung",
+ get_mesh_first_element_with_node(
+ group.getMeshGroup(mesh), coordinates, nodes.findNodeByIdentifier(nid)),
+ [0.0, 1.0, 1.0]))
+ # need to find element xi where medial base edge crosses 0.0
+ group = lower_left_lung_group if (lung == left_lung) else lower_right_lung_group
+ if group:
+ n1 = elements_count_lateral if (lung == left_lung) else 0
+ n3 = elements_count_oblique // 2 + elements_count_lower_extension
+ last_n2 = None
+ last_nid = None
+ last_nx = None
+ last_nd1 = None
+ for n2 in range(0, elements_count_oblique // 2 + 1):
+ nid = lower_ellipsoid_mesh.get_node_identifier(n1, n2, n3)
+ if nid:
+ nx, nd1 = lower_ellipsoid_mesh.get_node_parameters(n1, n2, n3)[:2]
+ if last_nid:
+ if (last_nx[1] <= 0.0) and (nx[1] >= 0.0):
+ element = get_mesh_first_element_with_node(
+ group.getMeshGroup(mesh), coordinates, nodes.findNodeByIdentifier(
+ nid if ((lung == left_lung) or (last_n2 == 0)) else last_nid))
+ if element.isValid():
+ # xi1 goes counter-clockwise from above so different for left and right
+ if lung == left_lung:
+ curve_location, x = getNearestLocationOnCurve(
+ [last_nx[1:2], nx[1:2]], [last_nd1[1:2], nd1[1:2]], [0.0],
+ startLocation=(0, -last_nx[1] / (nx[1] - last_nx[1])))
+ else:
+ curve_location, x = getNearestLocationOnCurve(
+ [nx[1:2], last_nx[1:2]], [nd1[1:2], last_nd1[1:2]], [0.0],
+ startLocation=(0, nx[1] / (nx[1] - last_nx[1])))
+ marker_name_element_xi.append((
+ "medial base of left lung" if (lung == left_lung) else
+ "medial base of right lung",
+ element, [curve_location[1], 0.0, 1.0]))
+ break
+ last_n2 = n2
+ last_nid = nid
+ last_nx = nx
+ last_nd1 = nd1
+ ventral_ellipsoid = upper_ellipsoid if (lung == left_lung) else middle_ellipsoid
+ nid = ventral_ellipsoid.get_node_identifier(
+ elements_count_lateral // 2, elements_count_oblique // 2, 0)
+ group = upper_left_lung_group if (lung == left_lung) else middle_right_lung_group
+ if nid and group:
+ marker_name_element_xi.append((
+ "ventral base of left lung" if (lung == left_lung) else "ventral base of right lung",
+ get_mesh_first_element_with_node(
+ group.getMeshGroup(mesh), coordinates, nodes.findNodeByIdentifier(nid)),
+ [0.0, 0.0, 1.0]))
+ else:
+ # need to make sure marker node identifiers are consistent if left/right not used
+ for i in range(4 if (lung == left_lung) else 5):
+ marker_name_element_xi.append(("dummy", None, None))
+
+ for lung in lungs:
+ is_left = lung == left_lung
+ lung_nodeset = (left_lung_group if is_left else right_lung_group).getNodesetGroup(nodes)
+
+ if (lower_lobe_base_concavity > 0.0) and (lower_lobe_extension > 0.0):
+ lower_lobe_group = lower_left_lung_group if (lung == left_lung) else lower_right_lung_group
+ form_lower_lobe_base_concavity(
+ lower_lobe_base_concavity, lower_lobe_extension, half_ml_size, half_dv_size, half_height,
+ fieldmodule, coordinates, lower_lobe_group.getNodesetGroup(nodes))
+
+ if ventral_sharpness_factor != 0.0:
+ taper_lung_edge(ventral_sharpness_factor, fieldmodule, coordinates, lung_nodeset, half_dv_size)
+
+ if cardiac_curvature > 0.0:
+ curve_cardiac_anterior(-cardiac_curvature if is_left else cardiac_curvature,
+ half_ml_size, half_dv_size, half_height, fieldmodule, coordinates, lung_nodeset)
+
+ translate_nodeset_coordinates(lung_nodeset, coordinates,
+ [-lung_spacing if is_left else lung_spacing, 0.0, 0.0])
+
+ # marker points; make after regular nodes so higher node numbers
+ lung_nodeset = lung_group.getNodesetGroup(nodes)
+ for marker_name, element, xi in marker_name_element_xi:
+ if element: # skip dummy markers for missing left/right, but increment node_identifier always
+ annotation_group = findOrCreateAnnotationGroupForTerm(
+ annotation_groups, region, get_lung_term(marker_name), isMarker=True)
+ marker_node = annotation_group.createMarkerNode(node_identifier, element=element, xi=xi)
+ if 'left' in marker_name:
+ left_lung_group.getNodesetGroup(nodes).addNode(marker_node)
+ else:
+ right_lung_group.getNodesetGroup(nodes).addNode(marker_node)
+ lung_nodeset.addNode(marker_node)
+ node_identifier += 1
+
+ return annotation_groups, None
+
+ @classmethod
+ def refineMesh(cls, meshRefinement, options):
+ """
+ Refine source mesh into separate region, with change of basis.
+ :param meshRefinement: MeshRefinement, which knows source and target region.
+ :param options: Dict containing options. See getDefaultOptions().
+ """
+ assert isinstance(meshRefinement, MeshRefinement)
+ refine_elements_count = options['Refine number of elements']
+ meshRefinement.refineAllElementsCubeStandard3d(
+ refine_elements_count, refine_elements_count, refine_elements_count)
+
+ @classmethod
+ def defineFaceAnnotations(cls, region, options, annotation_groups):
+ """
+ Add face annotation groups from the highest dimension mesh.
+ Must have defined faces and added subelements for highest dimension groups.
+ :param region: Zinc region containing model.
+ :param options: Dict containing options. See getDefaultOptions().
+ :param annotation_groups: List of annotation groups for top-level elements.
+ New face annotation groups are appended to this list.
+ """
+ fm = region.getFieldmodule()
+ mesh1d = fm.findMeshByDimension(1)
+ mesh2d = fm.findMeshByDimension(2)
+
+ has_left_lung = options["Left lung"]
+ has_right_lung = options["Right lung"]
+
+ if (not has_left_lung) or (not has_right_lung):
+ # destroy elements, faces, lines and nodes now to ensure persistent identifiers when only one side
+ # these are any objects not in the "lung" group:
+ not_lung = fm.createFieldNot(getAnnotationGroupForTerm(annotation_groups, get_lung_term("lung")).getGroup())
+ for mesh in [fm.findMeshByDimension(3), mesh2d, mesh1d]:
+ mesh.destroyElementsConditional(not_lung)
+ nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+ nodes.destroyNodesConditional(not_lung)
+ del nodes
+
+ is_exterior = fm.createFieldIsExterior()
+ is_face_xi1_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI1_0)
+ is_face_xi1_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI1_1)
+ is_face_xi2_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0)
+ is_face_xi2_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_1)
+ is_face_xi3_0 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)
+ is_face_xi3_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)
+ is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, is_face_xi3_1)
+
+ box_group = findAnnotationGroupByName(annotation_groups, "box")
+ is_box = box_group.getGroup()
+ transition_group = findAnnotationGroupByName(annotation_groups, "transition")
+ is_trans = transition_group.getGroup()
+ # this works correctly for faces, but gets extra layers for lines on base and fissures which are exterior:
+ is_on_ellipsoid = fm.createFieldAnd(fm.createFieldAnd(is_exterior_face_xi3_1, is_trans),
+ fm.createFieldNot(is_box))
+
+ is_face_xi1_0_or_xi1_1_or_xi2_0 = fm.createFieldOr(
+ fm.createFieldOr(is_face_xi1_0, is_face_xi1_1), is_face_xi2_0)
+ is_face_xi1_0_or_xi1_1_or_xi2_1 = fm.createFieldOr(
+ fm.createFieldOr(is_face_xi1_0, is_face_xi1_1), is_face_xi2_1)
+ is_face_box30_trans20 = fm.createFieldOr(
+ fm.createFieldAnd(is_box, is_face_xi3_0), fm.createFieldAnd(is_trans, is_face_xi2_0))
+ is_face_box31_trans21 = fm.createFieldOr(
+ fm.createFieldAnd(is_box, is_face_xi3_1), fm.createFieldAnd(is_trans, is_face_xi2_1))
+
+ face_term_conditionals_map = {}
+
+ if has_left_lung:
+ is_lower_left = getAnnotationGroupForTerm(annotation_groups, get_lung_term("lower lobe of left lung")).getGroup()
+ is_upper_left = getAnnotationGroupForTerm(annotation_groups, get_lung_term("upper lobe of left lung")).getGroup()
+ is_anterior_left = findAnnotationGroupByName(annotation_groups, "anterior left lung").getGroup()
+ is_lateral_left = findAnnotationGroupByName(annotation_groups, "lateral left lung").getGroup()
+ is_medial_left = findAnnotationGroupByName(annotation_groups, "medial left lung").getGroup()
+
+ left_face_term_conditionals_map = {
+ "base of lower lobe of left lung surface": (is_lower_left, is_exterior, is_face_box30_trans20),
+ "base of upper lobe of left lung surface": (is_upper_left, is_exterior, is_face_xi1_0_or_xi1_1_or_xi2_0,
+ is_anterior_left),
+ "lateral surface of left lung": (is_lateral_left, is_on_ellipsoid),
+ "lateral surface of lower lobe of left lung": (is_lower_left, is_on_ellipsoid, is_lateral_left),
+ "lateral surface of upper lobe of left lung": (is_upper_left, is_on_ellipsoid, is_lateral_left),
+ "lower lobe of left lung surface": (is_lower_left, is_exterior),
+ "upper lobe of left lung surface": (is_upper_left, is_exterior),
+ "medial surface of left lung": (is_medial_left, is_on_ellipsoid),
+ "medial surface of lower lobe of left lung": (is_lower_left, is_on_ellipsoid, is_medial_left),
+ "medial surface of upper lobe of left lung": (is_upper_left, is_on_ellipsoid, is_medial_left),
+ "oblique fissure of lower lobe of left lung": (is_lower_left, is_exterior,
+ is_face_xi1_0_or_xi1_1_or_xi2_1),
+ "oblique fissure of upper lobe of left lung": (is_upper_left, is_exterior,
+ fm.createFieldOr(is_face_xi1_0_or_xi1_1_or_xi2_0,
+ is_face_box30_trans20)),
+ }
+ face_term_conditionals_map.update(left_face_term_conditionals_map)
+ else:
+ is_lower_left = is_upper_left = is_lateral_left = is_medial_left = None
+
+ if has_right_lung:
+ is_lower_right = getAnnotationGroupForTerm(annotation_groups, get_lung_term("lower lobe of right lung")).getGroup()
+ is_middle_right = getAnnotationGroupForTerm(annotation_groups, get_lung_term("middle lobe of right lung")).getGroup()
+ is_upper_right = getAnnotationGroupForTerm(annotation_groups, get_lung_term("upper lobe of right lung")).getGroup()
+ is_anterior_right = findAnnotationGroupByName(annotation_groups, "anterior right lung").getGroup()
+ is_lateral_right = findAnnotationGroupByName(annotation_groups, "lateral right lung").getGroup()
+ is_medial_right = findAnnotationGroupByName(annotation_groups, "medial right lung").getGroup()
+
+ right_face_term_conditionals_map = {
+ "base of lower lobe of right lung surface": (is_lower_right, is_exterior, is_face_box30_trans20),
+ "base of middle lobe of right lung surface": (is_middle_right, is_exterior, is_face_xi1_0_or_xi1_1_or_xi2_0),
+ "horizontal fissure of middle lobe of right lung": (is_middle_right, is_exterior, is_face_box31_trans21),
+ "horizontal fissure of upper lobe of right lung": (is_upper_right, is_exterior, is_face_box30_trans20, is_anterior_right),
+ "lateral surface of lower lobe of right lung": (is_lower_right, is_on_ellipsoid, is_lateral_right),
+ "lateral surface of middle lobe of right lung": (is_middle_right, is_on_ellipsoid, is_lateral_right),
+ "lateral surface of right lung": (is_lateral_right, is_on_ellipsoid),
+ "lateral surface of upper lobe of right lung": (is_upper_right, is_on_ellipsoid, is_lateral_right),
+ "lower lobe of right lung surface": (is_lower_right, is_exterior),
+ "middle lobe of right lung surface": (is_middle_right, is_exterior),
+ "upper lobe of right lung surface": (is_upper_right, is_exterior),
+ "medial surface of lower lobe of right lung": (is_lower_right, is_on_ellipsoid, is_medial_right),
+ "medial surface of middle lobe of right lung": (is_middle_right, is_on_ellipsoid, is_medial_right),
+ "medial surface of right lung": (is_medial_right, is_on_ellipsoid),
+ "medial surface of upper lobe of right lung": (is_upper_right, is_on_ellipsoid, is_medial_right),
+ "oblique fissure of lower lobe of right lung": (is_lower_right, is_exterior, is_face_xi1_0_or_xi1_1_or_xi2_1),
+ "oblique fissure of middle lobe of right lung": (is_middle_right, is_exterior, is_face_xi1_0_or_xi1_1_or_xi2_0),
+ "oblique fissure of upper lobe of right lung": (is_upper_right, is_exterior, fm.createFieldAnd(is_face_box30_trans20, fm.createFieldNot(is_anterior_right))),
+ }
+ face_term_conditionals_map.update(right_face_term_conditionals_map)
+ else:
+ is_lower_right = is_middle_right = is_upper_right = is_lateral_right = is_medial_right = None
+
+ is_face_conditional = {}
+ for face_term, conditionals in face_term_conditionals_map.items():
+ annotation_group = findOrCreateAnnotationGroupForTerm(annotation_groups, region, get_lung_term(face_term))
+ group = annotation_group.getGroup()
+ conditional = conditionals[0]
+ for add_conditional in conditionals[1:]:
+ conditional = fm.createFieldAnd(conditional, add_conditional)
+ annotation_group.getMeshGroup(mesh2d).addElementsConditional(conditional)
+ annotation_group.addSubelements()
+ is_face_conditional[face_term] = group
+
+ line_term_conditionals_map = {}
+
+ if has_left_lung:
+ is_left_lateral_surface = is_face_conditional["lateral surface of left lung"]
+ is_left_medial_surface = is_face_conditional["medial surface of left lung"]
+ left_line_term_conditionals_map = {
+ "antero-posterior edge of upper lobe of left lung":
+ (is_upper_left, is_left_lateral_surface, is_left_medial_surface),
+ "base edge of oblique fissure of lower lobe of left lung": (
+ is_face_conditional["base of lower lobe of left lung surface"], is_face_conditional["oblique fissure of lower lobe of left lung"]),
+ "lateral edge of base of lower lobe of left lung": (
+ is_face_conditional["base of lower lobe of left lung surface"], is_left_lateral_surface),
+ "lateral edge of base of upper lobe of left lung": (
+ is_face_conditional["base of upper lobe of left lung surface"], is_left_lateral_surface),
+ "lateral edge of oblique fissure of lower lobe of left lung": (
+ is_face_conditional["oblique fissure of lower lobe of left lung"], is_left_lateral_surface),
+ "lateral edge of oblique fissure of upper lobe of left lung": (
+ is_face_conditional["oblique fissure of upper lobe of left lung"], is_left_lateral_surface),
+ "medial edge of base of lower lobe of left lung": (
+ is_face_conditional["base of lower lobe of left lung surface"], is_left_medial_surface),
+ "medial edge of base of upper lobe of left lung": (
+ is_face_conditional["base of upper lobe of left lung surface"], is_left_medial_surface),
+ "medial edge of oblique fissure of lower lobe of left lung": (
+ is_face_conditional["oblique fissure of lower lobe of left lung"], is_left_medial_surface),
+ "medial edge of oblique fissure of upper lobe of left lung": (
+ is_face_conditional["oblique fissure of upper lobe of left lung"], is_left_medial_surface),
+ "posterior edge of lower lobe of left lung":
+ (is_lower_left, is_left_lateral_surface, is_left_medial_surface),
+ }
+ line_term_conditionals_map.update(left_line_term_conditionals_map)
+
+ if has_right_lung:
+ is_right_lateral_surface = is_face_conditional["lateral surface of right lung"]
+ is_right_medial_surface = is_face_conditional["medial surface of right lung"]
+ right_line_term_conditionals_map = {
+ "anterior edge of middle lobe of right lung":
+ (is_middle_right, is_right_lateral_surface, is_right_medial_surface),
+ "antero-posterior edge of upper lobe of right lung":
+ (is_upper_right, is_right_lateral_surface, is_right_medial_surface),
+ "base edge of oblique fissure of lower lobe of right lung": (
+ is_face_conditional["base of lower lobe of right lung surface"], is_face_conditional["oblique fissure of lower lobe of right lung"]),
+ "lateral edge of base of lower lobe of right lung": (
+ is_face_conditional["base of lower lobe of right lung surface"], is_right_lateral_surface),
+ "lateral edge of base of middle lobe of right lung": (
+ is_face_conditional["base of middle lobe of right lung surface"], is_right_lateral_surface),
+ "lateral edge of horizontal fissure of middle lobe of right lung": (
+ is_face_conditional["horizontal fissure of middle lobe of right lung"], is_right_lateral_surface),
+ "lateral edge of horizontal fissure of upper lobe of right lung": (
+ is_face_conditional["horizontal fissure of upper lobe of right lung"], is_right_lateral_surface),
+ "lateral edge of oblique fissure of lower lobe of right lung": (
+ is_face_conditional["oblique fissure of lower lobe of right lung"], is_right_lateral_surface),
+ "lateral edge of oblique fissure of middle lobe of right lung": (
+ is_face_conditional["oblique fissure of middle lobe of right lung"], is_right_lateral_surface),
+ "lateral edge of oblique fissure of upper lobe of right lung": (
+ is_face_conditional["oblique fissure of upper lobe of right lung"], is_right_lateral_surface),
+ "medial edge of base of lower lobe of right lung": (
+ is_face_conditional["base of lower lobe of right lung surface"], is_right_medial_surface),
+ "medial edge of base of middle lobe of right lung": (
+ is_face_conditional["base of middle lobe of right lung surface"], is_right_medial_surface),
+ "medial edge of horizontal fissure of middle lobe of right lung": (
+ is_face_conditional["horizontal fissure of middle lobe of right lung"], is_right_medial_surface),
+ "medial edge of horizontal fissure of upper lobe of right lung": (
+ is_face_conditional["horizontal fissure of upper lobe of right lung"], is_right_medial_surface),
+ "medial edge of oblique fissure of lower lobe of right lung": (
+ is_face_conditional["oblique fissure of lower lobe of right lung"], is_right_medial_surface),
+ "medial edge of oblique fissure of middle lobe of right lung": (
+ is_face_conditional["oblique fissure of middle lobe of right lung"], is_right_medial_surface),
+ "medial edge of oblique fissure of upper lobe of right lung": (
+ is_face_conditional["oblique fissure of upper lobe of right lung"], is_right_medial_surface),
+ "posterior edge of lower lobe of right lung":
+ (is_lower_right, is_right_lateral_surface, is_right_medial_surface)
+ }
+ line_term_conditionals_map.update(right_line_term_conditionals_map)
+
+ for line_term, conditionals in line_term_conditionals_map.items():
+ annotation_group = findOrCreateAnnotationGroupForTerm(annotation_groups, region, (line_term, ""))
+ conditional = conditionals[0]
+ for add_conditional in conditionals[1:]:
+ conditional = fm.createFieldAnd(conditional, add_conditional)
+ annotation_group.getMeshGroup(mesh1d).addElementsConditional(conditional)
+
+ tweak_middle_lobe_tip_edges = True
+ if has_right_lung and tweak_middle_lobe_tip_edges:
+ # tweaks edges of base of middle lobe of right lung to include anterior edge of middle lobe
+ # as these edges of the base data can go sharply up the anterior edge
+ # but this is a difficult and undesirable fit for the middle lobe shape
+ is_anterior_edge = fm.findFieldByName("anterior edge of middle lobe of right lung").castGroup()
+ for base_edge_group_name in [
+ "lateral edge of base of middle lobe of right lung",
+ "medial edge of base of middle lobe of right lung"
+ ]:
+ base_edge_group = fm.findFieldByName(base_edge_group_name).castGroup()
+ base_edge_group.getMeshGroup(mesh1d).addElementsConditional(is_anterior_edge)
+
+ # remove temporary annotation groups
+ for group_name in [
+ "anterior left lung",
+ "anterior right lung",
+ "box",
+ "lateral left lung",
+ "lateral right lung",
+ "medial left lung",
+ "medial right lung",
+ "transition"
+ ]:
+ annotation_group = findAnnotationGroupByName(annotation_groups, group_name)
+ if annotation_group:
+ annotation_groups.remove(annotation_group)
+
+
+def curve_cardiac_anterior(curvature, half_ml_size, half_dv_size, half_height, fieldmodule, coordinates, nodeset):
+ """
+ Transform coordinates by bending with curvature about a centre point the radius in
+ x direction from stationaryPointXY.
+ :param curvature: 1/radius. Must be non-zero.
+ :param half_ml_size: Half medial-lateral ellipsoid size.
+ :param half_dv_size: Half dorsal-ventral ellipsoid size.
+ :param half_height: Half ellipsoid height.
+ :param fieldmodule: Field module being worked with.
+ :param coordinates: The coordinate field, initially circular in y-z plane.
+ :param nodeset: Zinc Nodeset containing nodes to transform.
+ """
+ radius_of_curvature = 1.0 / curvature
+ x_offset = fieldmodule.createFieldConstant(radius_of_curvature)
+
+ pi__3 = math.pi / 3.0
+ cos_pi__3 = math.cos(pi__3)
+ sin_pi__3 = math.sin(pi__3)
+ value_small = fieldmodule.createFieldConstant(1.0E-10)
+ value05 = fieldmodule.createFieldConstant(0.5)
+ x = fieldmodule.createFieldComponent(coordinates, 1)
+ y = fieldmodule.createFieldComponent(coordinates, 2)
+ z = fieldmodule.createFieldComponent(coordinates, 3)
+ yz = fieldmodule.createFieldConcatenate([y, z])
+ r = fieldmodule.createFieldMagnitude(yz)
+ # angle around hilum, not changing
+ theta = fieldmodule.createFieldAtan2(z, y)
+ cos_theta = fieldmodule.createFieldCos(theta)
+ sin_theta = fieldmodule.createFieldSin(theta)
+ # obliqueness factor as want maximum curvature at anterior oblique fissure, none opposite
+ tip_yz = getEllipsePointAtTrueAngle(half_dv_size, half_height, -pi__3)
+ double_max_o = fieldmodule.createFieldConstant(2.0 * magnitude(tip_yz))
+ oblique = fieldmodule.createFieldConstant([cos_pi__3, -sin_pi__3])
+ o = value05 + fieldmodule.createFieldDotProduct(yz, oblique) / double_max_o
+ phi_o = o * o # so curvature mostly at anterior edge
+ # radial curvature
+ alpha = phi_o * r / x_offset
+ cos_alpha = fieldmodule.createFieldCos(alpha)
+ sin_alpha = fieldmodule.createFieldSin(alpha)
+
+ mod_x_offset = x_offset / phi_o
+ mod_offset_x = x + mod_x_offset
+ new_x = mod_offset_x * cos_alpha - mod_x_offset
+
+ new_r = mod_offset_x * sin_alpha
+ new_y = new_r * cos_theta
+ new_z = new_r * sin_theta
+ new_coordinates = fieldmodule.createFieldIf(fieldmodule.createFieldLessThan(o, value_small), coordinates,
+ fieldmodule.createFieldConcatenate([new_x, new_y, new_z]))
+ fieldassignment = coordinates.createFieldassignment(new_coordinates)
+ fieldassignment.setNodeset(nodeset)
+ fieldassignment.assign()
+
+
+def form_lower_lobe_base_concavity(lower_lobe_base_concavity, lower_lobe_extension,
+ half_ml_size, half_dv_size, half_height,
+ fieldmodule, coordinates, nodeset):
+ """
+ Reshape the base of the lower lobe to be concave, tapering off to the oblique fissure.
+ :param lower_lobe_base_concavity: Lower lobe base concavity distance. Note this is reduced by the taper.
+ :param lower_lobe_extension: Distance along oblique fissure below origin that lower lobe extends.
+ :param half_ml_size: Half medial-lateral ellipsoid size.
+ :param half_dv_size: Half dorsal-ventral ellipsoid size.
+ :param half_height: Half ellipsoid height.
+ :param fieldmodule: Fieldmodule owning fields to modify.
+ :param coordinates: Coordinate field to modify.
+ :param nodeset: Nodeset to modify.
+ """
+ pi__3 = math.pi / 3.0
+ cos_pi__3 = math.cos(pi__3)
+ sin_pi__3 = math.sin(pi__3)
+ value0 = fieldmodule.createFieldConstant(0.0)
+ value1 = fieldmodule.createFieldConstant(1.0)
+ value2 = fieldmodule.createFieldConstant(2.0)
+ value3 = fieldmodule.createFieldConstant(3.0)
+ x = fieldmodule.createFieldComponent(coordinates, 1)
+ y = fieldmodule.createFieldComponent(coordinates, 2)
+ z = fieldmodule.createFieldComponent(coordinates, 3)
+ xx = x * x
+ yy = y * y
+ minus_c = fieldmodule.createFieldConstant(-half_height)
+ nz = z / minus_c
+ zfact = fieldmodule.createFieldSqrt(value1 - nz * nz)
+ a = fieldmodule.createFieldConstant(half_ml_size) * zfact
+ b = fieldmodule.createFieldConstant(half_dv_size) * zfact
+ aa = a * a
+ bb = b * b
+ r = fieldmodule.createFieldSqrt((xx / aa) + (yy / bb))
+ rr = r * r
+ phi_r = value1 - rr
+ z_ext = -lower_lobe_extension * sin_pi__3
+ e = z / fieldmodule.createFieldConstant(z_ext)
+ ee = e * e
+ eee = ee * e
+ phi_e = fieldmodule.createFieldIf(fieldmodule.createFieldLessThan(e, value0), value0,
+ value3 * ee - value2 * eee)
+ # taper concavity to zero at oblique fissure
+ y_ext = fieldmodule.createFieldConstant(lower_lobe_extension * cos_pi__3)
+ y_max = b + y_ext
+ d = value1 + (y - e * y_ext) / y_max
+ dd = d * d
+ phi_y = value1 - dd
+ phi = phi_e * phi_r * phi_y
+ delta_y = phi * fieldmodule.createFieldConstant(-lower_lobe_base_concavity * cos_pi__3 / sin_pi__3)
+ delta_z = phi * fieldmodule.createFieldConstant(lower_lobe_base_concavity)
+ # get span a at displaced z, to calculate delta_x
+ displ_nz = (z + delta_z) / minus_c
+ displ_zfact = fieldmodule.createFieldSqrt(value1 - displ_nz * displ_nz)
+ displ_a = fieldmodule.createFieldConstant(half_ml_size) * displ_zfact
+ delta_x = (x * displ_a / a - x)
+ displacement = fieldmodule.createFieldConcatenate([delta_x, delta_y, delta_z])
+ new_coordinates = coordinates + displacement
+ fieldassignment = coordinates.createFieldassignment(new_coordinates)
+ fieldassignment.setNodeset(nodeset)
+ fieldassignment.assign()
+
+
+def taper_lung_edge(sharpeningFactor, fieldmodule, coordinates, nodeset, halfValue, isBase=False):
+ """
+ Apply a tapering transformation to the lung geometry to sharpen the anterior edge or the base.
+ If isBase is False, it sharpens the anterior edge (along the y-axis).
+ If isBase is True, it sharpens the base (along the z-axis), but only for nodes below a certain height.
+ :param sharpeningFactor: A value between 0 and 1, where 1 represents the maximum sharpness.
+ :param fieldmodule: Field module being worked with.
+ :param coordinates: The coordinate field. The anterior edge is towards the +y-axis, and the base is towards the
+ -z-axis.
+ :param nodeset: Zinc NodesetGroup containing nodes to transform.
+ :param halfValue: Half value of lung breadth/height depending on isBase.
+ :param isBase: False if transforming the anterior edge, True if transforming the base of the lung.
+ """
+ x = fieldmodule.createFieldComponent(coordinates, 1)
+ y = fieldmodule.createFieldComponent(coordinates, 2)
+ z = fieldmodule.createFieldComponent(coordinates, 3)
+
+ coord_value = z if isBase else y
+ start_value = 0.5 * halfValue if isBase else -0.5 * halfValue
+ end_value = -1.1 * halfValue if isBase else 1.1 * halfValue
+
+ start_value_field = fieldmodule.createFieldConstant(start_value)
+ end_value_field = fieldmodule.createFieldConstant(end_value)
+
+ xi = (coord_value - start_value_field) / fieldmodule.createFieldConstant(end_value - start_value)
+ xi__2 = xi * xi
+ one = fieldmodule.createFieldConstant(1.0)
+ x_scale = one - fieldmodule.createFieldConstant(sharpeningFactor) * xi__2
+ if isBase:
+ new_x = fieldmodule.createFieldIf(fieldmodule.createFieldLessThan(coord_value, start_value_field), x * x_scale, x)
+ else:
+ new_x = fieldmodule.createFieldIf(fieldmodule.createFieldGreaterThan(coord_value, start_value_field), x * x_scale, x)
+ new_coordinates = fieldmodule.createFieldConcatenate([new_x, y, z])
+
+ fieldassignment = coordinates.createFieldassignment(new_coordinates)
+ fieldassignment.setNodeset(nodeset)
+ fieldassignment.assign()
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py
index 94330329..2fa698f0 100644
--- a/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py
+++ b/src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py
@@ -1931,8 +1931,6 @@ def generateBaseMesh(cls, region, options):
get_uterus_term(allMarkers[i]), isMarker=True)
markerNode = group.createMarkerNode(nodeIdentifier, element=element_junction, xi=xi_junction)
nodeIdentifier = markerNode.getIdentifier() + 1
- for group in annotationGroups:
- group.getNodesetGroup(nodes).addNode(markerNode)
annotationGroups.remove(fundusSerosa)
annotationGroups.remove(cervixSerosa)
diff --git a/src/scaffoldmaker/meshtypes/scaffold_base.py b/src/scaffoldmaker/meshtypes/scaffold_base.py
index c6d0a7a8..ad540ffa 100644
--- a/src/scaffoldmaker/meshtypes/scaffold_base.py
+++ b/src/scaffoldmaker/meshtypes/scaffold_base.py
@@ -147,6 +147,8 @@ def generateMesh(cls, region, options):
if options.get('Refine'):
baseRegion = region.createRegion()
annotationGroups = cls.generateBaseMesh(baseRegion, options)[0]
+ # need faces to determine shared or boundary nodes during mesh refinement
+ baseRegion.getFieldmodule().defineAllFaces()
meshrefinement = MeshRefinement(baseRegion, region, annotationGroups)
cls.refineMesh(meshrefinement, options)
annotationGroups = meshrefinement.getAnnotationGroups()
diff --git a/src/scaffoldmaker/scaffoldpackage.py b/src/scaffoldmaker/scaffoldpackage.py
index fdc887dd..3fda757a 100644
--- a/src/scaffoldmaker/scaffoldpackage.py
+++ b/src/scaffoldmaker/scaffoldpackage.py
@@ -291,6 +291,8 @@ def deleteElementsInRanges(self, region, deleteElementRanges):
self._region = region
fm = self._region.getFieldmodule()
mesh = get_highest_dimension_mesh(fm)
+ if not mesh:
+ return
meshDimension = mesh.getDimension()
nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
with ChangeManager(fm):
diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py
index f6bd8a5b..ebebfb91 100644
--- a/src/scaffoldmaker/scaffolds.py
+++ b/src/scaffoldmaker/scaffolds.py
@@ -40,6 +40,7 @@
from scaffoldmaker.meshtypes.meshtype_3d_lung1 import MeshType_3d_lung1
from scaffoldmaker.meshtypes.meshtype_3d_lung2 import MeshType_3d_lung2
from scaffoldmaker.meshtypes.meshtype_3d_lung3 import MeshType_3d_lung3
+from scaffoldmaker.meshtypes.meshtype_3d_lung4 import MeshType_3d_lung4
from scaffoldmaker.meshtypes.meshtype_3d_musclefusiform1 import MeshType_3d_musclefusiform1
from scaffoldmaker.meshtypes.meshtype_3d_nerve1 import MeshType_3d_nerve1
from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1
@@ -107,6 +108,7 @@ class Scaffolds(object):
MeshType_3d_lung1,
MeshType_3d_lung2,
MeshType_3d_lung3,
+ MeshType_3d_lung4,
MeshType_3d_musclefusiform1,
MeshType_3d_nerve1,
MeshType_3d_ostium1,
diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py
index 2594eaac..a23408b5 100644
--- a/src/scaffoldmaker/utils/eft_utils.py
+++ b/src/scaffoldmaker/utils/eft_utils.py
@@ -701,20 +701,20 @@ def __init__(self):
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]),
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]])]
self._nodeLayout3WayPoints13 = [
- HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
- HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
+ HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [-1.0, 0.0, -1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
+ HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, -1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]])]
self._nodeLayout3WayPoints23 = [
- HermiteNodeLayout([[0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
- HermiteNodeLayout([[0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
+ HermiteNodeLayout([[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.0, -1.0, -1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
+ HermiteNodeLayout([[0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 1.0, -1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
HermiteNodeLayout([[0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
HermiteNodeLayout([[0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])]
self._nodeLayout4WayPoints = [
- HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]),
- HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, -1.0, 1.0]]),
- HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]]),
- HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]),
+ HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [-1.0, -1.0, -1.0]]),
+ HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, -1.0, -1.0]]),
+ HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [-1.0, 1.0, -1.0]]),
+ HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 1.0, -1.0]]),
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, -1.0, 1.0]]),
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]),
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]),
diff --git a/src/scaffoldmaker/utils/ellipsoidmesh.py b/src/scaffoldmaker/utils/ellipsoidmesh.py
index be6b307e..4503963f 100644
--- a/src/scaffoldmaker/utils/ellipsoidmesh.py
+++ b/src/scaffoldmaker/utils/ellipsoidmesh.py
@@ -1,7 +1,7 @@
"""
Utilities for building solid ellipsoid meshes from hexahedral elements.
"""
-from cmlibs.maths.vectorops import add, cross, div, magnitude, mult, set_magnitude, sub
+from cmlibs.maths.vectorops import add, cross, div, dot, magnitude, mult, normalize, rejection, set_magnitude, sub
from cmlibs.zinc.element import Element, Elementbasis
from cmlibs.zinc.field import Field
from cmlibs.zinc.node import Node
@@ -9,10 +9,9 @@
from scaffoldmaker.utils.eft_utils import determineCubicHermiteSerendipityEft, HermiteNodeLayoutManager
from scaffoldmaker.utils.geometry import (
getEllipsePointAtTrueAngle, getEllipseTangentAtPoint, moveCoordinatesToEllipsoidSurface,
- moveDerivativeToEllipsoidSurface, sampleCurveOnEllipsoid)
-from scaffoldmaker.utils.interpolation import (
- DerivativeScalingMode, get_nway_point, linearlyInterpolateVectors, sampleHermiteCurve,
- smoothCubicHermiteDerivativesLine)
+ moveDerivativeToEllipsoidSurface, moveDerivativeToEllipsoidSurfaceInPlane, sampleCurveOnEllipsoid)
+from scaffoldmaker.utils.interpolation import DerivativeScalingMode, linearlyInterpolateVectors, sampleHermiteCurve
+from scaffoldmaker.utils.hextetrahedronmesh import HexTetrahedronMesh
from scaffoldmaker.utils.quadtrianglemesh import QuadTriangleMesh
import copy
from enum import Enum
@@ -22,6 +21,8 @@
class EllipsoidSurfaceD3Mode(Enum):
SURFACE_NORMAL = 1 # surface D3 are exact surface normals to ellipsoid
OBLIQUE_DIRECTION = 2 # surface D3 are in direction of surface point on ellipsoid, gives flat oblique planes
+ # Surface D3 are surface normals to ellipsoid projected onto radial planes transitioning between axis2 and axis3
+ SURFACE_NORMAL_PLANE_PROJECTION = 3
class EllipsoidMesh:
@@ -29,16 +30,13 @@ class EllipsoidMesh:
Generates a solid ellipsoid of hexahedral elements with oblique cross axes suited to describing lung geometry.
"""
- def __init__(self, a, b, c, element_counts, transition_element_count,
- axis2_x_rotation_radians, axis3_x_rotation_radians, surface_only=False):
+ def __init__(self, a, b, c, element_counts, transition_element_count, surface_only=False):
"""
:param a: Axis length (radius) in x direction.
:param b: Axis length (radius) in y direction.
:param c: Axis length (radius) in z direction.
:param element_counts: Number of elements across full ellipse in a, b, c.
:param transition_element_count: Number of transition elements around outside >= 1.
- :param axis2_x_rotation_radians: Rotation of axis 2 about +x direction
- :param axis3_x_rotation_radians: Rotation of axis 3 about +x direction.
:param surface_only: Set to True to only make nodes and 2-D elements on the surface.
"""
assert all((count >= 4) and (count % 2 == 0) for count in element_counts)
@@ -48,8 +46,6 @@ def __init__(self, a, b, c, element_counts, transition_element_count,
self._c = c
self._element_counts = element_counts
self._trans_count = transition_element_count
- self._axis2_x_rotation_radians = axis2_x_rotation_radians
- self._axis3_x_rotation_radians = axis3_x_rotation_radians
self._surface_only = surface_only
self._nway_d_factor = 0.6
self._surface_d3_mode = EllipsoidSurfaceD3Mode.SURFACE_NORMAL
@@ -95,6 +91,8 @@ def __init__(self, a, b, c, element_counts, transition_element_count,
# print(s)
self._nx.append(nx_layer)
self._nids.append(nids_layer)
+ self._node_layout_manager = HermiteNodeLayoutManager()
+ self._prescribed_node_layouts = [] # list of (n1, n2, n3, node_layout)
def set_box_transition_groups(self, box_group, transition_group):
"""
@@ -132,177 +130,84 @@ def set_surface_d3_mode(self, surface_d3_mode: EllipsoidSurfaceD3Mode):
"""
self._surface_d3_mode = surface_d3_mode
- def build(self):
+ def build(self, axis2_x_rotation_radians, axis3_x_rotation_radians):
"""
- Determine coordinates and derivatives over and within ellipsoid.
+ Determine coordinates and derivatives over and within the full ellipsoid.
+ :param axis2_x_rotation_radians: Rotation of axis 2 about +x direction
+ :param axis3_x_rotation_radians: Rotation of axis 3 about +x direction.
"""
half_counts = [count // 2 for count in self._element_counts]
- box_counts = [half_counts[i] - self._trans_count for i in range(3)]
- octant1 = self._build_ellipsoid_octant(half_counts, self._axis2_x_rotation_radians, self._axis3_x_rotation_radians)
+ octant1 = self.build_octant(half_counts, axis2_x_rotation_radians, axis3_x_rotation_radians)
+ self.merge_octant(octant1, quadrant=0)
+ octant1.mirror_yz()
+ self.merge_octant(octant1, quadrant=2)
- # copy octant1 into ellipsoid
- octant_parameters = octant1.get_parameters()
- for n3 in range(half_counts[2] + 1):
- octant_nx_layer = octant_parameters[n3]
- nx_layer = self._nx[half_counts[2] + n3]
- for n2 in range(half_counts[1] + 1):
- octant_nx_row = octant_nx_layer[n2]
- nx_row = nx_layer[half_counts[1] + n2]
- for n1 in range(half_counts[0] + 1):
- nx_row[half_counts[0] + n1] = copy.deepcopy(octant_nx_row[n1])
+ octant2 = self.build_octant([half_counts[0], half_counts[2], half_counts[1]],
+ axis3_x_rotation_radians, axis2_x_rotation_radians + math.pi)
+ self.merge_octant(octant2, quadrant=1)
+ octant2.mirror_yz()
+ self.merge_octant(octant2, quadrant=3)
- octant2 = self._build_ellipsoid_octant([half_counts[0], half_counts[2], half_counts[1]],
- self._axis3_x_rotation_radians - math.pi, self._axis2_x_rotation_radians)
+ self.copy_to_negative_axis1()
- # transfer parameters on bottom plane of octant1 to target location of octant2 for blending
- n3 = half_counts[2]
- for i2 in range(1, half_counts[1] + 1):
- n2 = half_counts[1] + i2
- for i1 in range(half_counts[0] + 1):
- n1 = half_counts[0] + i1
- box = (i1 <= box_counts[0]) and (i2 <= box_counts[1])
- parameters = self._nx[n3][n2][n1]
- if not parameters or not parameters[0]:
- continue
- x, d1, d2, d3 = parameters
- x = [x[0], -x[1], -x[2]]
- d1 = [d1[0], -d1[1], -d1[2]] if box else [-d1[0], d1[1], d1[2]]
- d2 = [-d2[0], d2[1], d2[2]]
- d3 = ([-d3[0], d3[1], d3[2]] if box else [d3[0], -d3[1], -d3[2]]) if d3 else None
- self._nx[n3][self._element_counts[1] - n2][n1] = [x, d1, d2, d3]
-
- # copy and mirror in y and z octant2 into ellipsoid, blending existing derivatives
- octant_parameters = octant2.get_parameters()
- for o3 in range(half_counts[1] + 1):
- octant_nx_layer = octant_parameters[o3]
- for o2 in range(half_counts[2] + 1):
- octant_nx_row = octant_nx_layer[o2]
- nx_row = self._nx[half_counts[2] + o2][half_counts[1] - o3]
- box_row = (o2 <= box_counts[2]) and (o3 <= box_counts[1])
- top_transition = o2 > box_counts[2]
- for o1 in range(half_counts[0] + 1):
- octant_nx = octant_nx_row[o1]
- if octant_nx and octant_nx[0]:
- box = box_row and (o1 <= box_counts[0])
- x, d1, d2, d3 = octant_nx
- x = [x[0], -x[1], -x[2]]
- if top_transition:
- if o3 > box_counts[1]:
- d1 = [d1[0], -d1[1], -d1[2]]
- d2 = [d2[0], -d2[1], -d2[2]]
- # fix 3-way point case on transition 'corner':
- if (o3 >= box_counts[1]) and (o1 >= box_counts[0]):
- d2 = add(d1, d2)
- else:
- d1 = [-d1[0], d1[1], d1[2]]
- d2 = [-d2[0], d2[1], d2[2]]
- d3 = [d3[0], -d3[1], -d3[2]] if d3 else None
- elif box:
- d1 = [d1[0], -d1[1], -d1[2]]
- d2, d3 = [-d3[0], d3[1], d3[2]], [d2[0], -d2[1], -d2[2]]
- else:
- if o3 > box_counts[1]:
- d1 = [d1[0], -d1[1], -d1[2]]
- d2 = [d2[0], -d2[1], -d2[2]]
- else:
- d1, d2 = [-d2[0], d2[1], d2[2]], [d1[0], -d1[1], -d1[2]]
- d3 = [d3[0], -d3[1], -d3[2]] if d3 else None
- new_nx = [x, d1, d2, d3]
- nx = nx_row[half_counts[0] + o1]
- # expect coordinates to be at the same location on boundaries
- nx[0] = copy.copy(x)
- for pix in range(1, 4):
- d = new_nx[pix]
- if nx[pix] and d:
- # blend derivatives with harmonic mean magnitude; should already be in same direction
- d = linearlyInterpolateVectors(
- nx[pix], d, 0.5, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)
- nx[pix] = copy.copy(d)
-
- # transfer blended d2 derivatives on bottom plane of octant2 back to octant1
- n3 = half_counts[2]
- for i2 in range(1, half_counts[1] + 1):
- n2 = half_counts[1] + i2
- for i1 in range(half_counts[0] + 1):
- n1 = half_counts[0] + i1
- box = (i1 <= box_counts[0]) and (i2 <= box_counts[1])
- parameters = self._nx[n3][self._element_counts[1] - n2][n1]
- if not parameters or not parameters[0]:
- continue
- x, d1, d2, d3 = parameters
- x = [x[0], -x[1], -x[2]]
- d1 = [d1[0], -d1[1], -d1[2]] if box else [-d1[0], d1[1], d1[2]]
- d2 = [-d2[0], d2[1], d2[2]]
- d3 = ([-d3[0], d3[1], d3[2]] if box else [d3[0], -d3[1], -d3[2]]) if d3 else None
- self._nx[n3][n2][n1] = [x, d1, d2, d3]
-
- # mirror octants over x = 0
- for n3 in range(half_counts[2], self._element_counts[2] + 1):
- for n2 in range(0, self._element_counts[1] + 1):
- for n1 in range(half_counts[0] + 1, self._element_counts[0] + 1):
- parameters = self._nx[n3][n2][n1]
- if not parameters or not parameters[0]:
- continue
- x, d1, d2, d3 = parameters
- x = [-x[0], x[1], x[2]]
- d1 = [d1[0], -d1[1], -d1[2]] if d1 else None
- d2 = [-d2[0], d2[1], d2[2]] if d2 else None
- d3 = [-d3[0], d3[1], d3[2]] if d3 else None
- self._nx[n3][n2][self._element_counts[0] - n1] = [x, d1, d2, d3]
-
- # flip top half about both y = 0 and z = 0
- for n3 in range(half_counts[2] + 1, self._element_counts[2] + 1):
- for n2 in range(0, self._element_counts[1] + 1):
- top_nx_row = self._nx[n3][n2]
- box2 = (self._trans_count <= n2 <=
- (self._element_counts[1] - self._trans_count))
- box_row = (n3 <= (half_counts[2] + box_counts[2])) and box2
- bottom_nx_row = self._nx[self._element_counts[2] - n3][self._element_counts[1] - n2]
- bottom_transition = box2 and (n3 >= (half_counts[2] + box_counts[2]))
- for n1 in range(self._element_counts[0] + 1):
- parameters = top_nx_row[n1]
- if parameters and parameters[0]:
- box = box_row and (self._trans_count <= n1 <=
- (self._element_counts[0] - self._trans_count))
- x, d1, d2, d3 = parameters
- x = [x[0], -x[1], -x[2]]
- d2 = [-d2[0], d2[1], d2[2]] if d2 else None
- if box and not bottom_transition:
- d1 = [d1[0], -d1[1], -d1[2]]
- d3 = [-d3[0], d3[1], d3[2]]
- else:
- d1 = [-d1[0], d1[1], d1[2]] if d1 else None
- d3 = [d3[0], -d3[1], -d3[2]] if d3 else None
- bottom_nx_row[n1] = [x, d1, d2, d3]
-
- def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x_rotation_radians):
+ def build_octant(self, half_counts, axis2_x_rotation_radians, axis3_x_rotation_radians,
+ axis2_extension=0.0, axis2_extension_elements_count=0):
"""
Get coordinates of top, right, front octant with supplied angles.
:param half_counts: Numbers of elements across octant 1, 2 and 3 directions.
:param axis2_x_rotation_radians: Rotation of axis 2 about +x direction
:param axis3_x_rotation_radians: Rotation of axis 3 about +x direction.
+ :param axis2_extension: Extension distance along axis2 beyond origin [0.0, 0.0, 0.0].
+ :param axis2_extension_elements_count: If axis2_extension: number of elements beyond origin.
+ Note: included in half_counts[1].
+ :return: HexTetrahedronMesh
"""
- elements_count_q12 = half_counts[0] + half_counts[1] - 2 * self._trans_count
- elements_count_q13 = half_counts[0] + half_counts[2] - 2 * self._trans_count
- elements_count_q23 = half_counts[1] + half_counts[2] - 2 * self._trans_count
+ assert ((axis2_extension == 0.0) and (axis2_extension_elements_count == 0)) or (
+ (axis2_extension > 0.0) and (0 < axis2_extension_elements_count))
box_counts = [half_counts[i] - self._trans_count for i in range(3)]
+ cos_axis2 = math.cos(axis2_x_rotation_radians)
+ sin_axis2 = math.sin(axis2_x_rotation_radians)
+ cos_axis3 = math.cos(axis3_x_rotation_radians)
+ sin_axis3 = math.sin(axis3_x_rotation_radians)
+
origin = [0.0, 0.0, 0.0]
- axis1 = [self._a, 0.0, 0.0]
+ ext_origin = [0.0, -axis2_extension * cos_axis2, -axis2_extension * sin_axis2]
+ ext_axis1 = axis1 = [self._a, 0.0, 0.0]
axis2 = [0.0] + getEllipsePointAtTrueAngle(self._b, self._c, axis2_x_rotation_radians)
- axis3 = [0.0] + getEllipsePointAtTrueAngle(self._b, self._c, axis3_x_rotation_radians)
+ axis2_mag = magnitude(axis2)
+ axis2_normal = normalize([0.0, axis2[2], -axis2[1]])
+ ext_axis3 = axis3 = [0.0] + getEllipsePointAtTrueAngle(self._b, self._c, axis3_x_rotation_radians)
+ axis3_normal = normalize([0.0, axis3[2], -axis3[1]])
+ if axis2_extension_elements_count:
+ assert axis2_extension < axis2_mag # extension must not go outside ellipsoid
+ xb, xa = getEllipsePointAtTrueAngle(axis2_mag, self._a, math.pi / 2.0, [-axis2_extension, 0.0])
+ ext_axis1 = [xa, xb * cos_axis2, xb * sin_axis2]
+ ext_axis3 = [0.0] + getEllipsePointAtTrueAngle(self._b, self._c, axis3_x_rotation_radians, ext_axis1[1:])
+ ext_axis3m = [0.0] + getEllipsePointAtTrueAngle(self._b, self._c, axis3_x_rotation_radians + math.pi, ext_axis1[1:])
+ centre_mod_axis3 = mult(add(ext_axis3, ext_axis3m), 0.5)
+ mod_axis3 = sub(ext_axis3, centre_mod_axis3)
+ mag_mod_axis3 = magnitude(mod_axis3)
+ else:
+ centre_mod_axis3 = origin
+ mag_mod_axis3 = magnitude(axis3)
+
axis_d1 = div(axis1, half_counts[0])
+ ext_axis_d1 = div(sub(ext_axis1, ext_origin), half_counts[0])
axis_d2 = div(axis2, half_counts[1])
axis_d3 = div(axis3, half_counts[2])
+ ext_axis_d3 = div(sub(ext_axis3, ext_origin), half_counts[2])
axis_md1 = [-d for d in axis_d1]
axis_md2 = [-d for d in axis_d2]
axis_md3 = [-d for d in axis_d3]
- # magnitude of most derivatives indicating only direction so magnitude not known
+ ext_axis_md1 = [-d for d in ext_axis_d1]
+
+ # most derivatives indicate only direction, so magnitude not known
dir_mag = min(magnitude(axis_d1), magnitude(axis_d2), magnitude(axis_d3))
axis2_dt = set_magnitude([0.0] + getEllipseTangentAtPoint(self._b, self._c, axis2[1:]), magnitude(axis_d3))
- axis3_dt = set_magnitude([0.0] + getEllipseTangentAtPoint(self._b, self._c, axis3[1:]), magnitude(axis_d2))
- axis3_mdt = [-d for d in axis3_dt]
+ ext_axis3_dt = set_magnitude([0.0] + getEllipseTangentAtPoint(self._b, self._c, ext_axis3[1:]), magnitude(ext_axis_d3))
+ ext_axis3_mdt = [-d for d in ext_axis3_dt]
axis2_mag = magnitude(axis2)
axis3_mag = magnitude(axis3)
@@ -314,33 +219,80 @@ def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x
start_weight, end_weight, overweighting, end_transition))
move_x_to_ellipsoid_surface = lambda x: moveCoordinatesToEllipsoidSurface(self._a, self._b, self._c, x)
move_d_to_ellipsoid_surface = lambda x, d: moveDerivativeToEllipsoidSurface(self._a, self._b, self._c, x, d)
+ def evaluate_surface_d3_ellipsoid_plane(tx, td1, td2):
+ """
+ Restrict d3 to be the ellipsoid normal constrained to be in radial planes from ext_origin through tx,
+ varying between axis_d2 and axis_d3.
+ :param tx: Coordinates of a point on the ellipsoid surface in the octant.
+ :param td1: Unused point d1.
+ :param td2: Unused point d2.
+ :return: Radial plane constrained ellipsoid normal d3 with magnitude dir_mag.
+ """
+ n = [tx[0] / (self._a * self._a), tx[1] / (self._b * self._b), tx[2] / (self._c * self._c)]
+ if dot(tx, axis3_normal) <= 1.0E-5:
+ if dot(tx, axis2_normal) >= -1.0E-5:
+ return set_magnitude(axis1, dir_mag)
+ else:
+ plane_normal = [0.0, axis3[2], -axis3[1]]
+ else:
+ plane_normal = [0.0, tx[2], -tx[1]]
+ normal = rejection(n, plane_normal)
+ return set_magnitude(normal, dir_mag)
if self._surface_d3_mode == EllipsoidSurfaceD3Mode.SURFACE_NORMAL:
evaluate_surface_d3_ellipsoid = lambda tx, td1, td2: set_magnitude(
[tx[0] / (self._a * self._a), tx[1] / (self._b * self._b), tx[2] / (self._c * self._c)], dir_mag)
- else:
- evaluate_surface_d3_ellipsoid = lambda tx, td1, td2: set_magnitude(tx, dir_mag)
-
- octant = EllipsoidOctantMesh(self._a, self._b, self._c, half_counts, self._trans_count,
- nway_d_factor=self._nway_d_factor)
+ elif self._surface_d3_mode == EllipsoidSurfaceD3Mode.OBLIQUE_DIRECTION:
+ evaluate_surface_d3_ellipsoid=lambda tx, td1, td2: set_magnitude(tx, dir_mag)
+ else: # EllipsoidSurfaceD3Mode.SURFACE_NORMAL_PLANE_PROJECTION
+ evaluate_surface_d3_ellipsoid = evaluate_surface_d3_ellipsoid_plane
+
+ ext_half_counts = [
+ half_counts[0],
+ half_counts[1] + axis2_extension_elements_count,
+ half_counts[2]
+ ]
+ diag_counts = [
+ half_counts[0] + half_counts[1] - 2 * self._trans_count,
+ half_counts[0] + half_counts[2] - 2 * self._trans_count,
+ half_counts[1] + half_counts[2] - 2 * self._trans_count
+ ]
+ ext_diag_counts = [
+ ext_half_counts[0] + ext_half_counts[1] - 2 * self._trans_count,
+ ext_half_counts[0] + ext_half_counts[2] - 2 * self._trans_count,
+ ext_half_counts[1] + ext_half_counts[2] - 2 * self._trans_count
+ ]
+ octant = HexTetrahedronMesh(ext_half_counts, ext_diag_counts, nway_d_factor=self._nway_d_factor)
# get outside curve from axis 1 to axis 2
abx, abd1, abd2 = sampleCurveOnEllipsoid(
self._a, self._b, self._c,
axis1, axis_d2, axis_d3,
axis2, axis_md1, axis2_dt,
- elements_count_q12)
+ diag_counts[0])
+ if axis2_extension_elements_count:
+ end_axis_d3 = moveDerivativeToEllipsoidSurfaceInPlane(
+ self._a, self._b, self._c, ext_axis1, [axis_d3[0], -axis_d3[2], axis_d3[1]], axis_d3)
+ ext_abx, ext_abd1, ext_abd2 = sampleCurveOnEllipsoid(
+ self._a, self._b, self._c,
+ axis1, [-d for d in axis_d2], axis_d3,
+ ext_axis1, None, end_axis_d3, # axis_d3
+ axis2_extension_elements_count)
+ for i in range(1, axis2_extension_elements_count + 1):
+ abx.insert(0, ext_abx[i])
+ abd1.insert(0, [-d for d in ext_abd1[i]])
+ abd2.insert(0, ext_abd2[i])
# get outside curve from axis 1 to axis 3
acx, acd2, acd1 = sampleCurveOnEllipsoid(
self._a, self._b, self._c,
abx[0], abd2[0], abd1[0],
- axis3, axis_md1, axis3_mdt,
- elements_count_q13)
+ ext_axis3, ext_axis_md1, ext_axis3_mdt,
+ ext_diag_counts[1])
# get outside curve from axis 2 to axis 3
bcx, bcd2, bcd1 = sampleCurveOnEllipsoid(
self._a, self._b, self._c,
abx[-1], abd2[-1], abd1[-1],
acx[-1], [-d for d in acd1[-1]], acd2[-1],
- elements_count_q23)
+ ext_diag_counts[2])
# fix first/last derivatives
abd2[0] = acd2[0]
abd2[-1] = bcd2[0]
@@ -348,7 +300,7 @@ def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x
# make outer surface triangle of octant 1
triangle_abc = QuadTriangleMesh(
- box_counts[0], box_counts[1], box_counts[2],
+ box_counts[0], box_counts[1] + axis2_extension_elements_count, box_counts[2],
sample_curve_on_ellipsoid, move_x_to_ellipsoid_surface, move_d_to_ellipsoid_surface, self._nway_d_factor)
triangle_abc.set_edge_parameters12(abx, abd1, abd2)
triangle_abc.set_edge_parameters13(acx, acd1, acd2)
@@ -366,25 +318,49 @@ def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x
# build interior lines from axis1, axis2, axis3 to origin
aox, aod2, aod1 = sampleHermiteCurve(
- axis1, axis_md1, abd1[0], origin, axis_md1, axis_d2, elements_count=half_counts[0])
+ ext_axis1, ext_axis_md1, abd1[0], ext_origin, ext_axis_md1, axis_d2, elements_count=half_counts[0])
box, bod2, bod3 = sampleHermiteCurve(
axis2, axis_md2, bcd2[0], origin, axis_md2, axis_d3, elements_count=half_counts[1])
+ if axis2_extension_elements_count:
+ ext_box, ext_bod2, ext_bod3 = sampleHermiteCurve(
+ box[-1], bod2[-1], bod3[-1], ext_origin, None, ext_axis_d3,
+ elements_count=axis2_extension_elements_count)
+ for i in range(1, axis2_extension_elements_count + 1):
+ box.append(ext_box[i])
+ bod2.append(ext_bod2[i])
+ bod3.append(ext_bod3[i])
bod1 = [abd1[-1]] + [axis_md1] * (len(box) - 1)
cox, cod2, cod1 = sampleHermiteCurve(
- axis3, axis_md3, [-d for d in acd1[-1]], origin, axis_md3, axis_md2, elements_count=half_counts[2])
+ ext_axis3, axis_md3, [-d for d in acd1[-1]], ext_origin, axis_md3, bod2[-1],
+ elements_count=half_counts[2])
# make inner surface triangle 1-2-origin
triangle_abo = QuadTriangleMesh(
- box_counts[0], box_counts[1], self._trans_count, sampleHermiteCurve,
+ box_counts[0], box_counts[1] + axis2_extension_elements_count, self._trans_count, sampleHermiteCurve,
nway_d_factor=self._nway_d_factor)
abd3 = [[-d for d in evaluate_surface_d3_ellipsoid(x, None, None)] for x in abx]
triangle_abo.set_edge_parameters12(abx, abd1, abd3, abd2)
- aod3 = [abd2[0]] + [axis_d3] * (len(aox) - 1)
+ count = len(aox) - 1
+ aod3 = [linearlyInterpolateVectors(abd2[0], axis_d3, i / count) for i in range(count + 1)]
+
triangle_abo.set_edge_parameters13(aox, aod1, aod2, aod3)
triangle_abo.set_edge_parameters23(box, bod1, bod2, bod3)
- triangle_abo.build()
- triangle_abo.assign_d3(lambda tx, td1, td2:
- linearlyInterpolateVectors(axis_d3, axis2_dt, magnitude(tx[1:]) / axis2_mag))
+ triangle_abo.build(regular_count2=axis2_extension_elements_count)
+ aa = self._a * self._a
+ bb = axis2_mag * axis2_mag # of axis2 ellipse
+ def evaluate_surface_d3_abo(tx, td1, td2):
+ y = tx[1] * cos_axis2 + tx[2] * sin_axis2
+ yy = y * y
+ if yy >= bb:
+ return axis2_dt # tip of ellipse
+ centre_d3 = (linearlyInterpolateVectors(axis_d3, axis2_dt, magnitude(tx[1:]) / axis2_mag)
+ if (y >= 0.0) else axis_d3)
+ xx = aa * (1.0 - yy / bb)
+ x = math.sqrt(xx)
+ side_x = [x, tx[1], tx[2]]
+ side_d3 = moveDerivativeToEllipsoidSurface(self._a, self._b, self._c, side_x, centre_d3)
+ return linearlyInterpolateVectors(centre_d3, side_d3, tx[0] / x)
+ triangle_abo.assign_d3(evaluate_surface_d3_abo)
octant.set_triangle_abo(triangle_abo)
# extract exact derivatives
aod1 = triangle_abo.get_edge_parameters13()[1]
@@ -402,15 +378,29 @@ def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x
cod3 = [acd2[-1]] + [axis_md1] * (len(cox) - 1)
triangle_aco.set_edge_parameters23(cox, cod3, cod2, cod1)
triangle_aco.build()
- triangle_aco.assign_d3(lambda tx, td1, td2:
- linearlyInterpolateVectors(axis_md2, axis3_dt, magnitude(tx[1:]) / axis3_mag))
+ aa = self._a * self._a
+ bb = mag_mod_axis3 * mag_mod_axis3 # mod_axis3 ellipse
+ def evaluate_surface_d3_aco(tx, td1, td2):
+ mx = sub(tx, centre_mod_axis3)
+ y = mx[1] * cos_axis3 + mx[2] * sin_axis3
+ yy = y * y
+ if yy >= bb:
+ return ext_axis3_dt # tip of ellipse
+ centre_d3 = (linearlyInterpolateVectors(axis_md2, ext_axis3_dt, magnitude(mx[1:]) / axis3_mag)
+ if (y >= 0.0) else axis_md2)
+ xx = aa * (1.0 - yy / bb)
+ x = math.sqrt(xx)
+ side_x = [x, tx[1], tx[2]]
+ side_d3 = moveDerivativeToEllipsoidSurface(self._a, self._b, self._c, side_x, centre_d3)
+ return linearlyInterpolateVectors(centre_d3, side_d3, tx[0] / x)
+ triangle_aco.assign_d3(evaluate_surface_d3_aco)
octant.set_triangle_aco(triangle_aco)
# extract exact derivatives
cod3 = triangle_aco.get_edge_parameters23()[1]
# make inner surface 2-3-origin
triangle_bco = QuadTriangleMesh(
- box_counts[1], box_counts[2], self._trans_count, sampleHermiteCurve,
+ box_counts[1] + axis2_extension_elements_count, box_counts[2], self._trans_count, sampleHermiteCurve,
nway_d_factor=self._nway_d_factor)
bcd3 = [bod2[0]] + [[-d for d in evaluate_surface_d3_ellipsoid(x, None, None)] for x in bcx[1:-1]] \
+ [cod2[-1]]
@@ -428,158 +418,121 @@ def _build_ellipsoid_octant(self, half_counts, axis2_x_rotation_radians, axis3_x
return octant
- def _next_increment_out_of_bounds(self, indexes, index_increment):
- for c in range(3):
- index = indexes[c] + index_increment[c]
- if (index < 0) or (index > self._element_counts[c]):
- return True
- return False
-
- def _set_coordinates_around_rim(self, parameters, parameter_indexes, start_indexes, index_increments,
- skip_start=False, skip_end=False,
- blend_start=False, blend_middle=False, blend_end=False):
+ def merge_octant(self, octant: HexTetrahedronMesh, quadrant: int):
"""
- Insert parameters around the rim into the coordinates array.
- :param parameters: List of lists of N node parameters e.g. [px, pd1, pd2]
- :param parameter_indexes: Lists of parameter indexes where x=0, d1=1, d2=2, d3=3. Starts with first and
- advances after each corner, then cycles back to first. Can be negative to invert vector.
- e.g. [[0, 1, 2], [0, -2, 1]] for [x, d1, d2] then [x, -d2, d1] after first corner.
- :param start_indexes: Index of first point.
- :param index_increments: List of increments in indexes. Starts with first and after at each corner, then
- cycles back to first.
- :param skip_start: Set to True to skip the first value.
- :param skip_end: Set to True to skip the last value.
- :param blend_start: Set to True to blend parameters with any old parameters at start location.
- :param blend_middle: Set to True to blend parameters with any old parameters at middle locations.
- :param blend_end: Set to True to blend parameters with any old parameters at end location.
+ Merge octant parameters into ellipsoid in one of the 4 quadrants on +axis1.
+ Octant can be extended into its axis2 over regular elements of ellipsoid.
+ :param octant: HexTetrahedronMesh
+ :param quadrant: 0 for +axis2, +axis3 increasing anticlockwise around +axis1 up to 3.
"""
- indexes = start_indexes
- parameter_number = 0
- parameter_index = parameter_indexes[0]
- increment_number = 0
- index_increment = index_increments[0]
- start_n = 1 if skip_start else 0
- last_n = len(parameters[0]) - 1
- limit_n = len(parameters[0]) - (1 if skip_end else 0)
- for n in range(start_n, limit_n):
- if n > 0:
- while True:
- indexes = [indexes[c] + index_increment[c] for c in range(3)]
- # skip over blank transition coordinates
- if self._nx[indexes[2]][indexes[1]][indexes[0]]:
- break
- if self._next_increment_out_of_bounds(indexes, index_increment):
- parameter_number += 1
- if parameter_number == len(parameter_indexes):
- parameter_number = 0
- parameter_index = parameter_indexes[parameter_number]
- increment_number += 1
- if increment_number == len(index_increments):
- increment_number = 0
- index_increment = index_increments[increment_number]
- nx = self._nx[indexes[2]][indexes[1]][indexes[0]]
- for parameter, spix in zip(parameters, parameter_index):
- new_parameter = [-d for d in parameter[n]] if (spix < 0) else copy.copy(parameter[n])
- pix = abs(spix)
- if nx[pix] and ((blend_start and (n == 0)) or (blend_middle and (0 < n < last_n)) or
- (blend_end and (n == last_n))):
- if pix == 0:
- # for fairness, move to surface before blending
- new_parameter = moveCoordinatesToEllipsoidSurface(self._a, self._b, self._c, new_parameter)
- new_parameter = [0.5 * (nx[pix][c] + new_parameter[c]) for c in range(3)]
- else:
- # harmonic mean to cope with significant element size differences on boundary
- new_parameter = linearlyInterpolateVectors(
- nx[pix], new_parameter, 0.5, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)
- nx[pix] = new_parameter
+ assert 0 <= quadrant <= 3
+ half_counts = [count // 2 for count in self._element_counts]
+ axis_counts = octant.get_axis_counts()
+ even_quadrant = quadrant in (0, 2)
+ assert half_counts[0] == axis_counts[0]
+ ext_count2 = (axis_counts[1] - half_counts[1]) if even_quadrant else 0
+ ext_count3 = 0 if even_quadrant else (axis_counts[1] - half_counts[2])
+ assert 0 <= ext_count2 < (half_counts[1] - self._trans_count)
+ assert 0 <= ext_count3 < (half_counts[2] - self._trans_count)
+ obox_counts = octant.get_box_counts()
+ # box_counts = [half_counts[i] - self._trans_count for i in range(3)]
+ octant_parameters = octant.get_parameters()
+
+ for o3 in range(axis_counts[2] + 1):
+ for o2 in range(axis_counts[1] + 1):
+ ox_row = octant_parameters[o3][o2]
+ if even_quadrant:
+ if quadrant == 0:
+ n3 = half_counts[2] + o3 - ext_count3
+ n2 = half_counts[1] + o2 - ext_count2
+ else: # quadrant == 2:
+ n3 = half_counts[2] - o3 + ext_count3
+ n2 = half_counts[1] - o2 + ext_count2
+ else:
+ if quadrant == 1:
+ n3 = half_counts[2] + o2 - ext_count3
+ n2 = half_counts[1] - o3 + ext_count2
+ else: # if quadrant == 3:
+ n3 = half_counts[2] - o2 + ext_count3
+ n2 = half_counts[1] + o3 - ext_count2
+ transition3 = (n3 < self._trans_count) or (n3 > (self._element_counts[2] - self._trans_count))
+ transition2 = (n2 < self._trans_count) or (n2 > (self._element_counts[1] - self._trans_count))
+ # bottom_transition = n3 < self._trans_count
+ nx_row = self._nx[n3][n2]
+ obox_row = (o3 <= obox_counts[2]) and (o2 <= obox_counts[1])
+ for o1 in range(axis_counts[0] + 1):
+ n1 = half_counts[0] + o1
+ transition1 = n1 > (self._element_counts[0] - self._trans_count)
+ obox = obox_row and (o1 <= obox_counts[0])
+ ox = ox_row[o1]
+ if ox and ox[0]:
+ x = copy.copy(ox[0])
+ perm = None
+ if even_quadrant:
+ if quadrant == 0:
+ perm = [1, 2, 3]
+ else: # quadrant == 2:
+ perm = [1, -2, -3] if obox else [-1, -2, 3]
+ else:
+ if obox:
+ perm = [1, -3, 2]
+ elif transition3:
+ if transition2:
+ if transition1:
+ # fix 3-way point
+ ox = [ox[0], ox[1], add(ox[1], ox[2]), ox[3]]
+ perm = [1, 2, 3]
+ else:
+ perm = [-1, -2, 3]
+ else:
+ if transition1 and not transition2:
+ perm = [-2, 1, 3]
+ else:
+ perm = [1, 2, 3]
+ if quadrant == 3:
+ perm = [perm[0], -perm[1], -perm[2]] if obox else [-perm[0], -perm[1], perm[2]]
+ d1, d2, d3 = [copy.copy(ox[i]) if (i > 0) else [-d for d in ox[-i]] for i in perm]
+ new_nx = [x, d1, d2, d3]
+ # merge:
+ nx = nx_row[n1]
+ for i in range(4):
+ d = new_nx[i]
+ if i and nx[i]:
+ # blend derivatives with harmonic mean magnitude; should already be in same direction
+ d = linearlyInterpolateVectors(
+ nx[i], d, 0.5, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)
+ nx[i] = d
- def _smooth_derivatives_around_rim(self, start_indexes, end_indexes, index_increments,
- derivative_indexes, end_derivative_index,
- fix_start_direction=True, fix_end_direction=True,
- blend_start=False, blend_end=False):
+ def copy_to_negative_axis1(self):
"""
- Smooth derivatives around the rim in the coordinates array.
- :param start_indexes: Indexes of first point.
- :param end_indexes: Indexes of last point.
- :param index_increments: List of increments in indexes. Starts with first and after at each corner, then
- cycles back to first.
- :param derivative_indexes: List of signed derivative parameter index to along where d1=1, d2=2, d3=3.
- Starts with first and advances after each corner, then cycles back to first. Can be negative to invert vector.
- e.g. [1, -2] for d1 then -d2 after first corner.
- :param end_derivative_index: List of signed derivative indexes to apply on the last point
- e.g. [1, -2] gives d1 - d2.
- :param fix_start_direction: Set to True to keep the start direction but scale its magnitude.
- :param fix_end_direction: Set to True to keep the end direction but scale its magnitude.
- :param blend_start: Set to True to 50:50 blend parameters with any old parameters at start location.
- :param blend_end: Set to True to 50:50 blend parameters with any old parameters at end location.
+ Copy parameters from +axis1 to -axis1.
"""
- indexes = start_indexes
- derivative_number = 0
- derivative_index = derivative_indexes[0]
- increment_number = 0
- index_increment = index_increments[0]
- indexes_list = []
- derivative_index_list = []
- px = []
- pd = []
- n = 0
- while True:
- if n > 0:
- if indexes == end_indexes:
- break
- while True:
- indexes = [indexes[c] + index_increment[c] for c in range(3)]
- # skip over blank transition coordinates
- if self._nx[indexes[2]][indexes[1]][indexes[0]]:
- break
- if self._next_increment_out_of_bounds(indexes, index_increment):
- derivative_number += 1
- if derivative_number == len(derivative_indexes):
- derivative_number = 0
- derivative_index = derivative_indexes[derivative_number]
- increment_number += 1
- if increment_number == len(index_increments):
- increment_number = 0
- index_increment = index_increments[increment_number]
- parameters = self._nx[indexes[2]][indexes[1]][indexes[0]]
- x = parameters[0]
- use_derivative_index = end_derivative_index if (indexes == end_indexes) else [derivative_index]
- indexes_list.append(copy.copy(indexes))
- derivative_index_list.append(copy.copy(use_derivative_index))
- d = [0.0, 0.0, 0.0]
- for i in range(len(use_derivative_index)):
- spix = use_derivative_index[i]
- pix = abs(spix)
- values = parameters[pix]
- if values:
- if spix < 0:
- values = [-ad for ad in values]
- d = add(d, values)
- px.append(x)
- pd.append(d)
- n += 1
- sd = smoothCubicHermiteDerivativesLine(
- px, pd, fixStartDirection=fix_start_direction, fixEndDirection=fix_end_direction,
- fixEndDerivative=len(end_derivative_index) > 1)
- for n in range(len(sd)):
- sd[n] = moveDerivativeToEllipsoidSurface(self._a, self._b, self._c, px[n], sd[n])
- sd = smoothCubicHermiteDerivativesLine(
- px, sd, fixAllDirections=True, fixEndDerivative=len(end_derivative_index) > 1)
- last_n = len(sd) - 1
- for n in range(len(sd)):
- indexes = indexes_list[n]
- derivative_index = derivative_index_list[n]
- parameters = self._nx[indexes[2]][indexes[1]][indexes[0]]
- if len(derivative_index) == 1:
- spix = derivative_index[0]
- new_derivative = [-d for d in sd[n]] if (spix < 0) else sd[n]
- pix = abs(spix)
- if parameters[pix] and (blend_start and (n == 0)) or (blend_end and (n == last_n)):
- new_derivative = linearlyInterpolateVectors(parameters[pix], new_derivative, 0.5)
- parameters[pix] = new_derivative
- # else:
- # # not putting back values if summed parameters
+ half_counts0 = self._element_counts[0] // 2
+ for n3 in range(self._element_counts[2] + 1):
+ nx_layer = self._nx[n3]
+ for n2 in range(self._element_counts[1] + 1):
+ nx_row = nx_layer[n2]
+ for n1 in range(half_counts0 + 1, self._element_counts[0] + 1):
+ nx = self._nx[n3][n2][n1]
+ if nx and nx[0]:
+ x, d1, d2, d3 = nx
+ x = [-x[0], x[1], x[2]]
+ d1 = [d1[0], -d1[1], -d1[2]] if d1 else None
+ d2 = [-d2[0], d2[1], d2[2]] if d2 else None
+ d3 = [-d3[0], d3[1], d3[2]] if d3 else None
+ nx_row[self._element_counts[0] - n1] = [x, d1, d2, d3]
+ def _next_increment_out_of_bounds(self, indexes, index_increment):
+ """
+ :param indexes: List of 3 current indexes in the EllipsoidMesh nodes block: [n1, n2, n3].
+ :param index_increment: Increments to the 3 indexes.
+ :return: True if adding index_increment to indexes puts the index outside the nodes block.
+ """
+ for c in range(3):
+ index = indexes[c] + index_increment[c]
+ if (index < 0) or (index > self._element_counts[c]):
+ return True
+ return False
def _get_nid_to_node_layout_map_3d(self, node_layout_manager):
"""
@@ -688,8 +641,37 @@ def _get_nid_to_node_layout_map_3d(self, node_layout_manager):
nid_to_node_layout[nid] = node_layout_3way12[2]
nid = self._nids[self._element_counts[2] - nt][self._element_counts[1] - nt][self._element_counts[0] - nt]
nid_to_node_layout[nid] = node_layout_3way12[3]
+ # add prescribed node layouts
+ for n1, n2, n3, node_layout in self._prescribed_node_layouts:
+ nid = self._nids[n3][n2][n1]
+ nid_to_node_layout[nid] = node_layout
return nid_to_node_layout
+ def get_node_layout_manager(self):
+ return self._node_layout_manager
+
+ def get_node_identifier(self, n1, n2, n3):
+ assert 0 <= n1 <= self._element_counts[0]
+ assert 0 <= n2 <= self._element_counts[1]
+ assert 0 <= n3 <= self._element_counts[2]
+ return self._nids[n3][n2][n1]
+
+ def get_node_parameters(self, n1, n2, n3):
+ assert 0 <= n1 <= self._element_counts[0]
+ assert 0 <= n2 <= self._element_counts[1]
+ assert 0 <= n3 <= self._element_counts[2]
+ return self._nx[n3][n2][n1]
+
+ def set_node_parameters(self, n1, n2, n3, parameters, nid=None, node_layout=None):
+ assert 0 <= n1 <= self._element_counts[0]
+ assert 0 <= n2 <= self._element_counts[1]
+ assert 0 <= n3 <= self._element_counts[2]
+ assert self._nx[n3][n2][n1] is not None
+ assert self._nids[n3][n2][n1] is None
+ self._nx[n3][n2][n1] = copy.deepcopy(parameters)
+ self._nids[n3][n2][n1] = nid
+ self._prescribed_node_layouts.append((n1, n2, n3, node_layout))
+
def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start_element_identifier=1):
"""
After build() has been called, generate nodes and elements of ellipsoid.
@@ -718,6 +700,8 @@ def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start
for n3 in range(self._element_counts[2] + 1):
for n2 in range(self._element_counts[1] + 1):
for n1 in range(self._element_counts[0] + 1):
+ if self._nids[n3][n2][n1] is not None:
+ continue # prescribed node
parameters = self._nx[n3][n2][n1]
if not parameters:
continue
@@ -741,8 +725,8 @@ def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start
mesh_dimension = 2 if self._surface_only else 3
mesh = fieldmodule.findMeshByDimension(mesh_dimension)
element_identifier = start_element_identifier
+ # return node_identifier, element_identifier
half_counts = [count // 2 for count in self._element_counts]
- node_layout_manager = HermiteNodeLayoutManager()
octant_mesh_group_lists = None
if self._octant_group_lists:
octant_mesh_group_lists = []
@@ -797,8 +781,8 @@ def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start
element_identifier += 1
last_nids_row = nids_row
# around sides
- node_layout_permuted = node_layout_manager.getNodeLayoutRegularPermuted(d3Defined=False)
- node_layout_triple_points = node_layout_manager.getNodeLayoutTriplePoint2D()
+ node_layout_permuted = self._node_layout_manager.getNodeLayoutRegularPermuted(d3Defined=False)
+ node_layout_triple_points = self._node_layout_manager.getNodeLayoutTriplePoint2D()
index_increments = [[0, 1, 0], [-1, 0, 0], [0, -1, 0], [1, 0, 0]]
increment_number = 0
index_increment = index_increments[0]
@@ -912,7 +896,7 @@ def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start
elementtemplate_special.setElementShapeType(Element.SHAPE_TYPE_CUBE)
box_counts = [half_counts[i] - self._trans_count for i in range(3)]
dbox_counts = [2 * box_counts[i] for i in range(3)]
- nid_to_node_layout = self._get_nid_to_node_layout_map_3d(node_layout_manager)
+ nid_to_node_layout = self._get_nid_to_node_layout_map_3d(self._node_layout_manager)
# bottom transition
last_nids_layer = None
last_nx_layer = None
@@ -1185,617 +1169,3 @@ def generate_mesh(self, fieldmodule, coordinates, start_node_identifier=1, start
last_nx_layer = nx_layer
return node_identifier, element_identifier
-
-class EllipsoidOctantMesh:
- """
- Generates one octant of an ellipsoid, 2-D surface or full 3-D volume.
- 2-D outer surface is merely added from a QuadTriangleMesh.
- 3-D volume requires 3 axis-aligned and outer surfaces to each be set from a QuadTriangleMesh,
- then the interior can be built.
- """
-
- def __init__(self, a, b, c, element_counts, transition_element_count, nway_d_factor=0.6):
- """
- A 2-D or 3-D Octant of an ellipsoid.
- Coordinates nx are indexed in 1, 2, 3 directions from origin at index 0, 0, 0
- with holes around the corners and 3-way points.
- :param a: Axis length (radius) in x direction.
- :param b: Axis length (radius) in y direction.
- :param c: Axis length (radius) in z direction.
- :param element_counts: Number of elements across octant only in 1, 2, 3 axes.
- :param transition_element_count: Number of transition elements around outside >= 1.
- :param nway_d_factor: Value, normally from 0.5 to 1.0 giving n-way derivative magnitude as a proportion
- of the minimum regular magnitude sampled to the n-way point. This reflects that distances from the mid-side
- of a triangle to the centre are shorter, so the derivative in the middle must be smaller.
- """
- assert all((count >= 2) for count in element_counts)
- assert 1 <= transition_element_count <= (min(element_counts) - 1)
- self._a = a
- self._b = b
- self._c = c
- self._element_counts = element_counts
- self._trans_count = transition_element_count
- self._nway_d_factor = nway_d_factor
- self._element_count12 = element_counts[0] + element_counts[1] - 2 * transition_element_count
- self._element_count13 = element_counts[0] + element_counts[2] - 2 * transition_element_count
- self._element_count23 = element_counts[1] + element_counts[2] - 2 * transition_element_count
- # counts of elements to 3-way point opposite to 3 node at axis 1, axis 2, axis 3
- self._box_counts = [self._element_counts[i] - self._trans_count for i in range(3)]
- none_parameters = [None] * 4 # x, d1, d2, d3
- self._nx = [] # shield mesh with holes over n3, n2, n1, d
- for n3 in range(element_counts[2] + 1):
- # index into transition zone
- trans3 = transition_element_count + n3 - element_counts[2]
- nx_layer = []
- for n2 in range(element_counts[1] + 1):
- # index into transition zone
- trans2 = transition_element_count + n2 - element_counts[1]
- nx_row = []
- # s = ""
- for n1 in range(element_counts[0] + 1):
- # index into transition zone
- trans1 = transition_element_count + n1 - element_counts[0]
- if (((trans1 <= 0) and (trans2 <= 0) and (trans3 <= 0)) or
- (trans1 == trans2 == trans3) or
- ((trans1 < 0) and ((trans2 == trans3) or (trans2 < 0))) or
- ((trans2 < 0) and ((trans3 == trans1) or (trans3 < 0))) or
- ((trans3 < 0) and ((trans1 == trans2) or (trans1 < 0)))):
- parameters = copy.copy(none_parameters)
- # s += "[]"
- else:
- parameters = None
- # s += " "
- nx_row.append(parameters)
- nx_layer.append(nx_row)
- # print(s)
- self._nx.append(nx_layer)
-
- def get_parameters(self):
- """
- Get parameters array e.g. for copying to ellipsoid.
- :return: Internal parameters array self._nx. Not to be modified.
- """
- return self._nx
-
- def set_triangle_abc(self, trimesh: QuadTriangleMesh):
- """
- Set parameters on the outer abc surface triangle of octant.
- :param trimesh: Coordinates to set on outer surface.
- """
- assert trimesh.get_element_count12() == self._element_count12
- assert trimesh.get_element_count13() == self._element_count13
- assert trimesh.get_element_count23() == self._element_count23
- start_indexes = [self._element_counts[0], 0, 0]
- for n3 in range(self._box_counts[2]):
- px, pd1, pd2, pd3 = trimesh.get_parameters12(n3)
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[0, 1, 0], [-1, 0, 0]])
- start_indexes[2] += 1
- start_indexes = [0, 0, self._element_counts[2]]
- for n2 in range(self._box_counts[1]):
- px, pd1, pd2, pd3 = trimesh.get_parameters31(n2, self._box_counts[0] + 1)
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
- start_indexes[1] += 1
- start_indexes = [0, self._element_counts[1], self._element_counts[2]]
- px, pd1, pd2, pd3 = trimesh.get_parameters_diagonal()
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
-
- def set_triangle_abo(self, trimesh: QuadTriangleMesh):
- """
- Set parameters on triangle 1-2-origin, an inner surface of octant.
- :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
- """
- assert trimesh.get_element_count12() == self._element_count12
- assert trimesh.get_element_count13() == self._element_counts[0]
- assert trimesh.get_element_count23() == self._element_counts[1]
- start_indexes = [self._element_counts[0], 0, 0]
- for n0 in range(self._trans_count):
- px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, -3, 2]], start_indexes, [[0, 1, 0], [-1, 0, 0]])
- start_indexes[0] -= 1
- start_indexes = [0, 0, 0]
- for n2 in range(self._box_counts[1] + 1):
- px, pd1, pd2, pd3 = (trimesh.get_parameters31(n2, self._box_counts[0] + 1) if (n2 < self._box_counts[1])
- else trimesh.get_parameters_diagonal())
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
- start_indexes[1] += 1
-
- def set_triangle_aco(self, trimesh: QuadTriangleMesh):
- """
- Set parameters on triangle 1-3-origin, an inner surface of octant.
- :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
- """
- assert trimesh.get_element_count12() == self._element_count13
- assert trimesh.get_element_count13() == self._element_counts[0]
- assert trimesh.get_element_count23() == self._element_counts[2]
- start_indexes = [self._element_counts[0], 0, 0]
- for n0 in range(self._trans_count):
- px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
- self._set_coordinates_across(
- [px, pd1, pd2, pd3], [[0, 2, -3, -1], [0, -1, -3, -2]], start_indexes, [[0, 0, 1], [-1, 0, 0]])
- start_indexes[0] -= 1
- start_indexes = [0, 0, 0]
- for n3 in range(self._box_counts[2] + 1):
- px, pd1, pd2, pd3 = (trimesh.get_parameters31(n3, self._box_counts[0] + 1) if (n3 < self._box_counts[2])
- else trimesh.get_parameters_diagonal())
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 3, -2]], start_indexes, [[1, 0, 0]])
- start_indexes[2] += 1
-
- def set_triangle_bco(self, trimesh: QuadTriangleMesh):
- """
- Set parameters on triangle 2-3-origin, an inner surface of octant.
- :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
- """
- assert trimesh.get_element_count12() == self._element_count23
- assert trimesh.get_element_count13() == self._element_counts[1]
- assert trimesh.get_element_count23() == self._element_counts[2]
- start_indexes = [0, self._element_counts[1], 0]
- for n0 in range(self._trans_count):
- px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
- self._set_coordinates_across(
- [px, pd1, pd2, pd3], [[0, 2, -3, -1], [0, -2, -3, 1]], start_indexes, [[0, 0, 1], [0, -1, 0]])
- start_indexes[1] -= 1
- start_indexes = [0, 0, 0]
- for n3 in range(self._box_counts[2] + 1):
- px, pd1, pd2, pd3 = (trimesh.get_parameters31(n3, self._box_counts[1] + 1) if (n3 < self._box_counts[2])
- else trimesh.get_parameters_diagonal())
- self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 2, 3, 1]], start_indexes, [[0, 1, 0]])
- start_indexes[2] += 1
-
- def _get_transitions(self, indexes):
- """
- For each index direction, get False if in core or True if in transition zone.
- :param indexes: Location indexes in 1, 2, 3 directions.
- :return: Transition 1, 2, 3 directions.
- """
- return [(indexes[i] - self._box_counts[i]) > 0 for i in range(3)]
-
- def _set_coordinates_across(self, parameters, parameter_indexes, start_indexes, index_increments,
- skip_start=False, skip_end=False, blend=False):
- """
- Insert parameters across the coordinates array.
- :param parameters: List of lists of N node parameters e.g. [px, pd1, pd2, pd3]
- :param parameter_indexes: Lists of parameter indexes where x=0, d1=1, d2=2, d3=3. Starts with first and
- advances at transitions change, then stays on the last. Can be negative to invert vector.
- e.g. [[0, 1, 2], [0, -2, 1]] for [x, d1, d2] then [x, -d2, d1] from first corner.
- :param start_indexes: Indexes into nx array for start point.
- :param index_increments: List of increments in indexes. Starts with first and uses next at each transition
- change, then stays on the last.
- :param skip_start: Set to True to skip the first value.
- :param skip_end: Set to True to skip the last value.
- :param blend: Set to True to blend parameters with any old parameters at locations.
- """
- indexes = start_indexes
- parameter_number = 0
- parameter_index = parameter_indexes[0]
- increment_number = 0
- index_increment = index_increments[0]
- start_n = 1 if skip_start else 0
- last_n = len(parameters[0]) - 1
- limit_n = len(parameters[0]) - (1 if skip_end else 0)
- last_trans = self._get_transitions(indexes)
- for n in range(start_n, limit_n):
- if n > 0:
- while True:
- indexes = [indexes[c] + index_increment[c] for c in range(3)]
- # skip over blank transition coordinates
- if self._nx[indexes[2]][indexes[1]][indexes[0]]:
- break
- trans = self._get_transitions(indexes)
- if last_trans and (trans != last_trans):
- if parameter_number < (len(parameter_indexes) - 1):
- parameter_number += 1
- parameter_index = parameter_indexes[parameter_number]
- if increment_number < (len(index_increments) - 1):
- increment_number += 1
- index_increment = index_increments[increment_number]
- nx = self._nx[indexes[2]][indexes[1]][indexes[0]]
- for parameter, spix in zip(parameters, parameter_index):
- if not parameter[n]:
- continue
- new_parameter = [-d for d in parameter[n]] if (spix < 0) else copy.copy(parameter[n])
- pix = abs(spix)
- if blend and nx[pix]:
- if pix == 0:
- # for fairness, move to surface before blending
- if any(indexes[i] == self._element_counts[i] for i in range(3)):
- new_parameter = moveCoordinatesToEllipsoidSurface(self._a, self._b, self._c, new_parameter)
- new_parameter = [0.5 * (nx[pix][c] + new_parameter[c]) for c in range(3)]
- else:
- # harmonic mean to cope with significant element size differences on boundary
- new_parameter = linearlyInterpolateVectors(
- nx[pix], new_parameter, 0.5, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)
- nx[pix] = new_parameter
- last_trans = trans
-
- def _smooth_derivative_across(self, start_indexes, end_indexes, index_increments, derivative_indexes,
- fix_start_direction=True, fix_end_direction=True, move_d_to_surface=False):
- """
- Smooth derivatives across octant.
- :param start_indexes: Indexes of first point.
- :param end_indexes: Indexes of last point.
- :param index_increments: List of increments in indexes. Starts with first and advances at transitions change,
- then stays on the last.
- :param derivative_indexes: List of signed derivative parameter index to set along where 1=d1, 2=d2, 3=d3.
- Starts with first and advances at transitions change, then stays on the last. Can be negative to invert vector.
- e.g. [1, -2] for d1 then -d2 from first transition change.
- :param fix_start_direction: Set to True to keep the start direction but scale its magnitude.
- :param fix_end_direction: Set to True to keep the end direction but scale its magnitude.
- :param move_d_to_surface: Set to True to force derivatives to surface tangents at their current location.
- Only use for smoothing over the outer surface of ellipsoid.
- """
- indexes = start_indexes
- derivative_number = 0
- derivative_index = derivative_indexes[0]
- increment_number = 0
- index_increment = index_increments[0]
- indexes_list = []
- derivative_index_list = []
- px = []
- pd = []
- n = 0
- last_trans = self._get_transitions(indexes)
- while True:
- if n > 0:
- if indexes == end_indexes:
- break
- while True:
- indexes = [indexes[i] + index_increment[i] for i in range(3)]
- # skip over blank coordinates in transition zone
- if self._nx[indexes[2]][indexes[1]][indexes[0]]:
- break
- trans = self._get_transitions(indexes)
- if last_trans and (trans != last_trans):
- if derivative_number < (len(derivative_indexes) - 1):
- derivative_number += 1
- derivative_index = derivative_indexes[derivative_number]
- if increment_number < (len(index_increments) - 1):
- increment_number += 1
- index_increment = index_increments[increment_number]
- parameters = self._nx[indexes[2]][indexes[1]][indexes[0]]
- x = parameters[0]
- indexes_list.append(copy.copy(indexes))
- spix = derivative_index
- derivative_index_list.append(spix)
- pix = abs(spix)
- if parameters[pix]:
- d = [-ad for ad in parameters[pix]] if (spix < 0) else parameters[pix]
- else:
- d = [0.0, 0.0, 0.0]
- px.append(x)
- pd.append(d)
- n += 1
- last_trans = trans
- sd = smoothCubicHermiteDerivativesLine(
- px, pd, fixStartDirection=fix_start_direction, fixEndDirection=fix_end_direction)
- if move_d_to_surface:
- for n in range(1, len(sd) - 1):
- sd[n] = self._move_d_to_surface(px[n], sd[n])
- sd = smoothCubicHermiteDerivativesLine(px, sd, fixAllDirections=True)
- for n in range(len(sd)):
- indexes = indexes_list[n]
- spix = derivative_index_list[n]
- new_derivative = [-d for d in sd[n]] if (spix < 0) else sd[n]
- pix = abs(spix)
- self._nx[indexes[2]][indexes[1]][indexes[0]][pix] = new_derivative
-
- def build_interior(self):
- """
- Determine interior coordinates from surface coordinates.
- """
- # determine 4-way point location from mean curves between side points linking to it
- point12 = self._nx[0][self._box_counts[1]][self._box_counts[0]]
- point13 = self._nx[self._box_counts[2]][0][self._box_counts[0]]
- point23 = self._nx[self._box_counts[2]][self._box_counts[1]][0]
- point123 = self._nx[self._element_counts[2]][self._element_counts[1]][self._element_counts[0]]
-
- x_4way, d_4way = get_nway_point(
- [point23[0], point13[0], point12[0], point123[0]],
- [point23[1], point13[2], point12[3], [-d for d in point123[3]]],
- [self._box_counts[0], self._box_counts[1], self._box_counts[2], self._trans_count],
- sampleHermiteCurve, nway_d_factor=self._nway_d_factor)
-
- # smooth sample from sides to 3-way points using end derivatives
- min_weight = 1 # GRC revisit, remove?
- ax, ad1 = sampleHermiteCurve(
- point23[0], point23[1], None, x_4way, d_4way[0], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- bx, bd2 = sampleHermiteCurve(
- point13[0], point13[2], None, x_4way, d_4way[1], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- cx, cd3 = sampleHermiteCurve(
- point12[0], point12[3], None, x_4way, d_4way[2], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- tx, td3 = sampleHermiteCurve(
- point123[0], [-d for d in point123[3]], None, x_4way, d_4way[3], None, self._trans_count,
- start_weight=self._trans_count + min_weight, end_weight=1.0 + min_weight, end_transition=True)
-
- self._set_coordinates_across([ax, ad1], [[0, 1]], [0, self._box_counts[1], self._box_counts[2]], [[1, 0, 0]])
- self._set_coordinates_across([bx, bd2], [[0, 2]], [self._box_counts[0], 0, self._box_counts[2]], [[0, 1, 0]])
- self._set_coordinates_across([cx, cd3], [[0, 3]], [self._box_counts[0], self._box_counts[1], 0], [[0, 0, 1]])
- self._set_coordinates_across([tx, td3], [[0, -3]],
- [self._element_counts[0], self._element_counts[1], self._element_counts[2]],
- [[-1, -1, -1]], skip_end=True)
-
- # sample up to 3-way lines connecting to 4-way point
- for n3 in range(1, self._box_counts[2]):
- point13 = self._nx[n3][0][self._box_counts[0]]
- point23 = self._nx[n3][self._box_counts[1]][0]
- point123 = self._nx[n3][self._element_counts[1]][self._element_counts[0]]
- point_3way = self._nx[n3][self._box_counts[1]][self._box_counts[0]]
-
- x_3way, d_3way = get_nway_point(
- [point23[0], point13[0], point123[0]],
- [point23[1], point13[2], [-d for d in point123[3]]],
- [self._box_counts[0], self._box_counts[1], self._trans_count],
- sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
-
- ax, ad1, ad2 = sampleHermiteCurve(
- point23[0], point23[1], point23[2], x_3way, d_3way[0], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- ad1[-1] = d_3way[0]
- ad2[-1] = d_3way[1]
- bx, bd2, bd1 = sampleHermiteCurve(
- point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- bd1[-1] = d_3way[0]
- bd2[-1] = d_3way[1]
- tx, td3, td1 = sampleHermiteCurve(
- point123[0], [-d for d in point123[3]], point123[1], x_3way, d_3way[2], None,
- self._trans_count, start_weight=self._trans_count + min_weight,
- end_weight=1.0 + min_weight, end_transition=True)
-
- self._set_coordinates_across([ax, ad1, ad2], [[0, 1, 2]], [0, self._box_counts[1], n3], [[1, 0, 0]])
- self._set_coordinates_across([bx, bd1, bd2], [[0, 1, 2]], [self._box_counts[0], 0, n3], [[0, 1, 0]])
- self._set_coordinates_across(
- [tx, td1, td3], [[0, 1, -3]], [self._element_counts[0], self._element_counts[1], n3], [[-1, -1, 0]],
- skip_end=True)
-
- for n2 in range(1, self._box_counts[1]):
- point12 = self._nx[0][n2][self._box_counts[0]]
- point23 = self._nx[self._box_counts[2]][n2][0]
- point123 = self._nx[self._element_counts[2]][n2][self._element_counts[0]]
- point_3way = self._nx[self._box_counts[2]][n2][self._box_counts[0]]
-
- x_3way, d_3way = get_nway_point(
- [point23[0], point12[0], point123[0]],
- [point23[1], point12[3], [-d for d in point123[3]]],
- [self._box_counts[0], self._box_counts[2], self._trans_count],
- sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
-
- ax, ad1, ad3 = sampleHermiteCurve(
- point23[0], point23[1], point23[3], x_3way, d_3way[0], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- ad1[-1] = d_3way[0]
- ad3[-1] = d_3way[1]
- bx, bd3, bd1 = sampleHermiteCurve(
- point12[0], point12[3], point12[1], x_3way, d_3way[1], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- bd1[-1] = d_3way[0]
- bd3[-1] = d_3way[1]
- tx, td3, td1 = sampleHermiteCurve(
- point123[0], [-d for d in point123[3]], point123[1], x_3way, d_3way[2], None,
- self._trans_count, start_weight=self._trans_count + min_weight,
- end_weight=1.0 + min_weight, end_transition=True)
-
- self._set_coordinates_across(
- [ax, ad1, ad3], [[0, 1, 3]], [0, n2, self._box_counts[2]], [[1, 0, 0]], blend=True)
- self._set_coordinates_across(
- [bx, bd1, bd3], [[0, 1, 3]], [self._box_counts[0], n2, 0], [[0, 0, 1]], blend=True)
- self._set_coordinates_across(
- [tx, td1, td3], [[0, 1, -3]], [self._element_counts[0], n2, self._element_counts[2]], [[-1, 0, -1]],
- skip_end=True, blend=True)
-
- for n1 in range(1, self._box_counts[0]):
- point12 = self._nx[0][self._box_counts[1]][n1]
- point13 = self._nx[self._box_counts[2]][0][n1]
- point123 = self._nx[self._element_counts[2]][self._element_counts[1]][n1]
- point_3way = self._nx[self._box_counts[2]][self._box_counts[1]][n1]
-
- x_3way, d_3way = get_nway_point(
- [point13[0], point12[0], point123[0]],
- [point13[2], point12[3], [-d for d in point123[3]]],
- [self._box_counts[0], self._box_counts[2], self._trans_count],
- sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
-
- ax, ad2, ad3 = sampleHermiteCurve(
- point13[0], point13[2], point13[3], x_3way, d_3way[0], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- ad2[-1] = d_3way[0]
- ad3[-1] = d_3way[1]
- bx, bd3, bd2 = sampleHermiteCurve(
- point12[0], point12[3], point12[2], x_3way, d_3way[1], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- bd2[-1] = d_3way[0]
- bd3[-1] = d_3way[1]
- tx, td3, td2 = sampleHermiteCurve(
- point123[0], [-d for d in point123[3]], point123[2], x_3way, d_3way[2], None,
- self._trans_count, start_weight=self._trans_count + min_weight,
- end_weight=1.0 + min_weight, end_transition=True)
-
- self._set_coordinates_across(
- [ax, ad2, ad3], [[0, 2, 3]], [n1, 0, self._box_counts[2]], [[0, 1, 0]], blend=True)
- self._set_coordinates_across(
- [bx, bd2, bd3], [[0, 2, 3]], [n1, self._box_counts[1], 0], [[0, 0, 1]], blend=True)
- self._set_coordinates_across(
- [tx, td2, td3], [[0, 2, -3]], [n1, self._element_counts[1], self._element_counts[2]],
- [[0, -1, -1]], skip_end=True, blend=True)
-
- for nt in range(1, self._trans_count):
- point12 = self._nx[0][self._element_counts[1] - nt][self._element_counts[0] - nt]
- point13 = self._nx[self._element_counts[2] - nt][0][self._element_counts[0] - nt]
- point23 = self._nx[self._element_counts[2] - nt][self._element_counts[1] - nt][0]
- point_3way = \
- self._nx[self._element_counts[2] - nt][self._element_counts[1] - nt][self._element_counts[0] - nt]
-
- x_3way, d_3way = get_nway_point(
- [point23[0], point13[0], point12[0]],
- [point23[1], point13[2], point12[2]],
- [self._box_counts[0], self._box_counts[1], self._box_counts[2]],
- sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
-
- ax, ad1, ad2 = sampleHermiteCurve(
- point23[0], point23[1], point23[2], x_3way, d_3way[0], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- ad1[-1] = d_3way[0]
- ad2[-1] = d_3way[1]
- bx, bd2, bd1 = sampleHermiteCurve(
- point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
- bd1[-1] = d_3way[0]
- bd2[-1] = d_3way[1]
- cx, cd2, cd1 = sampleHermiteCurve(
- point12[0], point12[2], point12[1], x_3way, d_3way[2], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight,
- end_transition=True)
-
- self._set_coordinates_across(
- [ax, ad1, ad2], [[0, 1, 2]], [0, self._element_counts[1] - nt, self._element_counts[2] - nt],
- [[1, 0, 0]], blend=True)
- self._set_coordinates_across(
- [bx, bd1, bd2], [[0, 1, 2]], [self._element_counts[0] - nt, 0, self._element_counts[2] - nt],
- [[0, 1, 0]], blend=True)
- self._set_coordinates_across(
- [cx, cd1, cd2], [[0, 1, 2]], [self._element_counts[0] - nt, self._element_counts[1] - nt, 0],
- [[0, 0, 1]], skip_end=True, blend=True)
-
- # average point coordinates across 3 directions between side faces and surfaces to 4 3-way lines.
- min_weight = 1 # GRC revisit, remove?
- # 1-direction
- for n2 in range(1, self._box_counts[1]):
- for n3 in range(1, self._box_counts[2]):
- start_indexes = [0, n2, n3]
- corner_indexes = [self._box_counts[0], n2, n3]
- end_indexes = [self._element_counts[0], n2, n3]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[1], None, corner[0], corner[1], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], corner[1], None, end[0], end[3], None, self._trans_count,
- start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
- for nt in range(1, self._trans_count):
- start_indexes = [0, n2, self._box_counts[2] + nt]
- corner_indexes = [self._box_counts[0] + nt, n2, self._box_counts[2] + nt]
- end_indexes = [self._box_counts[0] + nt, n2, 0]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[1], None, corner[0], corner[1], None, self._box_counts[0],
- start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], corner[1], None, end[0], [-d for d in end[2]], None, self._box_counts[2],
- start_weight=1.0 + min_weight, end_weight=self._box_counts[2] + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[0, 0, -1]], skip_start=True, skip_end=True, blend=True)
- # 2-direction
- for n3 in range(1, self._box_counts[2]):
- for n1 in range(1, self._box_counts[0]):
- start_indexes = [n1, 0, n3]
- corner_indexes = [n1, self._box_counts[1], n3]
- end_indexes = [n1, self._element_counts[1], n3]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[2], None, corner[0], corner[2], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], corner[2], None, end[0], end[3], None, self._trans_count,
- start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
- for nt in range(1, self._trans_count):
- start_indexes = [self._box_counts[0] + nt, 0, n3]
- corner_indexes = [self._box_counts[0] + nt, self._box_counts[1] + nt, n3]
- end_indexes = [0, self._box_counts[1] + nt, n3]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[1], None, corner[0], corner[1], None, self._box_counts[1],
- start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], corner[1], None, end[0], end[1], None, self._box_counts[0],
- start_weight=1.0 + min_weight, end_weight=self._box_counts[0] + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[-1, 0, 0]], skip_start=True, skip_end=True, blend=True)
- # 3-direction
- for n1 in range(1, self._box_counts[0]):
- for n2 in range(1, self._box_counts[1]):
- start_indexes = [n1, n2, 0]
- corner_indexes = [n1, n2, self._box_counts[2]]
- end_indexes = [n1, n2, self._element_counts[2]]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[3], None, corner[0], corner[3], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], corner[3], None, end[0], end[3], None, self._trans_count,
- start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
- for nt in range(1, self._trans_count):
- start_indexes = [n1, self._box_counts[1] + nt, 0]
- corner_indexes = [n1, self._box_counts[1] + nt, self._box_counts[2] + nt]
- end_indexes = [n1, 0, self._box_counts[2] + nt]
- start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
- corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
- end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
- px, _ = sampleHermiteCurve(
- start[0], start[2], None, corner[0], [-d for d in corner[2]], None, self._box_counts[2],
- start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight)
- self._set_coordinates_across(
- [px], [[0]], start_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
- px, _ = sampleHermiteCurve(
- corner[0], [-d for d in corner[2]], None, end[0], [-d for d in end[2]], None, self._box_counts[1],
- start_weight=1.0 + min_weight, end_weight=self._box_counts[1] + min_weight)
- self._set_coordinates_across(
- [px], [[0]], corner_indexes, [[0, -1, 0]], skip_start=True, skip_end=True, blend=True)
-
- # smooth 1-direction
- for n2 in range(1, self._box_counts[1]):
- for n3 in range(1, self._box_counts[2]):
- self._smooth_derivative_across(
- [0, n2, n3], [self._element_counts[0], n2, n3],
- [[1, 0, 0]], [1, 3], fix_start_direction=True, fix_end_direction=True)
- for nt in range(1, self._trans_count):
- self._smooth_derivative_across(
- [0, n2, self._box_counts[2] + nt], [self._box_counts[0] + nt, n2, 0],
- [[1, 0, 0], [0, 0, -1]], [1, 1, -2], fix_start_direction=True, fix_end_direction=True)
- # smooth 2-direction
- for n3 in range(1, self._box_counts[2]):
- for n1 in range(1, self._box_counts[0]):
- self._smooth_derivative_across(
- [n1, 0, n3], [n1, self._element_counts[1], n3],
- [[0, 1, 0]], [2, 3], fix_start_direction=True, fix_end_direction=True)
- for nt in range(1, self._trans_count):
- self._smooth_derivative_across(
- [self._box_counts[0] + nt, 0, n3], [0, self._box_counts[1] + nt, n3],
- [[0, 1, 0], [-1, 0, 0]], [1], fix_start_direction=True, fix_end_direction=True)
- # smooth 3-direction
- for n1 in range(1, self._box_counts[0]):
- for n2 in range(1, self._box_counts[1]):
- self._smooth_derivative_across(
- [n1, n2, 0], [n1, n2, self._element_counts[2]],
- [[0, 0, 1]], [3], fix_start_direction=True, fix_end_direction=True)
- for nt in range(1, self._trans_count):
- self._smooth_derivative_across(
- [n1, self._box_counts[1] + nt, 0], [n1, 0, self._box_counts[2] + nt],
- [[0, 0, 1], [0, -1, 0]], [2, -2], fix_start_direction=True, fix_end_direction=True)
diff --git a/src/scaffoldmaker/utils/geometry.py b/src/scaffoldmaker/utils/geometry.py
index 7d4025eb..0b522f26 100644
--- a/src/scaffoldmaker/utils/geometry.py
+++ b/src/scaffoldmaker/utils/geometry.py
@@ -7,7 +7,8 @@
import copy
import math
-from cmlibs.maths.vectorops import add, distance, magnitude, mult, normalize, cross, set_magnitude, rejection
+from cmlibs.maths.vectorops import (
+ cross, distance, dot, magnitude, mult, normalize, projection, set_magnitude, rejection)
from scaffoldmaker.utils.interpolation import (
computeCubicHermiteDerivativeScaling, computeHermiteLagrangeDerivativeScaling, getCubicHermiteArcLength,
interpolateHermiteLagrangeDerivative, linearlyInterpolateVectors, sampleCubicHermiteCurves,
@@ -150,31 +151,57 @@ def getEllipseRadiansToX(ax, bx, dx, initialTheta):
return theta
-def getEllipsePointAtTrueAngle(a, b, angle_radians):
+def getEllipsePointAtTrueAngle(a, b, angle_radians, origin=[0.0, 0.0]):
"""
- Get coordinates of intersection point of ellipse centred at origin with line radiating from origin at an angle.
+ Get coordinates of intersection point of ellipse centred at [0.0, 0.0] with line radiating from origin at an angle.
:param a: x/major axis length.
:param b: y/minor axis length.
:param angle_radians: Angle in radians starting at x axis, increasing towards y axis.
- :return: [x, y]
+ :param origin: Origin angles radiate from. Must be inside ellipsoid. Default is ellipsoid centre [0.0, 0.0].
+ :return: [x, y] on ellipse.
"""
# ellipse equation: x ** 2 / a ** 2 + y ** 2 / b ** 2 - 1 = 0
- cos_angle = math.cos(angle_radians)
- sin_angle = math.sin(angle_radians)
+ aa = a * a
+ bb = b * b
+ assert (origin[0] * origin[0] / aa + (origin[1] * origin[1]) / bb) < 1.0
# normal to line direction:
- ni = sin_angle
- nj = -cos_angle
- # line equation: ni * x + nj * y = 0
+ # line equation: ni * x + nj * y - k = 0
+ ni = math.sin(angle_radians)
+ nj = -math.cos(angle_radians)
+ k = dot([ni, nj], origin)
+ ii = ni * ni
+ jj = nj * nj
+ kk = k * k
if math.fabs(nj) > math.fabs(ni):
- # substitute y and solve for x
- denominator = 1.0 / (a * a) + (ni * ni) / (nj * nj * b * b)
- x = math.copysign(math.sqrt(1.0 / denominator), cos_angle)
- y = (-ni / nj) * x
+ # substitute y and solve for x with quadratic equation
+ jj_bb = jj * bb
+ qa = 1.0 / aa + ii / jj_bb
+ qb = -2.0 * k * ni / jj_bb
+ qc = kk / jj_bb - 1.0
+ det = qb * qb - 4.0 * qa * qc
+ sqrt_det = math.sqrt(det)
+ x1 = (-qb + sqrt_det) / (2.0 * qa)
+ x2 = (-qb - sqrt_det) / (2.0 * qa)
+ if x1 * nj < 0.0:
+ x = x1
+ else:
+ x = x2
+ y = (k - ni * x) / nj
else:
- # substitute y and solve for x
- denominator = 1.0 / (b * b) + (nj * nj) / (ni * ni * a * a)
- y = math.copysign(math.sqrt(1.0 / denominator), sin_angle)
- x = (-nj / ni) * y
+ # substitute x and solve for y with quadratic equation
+ ii_aa = ii * aa
+ qa = 1.0 / bb + jj / ii_aa
+ qb = -2.0 * k * nj / ii_aa
+ qc = kk / ii_aa - 1.0
+ det = qb * qb - 4.0 * qa * qc
+ sqrt_det = math.sqrt(det)
+ y1 = (-qb + sqrt_det) / (2.0 * qa)
+ y2 = (-qb - sqrt_det) / (2.0 * qa)
+ if y1 * ni > 0.0:
+ y = y1
+ else:
+ y = y2
+ x = (k - nj * y) / ni
return [x, y]
@@ -492,6 +519,22 @@ def moveDerivativeToEllipsoidSurface(a, b, c, x, start_d):
return set_magnitude(rejection(start_d, n), magnitude(start_d))
+def moveDerivativeToEllipsoidSurfaceInPlane(a, b, c, x, pn, start_d):
+ """
+ Convert derivative at point on surface of ellipsoid to be tangential to it and normal to pn.
+ :param a: x-axis length.
+ :param b: y-axis length.
+ :param c: z-axis length.
+ :param x: Coordinates on surface of ellipsoid.
+ :param pn: Direction normal to plane to force result to be in plane.
+ :param start_d: Derivative near tangential to ellipsoid surface at x.
+ :return: Derivative made tangential to surface with same magnitude.
+ """
+ en = [2.0 * x[0] / (a * a), 2.0 * x[1] / (b * b), 2.0 * x[2] / (c * c)]
+ cp = cross(pn, en)
+ return set_magnitude(projection(start_d, cp), magnitude(start_d))
+
+
def sampleCurveOnEllipsoid(a, b, c, start_x, start_d1, start_d2, end_x, end_d1, end_d2, elements_count,
start_weight=None, end_weight=None, overweighting=1.0, end_transition=False):
"""
@@ -514,7 +557,7 @@ def sampleCurveOnEllipsoid(a, b, c, start_x, start_d1, start_d2, end_x, end_d1,
by distance from the other end.
:param overweighting: Multiplier of arc length to use with initial curve to exaggerate end derivatives.
:param end_transition: If supplied with end_d1, modify size of last element to fit end_d1.
- :return: x[], d1[], d2[]
+ :return: x[], d1[]{, d2[] if start_d2 supplied}
"""
assert (not end_transition) or end_d1
end_d1_mag = magnitude(end_d1) if end_d1 else None
diff --git a/src/scaffoldmaker/utils/hextetrahedronmesh.py b/src/scaffoldmaker/utils/hextetrahedronmesh.py
new file mode 100644
index 00000000..dc0be76f
--- /dev/null
+++ b/src/scaffoldmaker/utils/hextetrahedronmesh.py
@@ -0,0 +1,646 @@
+"""
+Utilities for building 3-D tetrahedron-shaped meshes out of hexahedral elements with cubic Hermite serendipity
+interpolation.
+"""
+from scaffoldmaker.utils.interpolation import (
+ DerivativeScalingMode, get_nway_point, linearlyInterpolateVectors, sampleHermiteCurve,
+ smoothCubicHermiteDerivativesLine)
+from scaffoldmaker.utils.quadtrianglemesh import QuadTriangleMesh
+import copy
+
+
+class HexTetrahedronMesh:
+ """
+ Generates a tetrahedron mesh from c1-continuous hex (cube) elements, with a 3-way or 4-way points inside.
+ Tetrahedron is defined from 4 corner points 0-3, origin and other axis ends as right-handed axes.
+ Parameters on faces are set from QuadTriangleMesh objects.
+ 3-D parameters are interpolated from face parameters.
+ Some limitations on currently supported number of elements in directions are asserted in the constructor.
+ """
+
+ def __init__(self, axis_counts, diag_counts, nway_d_factor=0.6):
+ """
+ :param axis_counts: Number of elements along 0-1, 0-2, 0-3 axes.
+ :param diag_counts: Number of elements along 1-2, 1-3, 2-3 diagonals.
+ Coordinates nx are indexed in 1, 2, 3 directions from origin at index 0, 0, 0
+ with holes around the corners and 3-way or 4-way points.
+ :param nway_d_factor: Value, normally from 0.5 to 1.0 giving n-way derivative magnitude as a proportion
+ of the minimum regular magnitude sampled to the n-way point. This reflects that distances from the mid-side
+ of a triangle to the centre are shorter, so the derivative in the middle must be smaller.
+ """
+ assert all((count >= 2) for count in axis_counts)
+ assert all((count >= 2) for count in diag_counts)
+ # check the faces have valid element counts around them
+ max_diag_count0 = axis_counts[0] + axis_counts[1] - 2
+ assert any((diag_counts[0] == diag_count) for diag_count in range(max_diag_count0, 1, -2))
+ max_diag_count1 = axis_counts[0] + axis_counts[2] - 2
+ assert any((diag_counts[1] == diag_count) for diag_count in range(max_diag_count1, 1, -2))
+ max_diag_count2 = axis_counts[1] + axis_counts[2] - 2
+ assert any((diag_counts[2] == diag_count) for diag_count in range(max_diag_count2, 1, -2))
+ max_diag_count3 = diag_counts[0] + diag_counts[1] - 2
+ assert any((diag_counts[2] == diag_count) for diag_count in range(max_diag_count3, 1, -2))
+ self._axis_counts = copy.copy(axis_counts)
+ self._diag_counts = copy.copy(diag_counts)
+ self._nway_d_factor = nway_d_factor
+
+ # current limitation:
+ # only supports 4-way point for now, consistent with constant transition count away from origin
+ trans_count0 = (axis_counts[0] + axis_counts[1] - diag_counts[0]) // 2
+ trans_count1 = (axis_counts[0] + axis_counts[2] - diag_counts[1]) // 2
+ trans_count2 = (axis_counts[1] + axis_counts[2] - diag_counts[2]) // 2
+ assert trans_count0 == trans_count1 == trans_count2
+ self._trans_count = trans_count0
+
+ # counts of elements to 3-way point opposite to 3 node at axis 1, axis 2, axis 3
+ self._box_counts = [self._axis_counts[i] - self._trans_count for i in range(3)]
+ none_parameters = [None] * 4 # x, d1, d2, d3
+ self._nx = [] # shield mesh with holes over n3, n2, n1, d
+ for n3 in range(axis_counts[2] + 1):
+ # index into transition zone
+ trans3 = self._trans_count + n3 - axis_counts[2]
+ nx_layer = []
+ for n2 in range(axis_counts[1] + 1):
+ # index into transition zone
+ trans2 = self._trans_count + n2 - axis_counts[1]
+ nx_row = []
+ # s = ""
+ for n1 in range(axis_counts[0] + 1):
+ # index into transition zone
+ trans1 = self._trans_count + n1 - axis_counts[0]
+ if (((trans1 <= 0) and (trans2 <= 0) and (trans3 <= 0)) or
+ (trans1 == trans2 == trans3) or
+ ((trans1 < 0) and ((trans2 == trans3) or (trans2 < 0))) or
+ ((trans2 < 0) and ((trans3 == trans1) or (trans3 < 0))) or
+ ((trans3 < 0) and ((trans1 == trans2) or (trans1 < 0)))):
+ parameters = copy.copy(none_parameters)
+ # s += "[]"
+ else:
+ parameters = None
+ # s += " "
+ nx_row.append(parameters)
+ nx_layer.append(nx_row)
+ # print(s)
+ self._nx.append(nx_layer)
+
+ def get_axis_counts(self):
+ return self._axis_counts
+
+ def get_box_counts(self):
+ return self._box_counts
+
+ def get_parameters(self):
+ """
+ Get parameters array e.g. for copying to ellipsoid.
+ :return: Internal parameters array self._nx. Not to be modified.
+ """
+ return self._nx
+
+ def set_triangle_abc(self, trimesh: QuadTriangleMesh):
+ """
+ Set parameters on the outer abc surface triangle of octant.
+ :param trimesh: Coordinates to set on outer surface.
+ """
+ assert trimesh.get_element_count12() == self._diag_counts[0]
+ assert trimesh.get_element_count13() == self._diag_counts[1]
+ assert trimesh.get_element_count23() == self._diag_counts[2]
+ start_indexes = [self._axis_counts[0], 0, 0]
+ for n3 in range(self._box_counts[2]):
+ px, pd1, pd2, pd3 = trimesh.get_parameters12(n3)
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[0, 1, 0], [-1, 0, 0]])
+ start_indexes[2] += 1
+ start_indexes = [0, 0, self._axis_counts[2]]
+ for n2 in range(self._box_counts[1]):
+ px, pd1, pd2, pd3 = trimesh.get_parameters31(n2, self._box_counts[0] + 1)
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
+ start_indexes[1] += 1
+ start_indexes = [0, self._axis_counts[1], self._axis_counts[2]]
+ px, pd1, pd2, pd3 = trimesh.get_parameters_diagonal()
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
+
+ def set_triangle_abo(self, trimesh: QuadTriangleMesh):
+ """
+ Set parameters on triangle 1-2-origin, an inner surface of octant.
+ :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
+ """
+ assert trimesh.get_element_count12() == self._diag_counts[0]
+ assert trimesh.get_element_count13() == self._axis_counts[0]
+ assert trimesh.get_element_count23() == self._axis_counts[1]
+ start_indexes = [self._axis_counts[0], 0, 0]
+ for n0 in range(self._trans_count):
+ px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, -3, 2]], start_indexes, [[0, 1, 0], [-1, 0, 0]])
+ start_indexes[0] -= 1
+ start_indexes = [0, 0, 0]
+ for n2 in range(self._box_counts[1] + 1):
+ px, pd1, pd2, pd3 = (trimesh.get_parameters31(n2, self._box_counts[0] + 1) if (n2 < self._box_counts[1])
+ else trimesh.get_parameters_diagonal())
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 2, 3]], start_indexes, [[1, 0, 0]])
+ start_indexes[1] += 1
+
+ def set_triangle_aco(self, trimesh: QuadTriangleMesh):
+ """
+ Set parameters on triangle 1-3-origin, an inner surface of octant.
+ :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
+ """
+ assert trimesh.get_element_count12() == self._diag_counts[1]
+ assert trimesh.get_element_count13() == self._axis_counts[0]
+ assert trimesh.get_element_count23() == self._axis_counts[2]
+ start_indexes = [self._axis_counts[0], 0, 0]
+ for n0 in range(self._trans_count):
+ px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
+ self._set_coordinates_across(
+ [px, pd1, pd2, pd3], [[0, 2, -3, -1], [0, -1, -3, -2]], start_indexes, [[0, 0, 1], [-1, 0, 0]])
+ start_indexes[0] -= 1
+ start_indexes = [0, 0, 0]
+ for n3 in range(self._box_counts[2] + 1):
+ px, pd1, pd2, pd3 = (trimesh.get_parameters31(n3, self._box_counts[0] + 1) if (n3 < self._box_counts[2])
+ else trimesh.get_parameters_diagonal())
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 1, 3, -2]], start_indexes, [[1, 0, 0]])
+ start_indexes[2] += 1
+
+ def set_triangle_bco(self, trimesh: QuadTriangleMesh):
+ """
+ Set parameters on triangle 2-3-origin, an inner surface of octant.
+ :param trimesh: Triangle coordinate data with x, d1, d2, optional d3.
+ """
+ assert trimesh.get_element_count12() == self._diag_counts[2]
+ assert trimesh.get_element_count13() == self._axis_counts[1]
+ assert trimesh.get_element_count23() == self._axis_counts[2]
+ start_indexes = [0, self._axis_counts[1], 0]
+ for n0 in range(self._trans_count):
+ px, pd1, pd2, pd3 = trimesh.get_parameters12(n0)
+ self._set_coordinates_across(
+ [px, pd1, pd2, pd3], [[0, 2, -3, -1], [0, -2, -3, 1]], start_indexes, [[0, 0, 1], [0, -1, 0]])
+ start_indexes[1] -= 1
+ start_indexes = [0, 0, 0]
+ for n3 in range(self._box_counts[2] + 1):
+ px, pd1, pd2, pd3 = (trimesh.get_parameters31(n3, self._box_counts[1] + 1) if (n3 < self._box_counts[2])
+ else trimesh.get_parameters_diagonal())
+ self._set_coordinates_across([px, pd1, pd2, pd3], [[0, 2, 3, 1]], start_indexes, [[0, 1, 0]])
+ start_indexes[2] += 1
+
+ def _get_transitions(self, indexes):
+ """
+ For each index direction, get False if in core or True if in transition zone.
+ :param indexes: Location indexes in 1, 2, 3 directions.
+ :return: Transition 1, 2, 3 directions.
+ """
+ return [(indexes[i] - self._box_counts[i]) > 0 for i in range(3)]
+
+ def _set_coordinates_across(self, parameters, parameter_indexes, start_indexes, index_increments,
+ skip_start=False, skip_end=False, blend=False):
+ """
+ Insert parameters across the coordinates array.
+ :param parameters: List of lists of N node parameters e.g. [px, pd1, pd2, pd3]
+ :param parameter_indexes: Lists of parameter indexes where x=0, d1=1, d2=2, d3=3. Starts with first and
+ advances at transitions change, then stays on the last. Can be negative to invert vector.
+ e.g. [[0, 1, 2], [0, -2, 1]] for [x, d1, d2] then [x, -d2, d1] from first corner.
+ :param start_indexes: Indexes into nx array for start point.
+ :param index_increments: List of increments in indexes. Starts with first and uses next at each transition
+ change, then stays on the last.
+ :param skip_start: Set to True to skip the first value.
+ :param skip_end: Set to True to skip the last value.
+ :param blend: Set to True to blend parameters with any old parameters at locations.
+ """
+ indexes = start_indexes
+ parameter_number = 0
+ parameter_index = parameter_indexes[0]
+ increment_number = 0
+ index_increment = index_increments[0]
+ start_n = 1 if skip_start else 0
+ last_n = len(parameters[0]) - 1
+ limit_n = len(parameters[0]) - (1 if skip_end else 0)
+ last_trans = self._get_transitions(indexes)
+ for n in range(start_n, limit_n):
+ if n > 0:
+ while True:
+ indexes = [indexes[c] + index_increment[c] for c in range(3)]
+ # skip over blank transition coordinates
+ if self._nx[indexes[2]][indexes[1]][indexes[0]]:
+ break
+ trans = self._get_transitions(indexes)
+ if last_trans and (trans != last_trans):
+ if parameter_number < (len(parameter_indexes) - 1):
+ parameter_number += 1
+ parameter_index = parameter_indexes[parameter_number]
+ if increment_number < (len(index_increments) - 1):
+ increment_number += 1
+ index_increment = index_increments[increment_number]
+ nx = self._nx[indexes[2]][indexes[1]][indexes[0]]
+ for parameter, spix in zip(parameters, parameter_index):
+ if not parameter[n]:
+ continue
+ new_parameter = [-d for d in parameter[n]] if (spix < 0) else copy.copy(parameter[n])
+ pix = abs(spix)
+ if blend and nx[pix]:
+ if pix == 0:
+ new_parameter = [0.5 * (nx[pix][c] + new_parameter[c]) for c in range(3)]
+ else:
+ # harmonic mean to cope with significant element size differences on boundary
+ new_parameter = linearlyInterpolateVectors(
+ nx[pix], new_parameter, 0.5, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)
+ nx[pix] = new_parameter
+ last_trans = trans
+
+ def _smooth_derivative_across(self, start_indexes, end_indexes, index_increments, derivative_indexes,
+ fix_start_direction=True, fix_end_direction=True):
+ """
+ Smooth derivatives across octant.
+ :param start_indexes: Indexes of first point.
+ :param end_indexes: Indexes of last point.
+ :param index_increments: List of increments in indexes. Starts with first and advances at transitions change,
+ then stays on the last.
+ :param derivative_indexes: List of signed derivative parameter index to set along where 1=d1, 2=d2, 3=d3.
+ Starts with first and advances at transitions change, then stays on the last. Can be negative to invert vector.
+ e.g. [1, -2] for d1 then -d2 from first transition change.
+ :param fix_start_direction: Set to True to keep the start direction but scale its magnitude.
+ :param fix_end_direction: Set to True to keep the end direction but scale its magnitude.
+ """
+ indexes = start_indexes
+ derivative_number = 0
+ derivative_index = derivative_indexes[0]
+ increment_number = 0
+ index_increment = index_increments[0]
+ indexes_list = []
+ derivative_index_list = []
+ px = []
+ pd = []
+ n = 0
+ last_trans = self._get_transitions(indexes)
+ while True:
+ if n > 0:
+ if indexes == end_indexes:
+ break
+ while True:
+ indexes = [indexes[i] + index_increment[i] for i in range(3)]
+ # skip over blank coordinates in transition zone
+ if self._nx[indexes[2]][indexes[1]][indexes[0]]:
+ break
+ trans = self._get_transitions(indexes)
+ if last_trans and (trans != last_trans):
+ if derivative_number < (len(derivative_indexes) - 1):
+ derivative_number += 1
+ derivative_index = derivative_indexes[derivative_number]
+ if increment_number < (len(index_increments) - 1):
+ increment_number += 1
+ index_increment = index_increments[increment_number]
+ parameters = self._nx[indexes[2]][indexes[1]][indexes[0]]
+ x = parameters[0]
+ indexes_list.append(copy.copy(indexes))
+ spix = derivative_index
+ derivative_index_list.append(spix)
+ pix = abs(spix)
+ if parameters[pix]:
+ d = [-ad for ad in parameters[pix]] if (spix < 0) else parameters[pix]
+ else:
+ d = [0.0, 0.0, 0.0]
+ px.append(x)
+ pd.append(d)
+ n += 1
+ last_trans = trans
+ sd = smoothCubicHermiteDerivativesLine(
+ px, pd, fixStartDirection=fix_start_direction, fixEndDirection=fix_end_direction)
+ for n in range(len(sd)):
+ indexes = indexes_list[n]
+ spix = derivative_index_list[n]
+ new_derivative = [-d for d in sd[n]] if (spix < 0) else sd[n]
+ pix = abs(spix)
+ self._nx[indexes[2]][indexes[1]][indexes[0]][pix] = new_derivative
+
+ def build_interior(self):
+ """
+ Determine interior coordinates from surface coordinates.
+ """
+ # determine 4-way point location from mean curves between side points linking to it
+ point12 = self._nx[0][self._box_counts[1]][self._box_counts[0]]
+ point13 = self._nx[self._box_counts[2]][0][self._box_counts[0]]
+ point23 = self._nx[self._box_counts[2]][self._box_counts[1]][0]
+ point123 = self._nx[self._axis_counts[2]][self._axis_counts[1]][self._axis_counts[0]]
+
+ x_4way, d_4way = get_nway_point(
+ [point23[0], point13[0], point12[0], point123[0]],
+ [point23[1], point13[2], point12[3], [-d for d in point123[3]]],
+ [self._box_counts[0], self._box_counts[1], self._box_counts[2], self._trans_count],
+ sampleHermiteCurve, nway_d_factor=self._nway_d_factor)
+
+ # smooth sample from sides to 3-way points using end derivatives
+ min_weight = 1 # GRC revisit, remove?
+ ax, ad1 = sampleHermiteCurve(
+ point23[0], point23[1], None, x_4way, d_4way[0], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ bx, bd2 = sampleHermiteCurve(
+ point13[0], point13[2], None, x_4way, d_4way[1], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ cx, cd3 = sampleHermiteCurve(
+ point12[0], point12[3], None, x_4way, d_4way[2], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ tx, td3 = sampleHermiteCurve(
+ point123[0], [-d for d in point123[3]], None, x_4way, d_4way[3], None, self._trans_count,
+ start_weight=self._trans_count + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+
+ self._set_coordinates_across([ax, ad1], [[0, 1]], [0, self._box_counts[1], self._box_counts[2]], [[1, 0, 0]])
+ self._set_coordinates_across([bx, bd2], [[0, 2]], [self._box_counts[0], 0, self._box_counts[2]], [[0, 1, 0]])
+ self._set_coordinates_across([cx, cd3], [[0, 3]], [self._box_counts[0], self._box_counts[1], 0], [[0, 0, 1]])
+ self._set_coordinates_across([tx, td3], [[0, -3]],
+ [self._axis_counts[0], self._axis_counts[1], self._axis_counts[2]],
+ [[-1, -1, -1]], skip_end=True)
+
+ # sample up to 3-way lines connecting to 4-way point
+ for n3 in range(1, self._box_counts[2]):
+ point13 = self._nx[n3][0][self._box_counts[0]]
+ point23 = self._nx[n3][self._box_counts[1]][0]
+ point123 = self._nx[n3][self._axis_counts[1]][self._axis_counts[0]]
+ point_3way = self._nx[n3][self._box_counts[1]][self._box_counts[0]]
+
+ x_3way, d_3way = get_nway_point(
+ [point23[0], point13[0], point123[0]],
+ [point23[1], point13[2], [-d for d in point123[3]]],
+ [self._box_counts[0], self._box_counts[1], self._trans_count],
+ sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
+
+ ax, ad1, ad2 = sampleHermiteCurve(
+ point23[0], point23[1], point23[2], x_3way, d_3way[0], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ ad1[-1] = d_3way[0]
+ ad2[-1] = d_3way[1]
+ bx, bd2, bd1 = sampleHermiteCurve(
+ point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ bd1[-1] = d_3way[0]
+ bd2[-1] = d_3way[1]
+ tx, td3, td1 = sampleHermiteCurve(
+ point123[0], [-d for d in point123[3]], point123[1], x_3way, d_3way[2], None,
+ self._trans_count, start_weight=self._trans_count + min_weight,
+ end_weight=1.0 + min_weight, end_transition=True)
+
+ self._set_coordinates_across([ax, ad1, ad2], [[0, 1, 2]], [0, self._box_counts[1], n3], [[1, 0, 0]])
+ self._set_coordinates_across([bx, bd1, bd2], [[0, 1, 2]], [self._box_counts[0], 0, n3], [[0, 1, 0]])
+ self._set_coordinates_across(
+ [tx, td1, td3], [[0, 1, -3]], [self._axis_counts[0], self._axis_counts[1], n3], [[-1, -1, 0]],
+ skip_end=True)
+
+ for n2 in range(1, self._box_counts[1]):
+ point12 = self._nx[0][n2][self._box_counts[0]]
+ point23 = self._nx[self._box_counts[2]][n2][0]
+ point123 = self._nx[self._axis_counts[2]][n2][self._axis_counts[0]]
+ point_3way = self._nx[self._box_counts[2]][n2][self._box_counts[0]]
+
+ x_3way, d_3way = get_nway_point(
+ [point23[0], point12[0], point123[0]],
+ [point23[1], point12[3], [-d for d in point123[3]]],
+ [self._box_counts[0], self._box_counts[2], self._trans_count],
+ sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
+
+ ax, ad1, ad3 = sampleHermiteCurve(
+ point23[0], point23[1], point23[3], x_3way, d_3way[0], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ ad1[-1] = d_3way[0]
+ ad3[-1] = d_3way[1]
+ bx, bd3, bd1 = sampleHermiteCurve(
+ point12[0], point12[3], point12[1], x_3way, d_3way[1], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ bd1[-1] = d_3way[0]
+ bd3[-1] = d_3way[1]
+ tx, td3, td1 = sampleHermiteCurve(
+ point123[0], [-d for d in point123[3]], point123[1], x_3way, d_3way[2], None,
+ self._trans_count, start_weight=self._trans_count + min_weight,
+ end_weight=1.0 + min_weight, end_transition=True)
+
+ self._set_coordinates_across(
+ [ax, ad1, ad3], [[0, 1, 3]], [0, n2, self._box_counts[2]], [[1, 0, 0]], blend=True)
+ self._set_coordinates_across(
+ [bx, bd1, bd3], [[0, 1, 3]], [self._box_counts[0], n2, 0], [[0, 0, 1]], blend=True)
+ self._set_coordinates_across(
+ [tx, td1, td3], [[0, 1, -3]], [self._axis_counts[0], n2, self._axis_counts[2]], [[-1, 0, -1]],
+ skip_end=True, blend=True)
+
+ for n1 in range(1, self._box_counts[0]):
+ point12 = self._nx[0][self._box_counts[1]][n1]
+ point13 = self._nx[self._box_counts[2]][0][n1]
+ point123 = self._nx[self._axis_counts[2]][self._axis_counts[1]][n1]
+ point_3way = self._nx[self._box_counts[2]][self._box_counts[1]][n1]
+
+ x_3way, d_3way = get_nway_point(
+ [point13[0], point12[0], point123[0]],
+ [point13[2], point12[3], [-d for d in point123[3]]],
+ [self._box_counts[0], self._box_counts[2], self._trans_count],
+ sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
+
+ ax, ad2, ad3 = sampleHermiteCurve(
+ point13[0], point13[2], point13[3], x_3way, d_3way[0], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ ad2[-1] = d_3way[0]
+ ad3[-1] = d_3way[1]
+ bx, bd3, bd2 = sampleHermiteCurve(
+ point12[0], point12[3], point12[2], x_3way, d_3way[1], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ bd2[-1] = d_3way[0]
+ bd3[-1] = d_3way[1]
+ tx, td3, td2 = sampleHermiteCurve(
+ point123[0], [-d for d in point123[3]], point123[2], x_3way, d_3way[2], None,
+ self._trans_count, start_weight=self._trans_count + min_weight,
+ end_weight=1.0 + min_weight, end_transition=True)
+
+ self._set_coordinates_across(
+ [ax, ad2, ad3], [[0, 2, 3]], [n1, 0, self._box_counts[2]], [[0, 1, 0]], blend=True)
+ self._set_coordinates_across(
+ [bx, bd2, bd3], [[0, 2, 3]], [n1, self._box_counts[1], 0], [[0, 0, 1]], blend=True)
+ self._set_coordinates_across(
+ [tx, td2, td3], [[0, 2, -3]], [n1, self._axis_counts[1], self._axis_counts[2]],
+ [[0, -1, -1]], skip_end=True, blend=True)
+
+ for nt in range(1, self._trans_count):
+ point12 = self._nx[0][self._axis_counts[1] - nt][self._axis_counts[0] - nt]
+ point13 = self._nx[self._axis_counts[2] - nt][0][self._axis_counts[0] - nt]
+ point23 = self._nx[self._axis_counts[2] - nt][self._axis_counts[1] - nt][0]
+ point_3way = \
+ self._nx[self._axis_counts[2] - nt][self._axis_counts[1] - nt][self._axis_counts[0] - nt]
+
+ x_3way, d_3way = get_nway_point(
+ [point23[0], point13[0], point12[0]],
+ [point23[1], point13[2], point12[2]],
+ [self._box_counts[0], self._box_counts[1], self._box_counts[2]],
+ sampleHermiteCurve, prescribed_x_nway=point_3way[0], nway_d_factor=self._nway_d_factor)
+
+ ax, ad1, ad2 = sampleHermiteCurve(
+ point23[0], point23[1], point23[2], x_3way, d_3way[0], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ ad1[-1] = d_3way[0]
+ ad2[-1] = d_3way[1]
+ bx, bd2, bd1 = sampleHermiteCurve(
+ point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ bd1[-1] = d_3way[0]
+ bd2[-1] = d_3way[1]
+ cx, cd2, cd1 = sampleHermiteCurve(
+ point12[0], point12[2], point12[1], x_3way, d_3way[2], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight,
+ end_transition=True)
+
+ self._set_coordinates_across(
+ [ax, ad1, ad2], [[0, 1, 2]], [0, self._axis_counts[1] - nt, self._axis_counts[2] - nt],
+ [[1, 0, 0]], blend=True)
+ self._set_coordinates_across(
+ [bx, bd1, bd2], [[0, 1, 2]], [self._axis_counts[0] - nt, 0, self._axis_counts[2] - nt],
+ [[0, 1, 0]], blend=True)
+ self._set_coordinates_across(
+ [cx, cd1, cd2], [[0, 1, 2]], [self._axis_counts[0] - nt, self._axis_counts[1] - nt, 0],
+ [[0, 0, 1]], skip_end=True, blend=True)
+
+ # average point coordinates across 3 directions between side faces and surfaces to 4 3-way lines.
+ min_weight = 1 # GRC revisit, remove?
+ # 1-direction
+ for n2 in range(1, self._box_counts[1]):
+ for n3 in range(1, self._box_counts[2]):
+ start_indexes = [0, n2, n3]
+ corner_indexes = [self._box_counts[0], n2, n3]
+ end_indexes = [self._axis_counts[0], n2, n3]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[1], None, corner[0], corner[1], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], corner[1], None, end[0], end[3], None, self._trans_count,
+ start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
+ for nt in range(1, self._trans_count):
+ start_indexes = [0, n2, self._box_counts[2] + nt]
+ corner_indexes = [self._box_counts[0] + nt, n2, self._box_counts[2] + nt]
+ end_indexes = [self._box_counts[0] + nt, n2, 0]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[1], None, corner[0], corner[1], None, self._box_counts[0],
+ start_weight=self._box_counts[0] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[1, 0, 0]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], corner[1], None, end[0], [-d for d in end[2]], None, self._box_counts[2],
+ start_weight=1.0 + min_weight, end_weight=self._box_counts[2] + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[0, 0, -1]], skip_start=True, skip_end=True, blend=True)
+ # 2-direction
+ for n3 in range(1, self._box_counts[2]):
+ for n1 in range(1, self._box_counts[0]):
+ start_indexes = [n1, 0, n3]
+ corner_indexes = [n1, self._box_counts[1], n3]
+ end_indexes = [n1, self._axis_counts[1], n3]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[2], None, corner[0], corner[2], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], corner[2], None, end[0], end[3], None, self._trans_count,
+ start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
+ for nt in range(1, self._trans_count):
+ start_indexes = [self._box_counts[0] + nt, 0, n3]
+ corner_indexes = [self._box_counts[0] + nt, self._box_counts[1] + nt, n3]
+ end_indexes = [0, self._box_counts[1] + nt, n3]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[1], None, corner[0], corner[1], None, self._box_counts[1],
+ start_weight=self._box_counts[1] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[0, 1, 0]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], corner[1], None, end[0], end[1], None, self._box_counts[0],
+ start_weight=1.0 + min_weight, end_weight=self._box_counts[0] + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[-1, 0, 0]], skip_start=True, skip_end=True, blend=True)
+ # 3-direction
+ for n1 in range(1, self._box_counts[0]):
+ for n2 in range(1, self._box_counts[1]):
+ start_indexes = [n1, n2, 0]
+ corner_indexes = [n1, n2, self._box_counts[2]]
+ end_indexes = [n1, n2, self._axis_counts[2]]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[3], None, corner[0], corner[3], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], corner[3], None, end[0], end[3], None, self._trans_count,
+ start_weight=1.0 + min_weight, end_weight=self._trans_count + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
+ for nt in range(1, self._trans_count):
+ start_indexes = [n1, self._box_counts[1] + nt, 0]
+ corner_indexes = [n1, self._box_counts[1] + nt, self._box_counts[2] + nt]
+ end_indexes = [n1, 0, self._box_counts[2] + nt]
+ start = self._nx[start_indexes[2]][start_indexes[1]][start_indexes[0]]
+ corner = self._nx[corner_indexes[2]][corner_indexes[1]][corner_indexes[0]]
+ end = self._nx[end_indexes[2]][end_indexes[1]][end_indexes[0]]
+ px, _ = sampleHermiteCurve(
+ start[0], start[2], None, corner[0], [-d for d in corner[2]], None, self._box_counts[2],
+ start_weight=self._box_counts[2] + min_weight, end_weight=1.0 + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], start_indexes, [[0, 0, 1]], skip_start=True, skip_end=True, blend=True)
+ px, _ = sampleHermiteCurve(
+ corner[0], [-d for d in corner[2]], None, end[0], [-d for d in end[2]], None, self._box_counts[1],
+ start_weight=1.0 + min_weight, end_weight=self._box_counts[1] + min_weight)
+ self._set_coordinates_across(
+ [px], [[0]], corner_indexes, [[0, -1, 0]], skip_start=True, skip_end=True, blend=True)
+
+ # smooth 1-direction
+ for n2 in range(1, self._box_counts[1]):
+ for n3 in range(1, self._box_counts[2]):
+ self._smooth_derivative_across(
+ [0, n2, n3], [self._axis_counts[0], n2, n3],
+ [[1, 0, 0]], [1, 3], fix_start_direction=True, fix_end_direction=True)
+ for nt in range(1, self._trans_count):
+ self._smooth_derivative_across(
+ [0, n2, self._box_counts[2] + nt], [self._box_counts[0] + nt, n2, 0],
+ [[1, 0, 0], [0, 0, -1]], [1, 1, -2], fix_start_direction=True, fix_end_direction=True)
+ # smooth 2-direction
+ for n3 in range(1, self._box_counts[2]):
+ for n1 in range(1, self._box_counts[0]):
+ self._smooth_derivative_across(
+ [n1, 0, n3], [n1, self._axis_counts[1], n3],
+ [[0, 1, 0]], [2, 3], fix_start_direction=True, fix_end_direction=True)
+ for nt in range(1, self._trans_count):
+ self._smooth_derivative_across(
+ [self._box_counts[0] + nt, 0, n3], [0, self._box_counts[1] + nt, n3],
+ [[0, 1, 0], [-1, 0, 0]], [1], fix_start_direction=True, fix_end_direction=True)
+ # smooth 3-direction
+ for n1 in range(1, self._box_counts[0]):
+ for n2 in range(1, self._box_counts[1]):
+ self._smooth_derivative_across(
+ [n1, n2, 0], [n1, n2, self._axis_counts[2]],
+ [[0, 0, 1]], [3], fix_start_direction=True, fix_end_direction=True)
+ for nt in range(1, self._trans_count):
+ self._smooth_derivative_across(
+ [n1, self._box_counts[1] + nt, 0], [n1, 0, self._box_counts[2] + nt],
+ [[0, 0, 1], [0, -1, 0]], [2, -2], fix_start_direction=True, fix_end_direction=True)
+
+ def mirror_yz(self):
+ """
+ Mirror coordinates and derivatives about both y = 0 and z = 0 planes.
+ """
+ for n3 in range(self._axis_counts[2] + 1):
+ nx_layer = self._nx[n3]
+ for n2 in range(self._axis_counts[1] + 1):
+ nx_row = nx_layer[n2]
+ for n1 in range(self._axis_counts[0] + 1):
+ nx = nx_row[n1]
+ if nx:
+ for d in nx:
+ if d:
+ d[1] = -d[1]
+ d[2] = -d[2]
diff --git a/src/scaffoldmaker/utils/interpolation.py b/src/scaffoldmaker/utils/interpolation.py
index 2b1c324a..2799f214 100644
--- a/src/scaffoldmaker/utils/interpolation.py
+++ b/src/scaffoldmaker/utils/interpolation.py
@@ -1791,7 +1791,7 @@ def track_curve_side_direction(cx, cd1, start_direction, start_location, end_loc
def linearlyInterpolateVectors(u, v, xi, magnitudeScalingMode=DerivativeScalingMode.ARITHMETIC_MEAN):
"""
- Linearly interpolate two vectors less than 180 degrees apart to get an in-between vector
+ Linearly interpolate two 3-component vectors less than 180 degrees apart to get an in-between vector
rotated from direction of u to direction of v proportionatly to xi, with magnitude linearly interpolated by xi.
:param u: First vector.
:param v: Second vector from 0 to less than 180 degrees away from u.
@@ -1820,7 +1820,10 @@ def linearlyInterpolateVectors(u, v, xi, magnitudeScalingMode=DerivativeScalingM
cos_phi = math.cos(phi)
sin_phi = math.sin(phi)
axis1 = dirn_u
- axis3 = normalize(cross(dirn_u, dirn_v))
+ cross_u_v = cross(dirn_u, dirn_v)
+ if magnitude(cross_u_v) == 0.0:
+ return [0.0, 0.0, 0.0]
+ axis3 = normalize(cross_u_v)
axis2 = cross(axis3, dirn_u)
dirn = add(mult(axis1, cos_phi), mult(axis2, sin_phi))
return set_magnitude(dirn, mag)
@@ -1841,7 +1844,7 @@ def sampleHermiteCurve(start_x, start_d1, start_d2, end_x, end_d1, end_d2, eleme
:param start_weight, end_weight: Optional relative weights for start/end d1.
:param overweighting: Multiplier of arc length to use with initial curve to exaggerate end derivatives.
:param end_transition: If supplied with end_d1, modify size of last element to fit end_d1.
- :return: x[], d1[], d2[]
+ :return: x[], d1[]{, d2[] if start_d2 supplied}
"""
assert (not end_transition) or end_d1
end_d1_mag = magnitude(end_d1) if end_d1 else None
diff --git a/src/scaffoldmaker/utils/meshrefinement.py b/src/scaffoldmaker/utils/meshrefinement.py
index 5b835e82..d5d2c0a5 100644
--- a/src/scaffoldmaker/utils/meshrefinement.py
+++ b/src/scaffoldmaker/utils/meshrefinement.py
@@ -1,19 +1,19 @@
"""
Class for refining a mesh from one region to another.
"""
-from __future__ import division
-
-import math
-
from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, \
findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString
+from cmlibs.utils.zinc.general import ChangeManager
from cmlibs.zinc.element import Element, Elementbasis
from cmlibs.zinc.field import Field
from cmlibs.zinc.node import Node
-from cmlibs.zinc.result import RESULT_OK as ZINC_OK
-from scaffoldmaker.annotation.annotationgroup import AnnotationGroup
+from cmlibs.zinc.result import RESULT_OK
+from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName
from scaffoldmaker.utils.octree import Octree
+import copy
+import math
+
class MeshRefinement:
"""
@@ -35,22 +35,26 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups=[]):
sourceNodes = self._sourceFm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
minimumsField = self._sourceFm.createFieldNodesetMinimum(self._sourceCoordinates, sourceNodes)
result, minimums = minimumsField.evaluateReal(self._sourceCache, 3)
- assert result == ZINC_OK, 'MeshRefinement failed to get minimum coordinates'
+ assert result == RESULT_OK, 'MeshRefinement failed to get minimum coordinates'
maximumsField = self._sourceFm.createFieldNodesetMaximum(self._sourceCoordinates, sourceNodes)
result, maximums = maximumsField.evaluateReal(self._sourceCache, 3)
- assert result == ZINC_OK, 'MeshRefinement failed to get maximum coordinates'
+ assert result == RESULT_OK, 'MeshRefinement failed to get maximum coordinates'
xrange = [(maximums[i] - minimums[i]) for i in range(3)]
edgeTolerance = 0.5 * (max(xrange))
if edgeTolerance == 0.0:
edgeTolerance = 1.0
minimums = [(minimums[i] - edgeTolerance) for i in range(3)]
maximums = [(maximums[i] + edgeTolerance) for i in range(3)]
- minimumsField = None
- maximumsField = None
+ del minimumsField
+ del maximumsField
self._sourceMesh = self._sourceFm.findMeshByDimension(3)
+ self._sourceFaceMesh = self._sourceFm.findMeshByDimension(2)
+ self._sourceLineMesh = self._sourceFm.findMeshByDimension(1)
self._sourceNodes = self._sourceFm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
self._sourceElementiterator = self._sourceMesh.createElementiterator()
self._octree = Octree(minimums, maximums)
+ self._tolerance = self._octree.getTolerance()
+ self._is_exterior_field = self._sourceFm.createFieldIsExterior()
self._targetRegion = targetRegion
self._targetFm = targetRegion.getFieldmodule()
@@ -75,15 +79,13 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups=[]):
self._sourceAnnotationGroups = sourceAnnotationGroups
self._annotationGroups = []
self._sourceAndTargetMeshGroups = []
- self._sourceAndTargetNodesetGroups = []
for sourceAnnotationGroup in sourceAnnotationGroups:
- targetAnnotationGroup = AnnotationGroup(self._targetRegion, sourceAnnotationGroup.getTerm())
+ targetAnnotationGroup = AnnotationGroup(
+ self._targetRegion, sourceAnnotationGroup.getTerm(), isMarker=sourceAnnotationGroup.isMarker())
self._annotationGroups.append(targetAnnotationGroup)
# assume have only highest dimension element or node/point annotation groups:
if sourceAnnotationGroup.hasMeshGroup(self._sourceMesh):
self._sourceAndTargetMeshGroups.append((sourceAnnotationGroup.getMeshGroup(self._sourceMesh), targetAnnotationGroup.getMeshGroup(self._targetMesh)))
- else:
- self._sourceAndTargetNodesetGroups.append((sourceAnnotationGroup.getNodesetGroup(self._sourceNodes), targetAnnotationGroup.getNodesetGroup(self._targetNodes)))
# prepare element -> marker point list map
self.elementMarkerMap = {}
@@ -121,62 +123,237 @@ def __del__(self):
def getAnnotationGroups(self):
return self._annotationGroups
- def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, numberInXi3,
- addNewNodesToOctree=True, shareNodeIds=None, shareNodeCoordinates=None):
+ def _face_add_line_ids_ending_in_x(self, face, x, line_ids: set):
+ """
+ Add identifiers of lines on face element with either end's coordinates within tolerance of x to supplied set.
+ :param face: Zinc 2-D face element.
+ :param x: 3 component coordinates list.
+ :param line_ids: Set of line identifiers.
+ """
+ for f in range(face.getNumberOfFaces()):
+ line = face.getFaceElement(f + 1)
+ if line.isValid():
+ line_id = line.getIdentifier()
+ if line_id not in line_ids:
+ # add line if it has coordinates within tolerance of x at either end
+ for xi in (0.0, 1.0):
+ self._sourceCache.setMeshLocation(line, [xi])
+ result, line_x = self._sourceCoordinates.evaluateReal(self._sourceCache, 3)
+ if result == RESULT_OK:
+ for c, line_c in zip(x, line_x):
+ if math.fabs(c - line_c) > self._tolerance:
+ break
+ else:
+ line_ids.add(line_id)
+ break
+
+ def _get_connected_exterior_face_ids(self, faces, x):
+ """
+ Get connected exterior face element identifiers with common lines to any pairs in faces list.
+ :param faces: List of at least 2 faces with common lines.
+ :param x: Coordinates of point; used only for > 2 faces.
+ :return: List of exterior boundary face identifiers from lowest to highest.
+ """
+ initial_face_count = len(faces)
+ assert initial_face_count > 1
+
+ # add faces passed in to set
+ # get common lines between them
+ face_ids = set()
+ face_line_ids = []
+ for face in faces:
+ face_ids.add(face.getIdentifier())
+ tmp_line_ids = []
+ for i in range(face.getNumberOfFaces()):
+ line = face.getFaceElement(i + 1)
+ if line.isValid():
+ tmp_line_ids.append(line.getIdentifier())
+ face_line_ids.append(tmp_line_ids)
+ new_line_ids = set()
+ # if there is a single line between 2 faces can do less work later, but not if there are collapsed faces
+ single_line = initial_face_count < 3
+ for f1 in range(len(faces) - 1):
+ for f2 in range(f1 + 1, len(faces)):
+ add_count = 0
+ for line_id in face_line_ids[f1]:
+ if line_id in face_line_ids[f2]:
+ new_line_ids.add(line_id)
+ add_count += 1
+ if add_count == 0:
+ # assume collapsed face, so add all lines from both faces ending in x at either end
+ single_line = False
+ for fi in (f1, f2):
+ self._face_add_line_ids_ending_in_x(faces[fi], x, new_line_ids)
+ line_ids = copy.copy(new_line_ids)
+
+ while True:
+ # ensure all parent elements of common lines are in face_ids set
+ new_face_ids = []
+ for line_id in new_line_ids:
+ line = self._sourceLineMesh.findElementByIdentifier(line_id)
+ for p in range(line.getNumberOfParents()):
+ face = line.getParentElement(p + 1)
+ face_id = face.getIdentifier()
+ if face_id not in face_ids:
+ new_face_ids.append(face_id)
+ face_ids.update(new_face_ids)
+
+ if single_line or (not new_face_ids):
+ break
+
+ new_line_ids.clear()
+ for face_id in new_face_ids:
+ face = self._sourceFaceMesh.findElementByIdentifier(face_id)
+ self._face_add_line_ids_ending_in_x(face, x, new_line_ids)
+ line_ids.update(new_line_ids)
+
+ exterior_face_ids = []
+ for face_id in face_ids:
+ face = self._sourceFaceMesh.findElementByIdentifier(face_id)
+ self._sourceCache.setElement(face)
+ result, value = self._is_exterior_field.evaluateReal(self._sourceCache, 1)
+ if (result == RESULT_OK) and (value != 0.0):
+ exterior_face_ids.append(face_id)
+
+ exterior_face_ids.sort()
+ return exterior_face_ids
+
+ cube_mid_face_xi = [
+ [0.0, 0.5, 0.5],
+ [1.0, 0.5, 0.5],
+ [0.5, 0.0, 0.5],
+ [0.5, 1.0, 0.5],
+ [0.5, 0.5, 0.0],
+ [0.5, 0.5, 1.0]
+ ]
+
+ square_mid_edge_xi = [
+ [0.0, 0.5],
+ [1.0, 0.5],
+ [0.5, 0.0],
+ [0.5, 1.0]
+ ]
+
+ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, numberInXi3):
"""
Refine cube sourceElement to numberInXi1*numberInXi2*numberInXi3 linear cube
sub-elements, evenly spaced in xi.
- :param addNewNodesToOctree: If True (default) add newly created nodes to
- octree to be found when refining later elements. Set to False when nodes are at the
- same location and not intended to be shared.
- :param shareNodeIds, shareNodeCoordinates: Arrays of identifiers and coordinates of
- nodes which may be shared in refining this element. If supplied, these are preferentially
- used ahead of points in the octree. Used to control merging with known nodes, e.g.
- those returned by this function for elements which used addNewNodesToOctree=False.
:return: Node identifiers, node coordinates used in refinement of sourceElement.
"""
- assert (shareNodeIds and shareNodeCoordinates) or (not shareNodeIds and not shareNodeCoordinates), \
- 'refineElementCubeStandard3d. Must supply both of shareNodeIds and shareNodeCoordinates, or neither'
- shareNodesCount = len(shareNodeIds) if shareNodeIds else 0
+ # element_identifier = sourceElement.getIdentifier()
+ # print("refine element", element_identifier)
meshGroups = []
for sourceAndTargetMeshGroup in self._sourceAndTargetMeshGroups:
if sourceAndTargetMeshGroup[0].containsElement(sourceElement):
meshGroups.append(sourceAndTargetMeshGroup[1])
+
+ faces = [None] * 6
+ exterior_faces = [False] * 6 # whether face is on exterior boundary of mesh
+ null_face_count = 0
+ if sourceElement.getNumberOfFaces() == 6:
+ for f in range(6):
+ face = sourceElement.getFaceElement(f + 1)
+ if face and face.isValid():
+ faces[f] = face
+ self._sourceCache.setElement(face)
+ result, value = self._is_exterior_field.evaluateReal(self._sourceCache, 1)
+ exterior_faces[f] = (result == RESULT_OK) and (value != 0.0)
+ else:
+ faces[f] = None
+ null_face_count += 1
+ # collapsed elements have no face, so get all faces which are adjacent to collapsed face
+ # check there is at least one valid face and one null face
+ if 0 < null_face_count < 6:
+ for f, face in enumerate(faces):
+ if not face:
+ # get coordinates at centre of face
+ self._sourceCache.setMeshLocation(sourceElement, self.cube_mid_face_xi[f])
+ result, mid_face_x = self._sourceCoordinates.evaluateReal(self._sourceCache, 3)
+ if result != RESULT_OK:
+ continue
+ # get list of all faces with mid-edge coordinate within tolerance of mid_face_x
+ adjacent_faces = []
+ for f2, face2 in enumerate(faces):
+ if face2 and not isinstance(face2, list):
+ for xi in self.square_mid_edge_xi:
+ self._sourceCache.setMeshLocation(face2, xi)
+ result, mid_edge_x = self._sourceCoordinates.evaluateReal(self._sourceCache, 3)
+ if result != RESULT_OK:
+ continue
+ for c in range(3):
+ if math.fabs(mid_edge_x[c] - mid_face_x[c]) > self._tolerance:
+ break
+ else:
+ adjacent_faces.append(face2)
+ if exterior_faces[f2]:
+ exterior_faces[f] = True
+ break
+ if adjacent_faces:
+ faces[f] = adjacent_faces
+
+ # 6 faces above + 1 extra face for not_a_face_index
+ not_a_face_index = 6
+ faces.append(None)
+ exterior_faces.append(False)
+
# create nodes
nids = []
nx = []
xi = [0.0, 0.0, 0.0]
tol = self._octree._tolerance
for k in range(numberInXi3 + 1):
- kExterior = (k == 0) or (k == numberInXi3)
+ k_face_index = 4 if (k == 0) else 5 if (k == numberInXi3) else not_a_face_index
+ k_face = faces[k_face_index]
+ k_exterior = exterior_faces[k_face_index]
xi[2] = k / numberInXi3
for j in range(numberInXi2 + 1):
- jExterior = kExterior or (j == 0) or (j == numberInXi2)
+ j_face_index = 2 if (j == 0) else 3 if (j == numberInXi2) else not_a_face_index
+ j_face = faces[j_face_index]
+ j_exterior = exterior_faces[j_face_index]
xi[1] = j / numberInXi2
for i in range(numberInXi1 + 1):
- iExterior = jExterior or (i == 0) or (i == numberInXi1)
+ i_face_index = 0 if (i == 0) else 1 if (i == numberInXi1) else not_a_face_index
+ i_face = faces[i_face_index]
+ i_exterior = exterior_faces[i_face_index]
xi[0] = i / numberInXi1
self._sourceCache.setMeshLocation(sourceElement, xi)
result, x = self._sourceCoordinates.evaluateReal(self._sourceCache, 3)
- # only exterior points are ever common:
+ shareable = False
+ surface_face_ids = True # since None is used for no extra data in Octree
+
+ connected_faces = []
+ for face in (i_face, j_face, k_face):
+ if face:
+ for tmp_face in face if isinstance(face, list) else [face]:
+ if tmp_face not in connected_faces:
+ connected_faces.append(tmp_face)
+ face_count = len(connected_faces)
+ if face_count > 0:
+ shareable = True
+ exterior_count = (i_exterior, j_exterior, k_exterior).count(True)
+ if face_count == 1:
+ if exterior_count == 1:
+ shareable = False # nodes only belong to this element
+ # else interior
+ else:
+ surface_face_ids = self._get_connected_exterior_face_ids(connected_faces, x)
+ if not surface_face_ids:
+ surface_face_ids = True
nodeId = None
- if iExterior:
- if shareNodeIds:
- for n in range(shareNodesCount):
- if (math.fabs(shareNodeCoordinates[n][0] - x[0]) <= tol) and \
- (math.fabs(shareNodeCoordinates[n][1] - x[1]) <= tol) and \
- (math.fabs(shareNodeCoordinates[n][2] - x[2]) <= tol):
- nodeId = shareNodeIds[n]
- break
- if nodeId is None:
- nodeId = self._octree.findObjectByCoordinates(x)
+ if shareable:
+ nodeId, extra_data = self._octree.findObjectByCoordinates(x, surface_face_ids)
+ # if nodeId:
+ # print("Found existing node", nodeId, extra_data, "at", x)
if nodeId is None:
node = self._targetNodes.createNode(self._nodeIdentifier, self._nodetemplate)
self._targetCache.setNode(node)
result = self._targetCoordinates.setNodeParameters(self._targetCache, -1, Node.VALUE_LABEL_VALUE, 1, x)
nodeId = self._nodeIdentifier
- if iExterior and addNewNodesToOctree:
- self._octree.addObjectAtCoordinates(x, nodeId)
+ if shareable:
+ # print("Add shareable node", nodeId, surface_face_ids, "at", x)
+ self._octree.addObjectAtCoordinates(x, (nodeId, surface_face_ids))
+ # else:
+ # print("Add unique node", nodeId, surface_face_ids, "at", x)
self._nodeIdentifier += 1
nids.append(nodeId)
nx.append(x)
@@ -192,7 +369,7 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n
enids = [nids[bni], nids[bni + 1], nids[bni + oj], nids[bni + oj + 1],
nids[bni + ok], nids[bni + ok + 1], nids[bni + ok + oj], nids[bni + ok + oj + 1]]
result = element.setNodesByIdentifier(self._targetEft, enids)
- # if result != ZINC_OK:
+ # if result != RESULT_OK:
# print('Element', self._elementIdentifier, result, enids)
self._elementIdentifier += 1
@@ -208,10 +385,10 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n
targetXi = [0.0] * 3
for marker in markerList:
markerName, sourceXi, sourceNodeIdentifier = marker
- sourceNode = self._sourceNodes.findNodeByIdentifier(sourceNodeIdentifier)
- node = self._targetMarkerNodes.createNode(self._nodeIdentifier, self._targetMarkerTemplate)
- self._targetCache.setNode(node)
- self._targetMarkerName.assignString(self._targetCache, markerName)
+ annotationGroup = findAnnotationGroupByName(self._annotationGroups, markerName)
+ if not annotationGroup:
+ print("Could not find annotation group", markerName)
+ continue
# determine which sub-element, targetXi that sourceXi maps to
targetElementIdentifier = startElementIdentifier
for i in range(3):
@@ -224,12 +401,19 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n
targetXi[i] = 1.0
targetElementIdentifier += el * elementOffset[i]
targetElement = self._targetMesh.findElementByIdentifier(targetElementIdentifier)
- result = self._targetMarkerLocation.assignMeshLocation(self._targetCache, targetElement, targetXi)
+ markerNode = annotationGroup.createMarkerNode(
+ self._nodeIdentifier, element=targetElement, xi=targetXi)
+ # add marker node to target annotation groups for any non-marker source annotation groups
+ # the source node was in. (Marker annotation groups all share the same 'marker' group.)
+ for sourceAnnotationGroup in self._sourceAnnotationGroups:
+ if not sourceAnnotationGroup.isMarker():
+ if sourceAnnotationGroup.getNodesetGroup(self._sourceNodes).findNodeByIdentifier(
+ sourceNodeIdentifier).isValid():
+ targetAnnotationGroup = findAnnotationGroupByName(
+ self._annotationGroups, sourceAnnotationGroup.getName())
+ if targetAnnotationGroup:
+ targetAnnotationGroup.getNodesetGroup(self._targetNodes).addNode(markerNode)
self._nodeIdentifier += 1
- # add new node to matching annotation groups the previous one was in
- for sourceAndTargetNodesetGroup in self._sourceAndTargetNodesetGroups:
- if sourceAndTargetNodesetGroup[0].containsNode(sourceNode):
- sourceAndTargetNodesetGroup[1].addNode(node)
return nids, nx
diff --git a/src/scaffoldmaker/utils/octree.py b/src/scaffoldmaker/utils/octree.py
index ebde76ec..f8cd87c8 100644
--- a/src/scaffoldmaker/utils/octree.py
+++ b/src/scaffoldmaker/utils/octree.py
@@ -1,6 +1,6 @@
-'''
+"""
Octree for searching for objects by coordinates
-'''
+"""
from __future__ import division
import copy
@@ -8,16 +8,16 @@
class Octree:
- '''
+ """
Octree for searching for objects by coordinates
- '''
+ """
def __init__(self, minimums, maximums, tolerance = None):
- '''
+ """
:param minimums: List of 3 minimum coordinate values. Caller to include any edge allowance.
:param maximums: List of 3 maximum coordinate values. Caller to include any edge allowance.
:param tolerance: If supplied, tolerance to use, or None to compute as 1.0E-6*diagonal.
- '''
+ """
self._dimension = 3
self._dimensionPower2 = 1 << self._dimension
self._maxObjects = 20
@@ -34,26 +34,29 @@ def __init__(self, minimums, maximums, tolerance = None):
# exactly 2^self._dimension children, cycling in lowest x index fastest
self._children = None
-
- def _findObjectByCoordinates(self, x):
- '''
+ def _findObjectByCoordinates(self, x, extra_data):
+ """
Find closest existing object with |x - ox| < tolerance.
:param x: 3 coordinates in a list.
+ :param extra_data: Extra data to compare with 2nd component of stored tuple (object, extra_data) or None if
+ not a tuple with extra data. Value must resolve to True, or be None.
:return: nearest distance, nearest object or None, None if none found.
- '''
+ """
nearestDistance = None
- nearestObject = None
+ nearestObject = (None, None) if extra_data else None
if self._coordinatesObjects is not None:
for coordinatesObject in self._coordinatesObjects:
# cheaply determine if in 2*tolerance sized box around object
- inBox = True
+ cox = coordinatesObject[0]
for c in range(self._dimension):
- if math.fabs(x[c] - coordinatesObject[0][c]) > self._tolerance:
- inBox = False
+ if math.fabs(x[c] - cox[c]) > self._tolerance:
break
- if inBox:
+ else:
+ if extra_data:
+ if extra_data != coordinatesObject[1][1]:
+ continue # extra data does not match
# now test exact distance
- distance = math.sqrt(sum((x[i] - coordinatesObject[0][i])*(x[i] - coordinatesObject[0][i]) for i in range(self._dimension)))
+ distance = math.sqrt(sum((x[c] - cox[c])*(x[c] - cox[c]) for c in range(self._dimension)))
if (distance < self._tolerance) and ((nearestDistance is None) or (distance < nearestDistance)):
nearestDistance = distance
nearestObject = coordinatesObject[1]
@@ -70,31 +73,33 @@ def _findObjectByCoordinates(self, x):
inBoundsPlusTolerance = False
break
if inBoundsPlusTolerance:
- distance, obj = self._children[i]._findObjectByCoordinates(x)
+ distance, obj = self._children[i]._findObjectByCoordinates(x, extra_data)
if (distance is not None) and ((nearestDistance is None) or (distance < nearestDistance)):
nearestDistance = distance
nearestObject = obj
return nearestDistance, nearestObject
-
- def findObjectByCoordinates(self, x):
- '''
+ def findObjectByCoordinates(self, x, extra_data=None):
+ """
Find closest existing object with |x - ox| < tolerance.
:param x: 3 coordinates in a list.
- :return: nearest object or None if not found.
- '''
- nearestDistance, nearestObject = self._findObjectByCoordinates(x)
+ :param extra_data: Optional extra data to compare with 2nd component of stored tuple (object, extra_data).
+ Default/None means no tuple, no extra data.
+ :return: nearest object (or object tuple) or None (or (None, None) if extra_data) if not found.
+ """
+ nearestDistance, nearestObject = self._findObjectByCoordinates(x, extra_data)
return nearestObject
def addObjectAtCoordinates(self, x, obj):
- '''
+ """
Add object at coordinates to octree.
Caller must have received None result for findObjectByCoordinates() first!
Assumes caller has verified x is within range of Octree.
:param x: 3 coordinates in a list.
- :param obj: object to store with coordinates.
- '''
+ :param obj: object to store with coordinates. Must be a tuple of (object, extra data) if needing to match
+ extra data when searching octree.
+ """
if self._coordinatesObjects is not None:
if len(self._coordinatesObjects) < self._maxObjects:
self._coordinatesObjects.append( (copy.deepcopy(x), obj) )
@@ -124,3 +129,6 @@ def addObjectAtCoordinates(self, x, obj):
if x[c] > centre[c]:
i += 1 << c
self._children[i].addObjectAtCoordinates(x, obj)
+
+ def getTolerance(self):
+ return self._tolerance
diff --git a/src/scaffoldmaker/utils/quadtrianglemesh.py b/src/scaffoldmaker/utils/quadtrianglemesh.py
index 012be6e6..ad7429e5 100644
--- a/src/scaffoldmaker/utils/quadtrianglemesh.py
+++ b/src/scaffoldmaker/utils/quadtrianglemesh.py
@@ -1,12 +1,9 @@
"""
-Utilities for building 3-D triangle-topology meshes out of quad elements
+Utilities for building 2-D triangle-shaped meshes out of quad elements with cubic Hermite serendipity interpolation.
"""
-from cmlibs.maths.vectorops import add, cross, div, dot, magnitude, mult, normalize, set_magnitude
from scaffoldmaker.utils.interpolation import (
- DerivativeScalingMode, get_nway_point, getCubicHermiteCurvesLength, linearlyInterpolateVectors,
- smoothCubicHermiteDerivativesLine)
+ DerivativeScalingMode, get_nway_point, linearlyInterpolateVectors, smoothCubicHermiteDerivativesLine)
import copy
-import math
class QuadTriangleMesh:
@@ -411,30 +408,50 @@ def _smooth_derivative_across(self, start_indexes, end_indexes, index_increments
pix = abs(spix)
self._nx[indexes[1]][indexes[0]][pix] = new_derivative
- def build(self):
+ def build(self, regular_count2=0):
"""
Determine interior coordinates from edge coordinates.
+ :param regular_count2: Optional number of rows of regular box elements in 2 direction.
"""
+ assert 0 <= regular_count2 <= (self._box_count2 - 1)
+ if regular_count2:
+ # sample regular row on common boundary with triangle
+ pointr1 = self._nx[0][regular_count2]
+ pointr2 = self._nx[self._element_count13][regular_count2]
+ rx, rd2, rd1 = self._sample_curve(
+ pointr1[0], pointr1[2], pointr1[1],
+ pointr2[0], [-d for d in pointr2[1]], pointr2[2],
+ self._element_count13)
+ self._set_coordinates_across([rx, rd1, rd2], [[0, 1, 2], [0, 2, 1]], [regular_count2, 0], [[0, 1]],
+ skip_start=True, skip_end=True)
+ # sample point coordinates of regular rows parallel to triangle
+ for r in range(1, regular_count2):
+ pointr1 = self._nx[0][r]
+ pointr2 = self._nx[self._element_count13][r]
+ px, _ = self._sample_curve(
+ pointr1[0], pointr1[2], None, pointr2[0], [-d for d in pointr2[1]], None, self._element_count13)
+ self._set_coordinates_across([px], [[0]], [r, 0], [[0, 1]], skip_start=True, skip_end=True)
+
# determine 3-way point location from mean curves between side points linking to it
point12 = self._nx[0][self._box_count2]
- point13 = self._nx[self._box_count3][0]
+ point13 = self._nx[self._box_count3][regular_count2]
point23 = self._nx[self._element_count13][self._element_count12]
-
x_3way, d_3way = get_nway_point(
[point23[0], point13[0], point12[0]],
[point23[1], point13[2], point12[2]],
- [self._box_count1, self._box_count2, self._box_count3],
+ [self._box_count1, self._box_count2 - regular_count2, self._box_count3],
self._sample_curve, self._move_x_to_surface,
nway_d_factor=self._3_way_d_factor)
# smooth sample from sides to 3-way points using end derivatives
- min_weight = 1 # GRC revisit, remove?
+ min_weight = 1.0 # arbitrary but looks good
ax, ad1, ad2 = self._sample_curve(
point23[0], point23[1], point23[2], x_3way, d_3way[0], None, self._box_count1,
start_weight=self._box_count1 + min_weight, end_weight=1.0 + min_weight, end_transition=True)
bx, bd2, bd1 = self._sample_curve(
- point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_count2,
- start_weight=self._box_count2 + min_weight, end_weight=1.0 + min_weight, end_transition=True)
+ point13[0], point13[2], point13[1], x_3way, d_3way[1], None, self._box_count2 - regular_count2,
+ start_weight=self._box_count2 - regular_count2 + min_weight, end_weight=1.0 + min_weight,
+ end_transition=True)
cx, cd2, cd1 = self._sample_curve(
point12[0], point12[2], point12[1], x_3way, d_3way[2], None, self._box_count3,
start_weight=self._box_count3 + min_weight, end_weight=1.0 + min_weight, end_transition=True)
@@ -442,13 +459,13 @@ def build(self):
bd1[-1] = ad1[-1]
self._set_coordinates_across([ax, ad1, ad2], [[0, 1, 2]], [self._element_count12, self._element_count13],
[[-1, -1]])
- self._set_coordinates_across([bx, bd1, bd2], [[0, 1, 2]], [0, self._box_count3], [[1, 0]])
+ self._set_coordinates_across([bx, bd1, bd2], [[0, 1, 2]], [regular_count2, self._box_count3], [[1, 0]])
self._set_coordinates_across([cx, cd1, cd2], [[0, 1, 2]], [self._box_count2, 0], [[0, 1]], skip_end=True)
- # average point coordinates across 2 directions between edges and 3-way lines.
+ # average point coordinates across 2 directions between triangle edges and 3-way lines.
# 1-2 curves
- min_weight = 1 # GRC revisit, remove?
- start_indexes = [0, 0]
+ min_weight = 1.0 # arbitrary but looks good
+ start_indexes = [regular_count2, 0]
corner_indexes = [self._box_count2, 0]
end_indexes = [self._element_count12, 0]
for i in range(1, self._box_count3):
@@ -459,8 +476,8 @@ def build(self):
corner = self._nx[corner_indexes[1]][corner_indexes[0]]
end = self._nx[end_indexes[1]][end_indexes[0]]
px, _ = self._sample_curve(
- start[0], start[1], None, corner[0], corner[1], None, self._box_count2,
- start_weight=self._box_count2 + min_weight, end_weight=1.0 + min_weight)
+ start[0], start[1], None, corner[0], corner[1], None, self._box_count2 - regular_count2,
+ start_weight=self._box_count2 - regular_count2 + min_weight, end_weight=1.0 + min_weight)
self._set_coordinates_across(
[px], [[0]], start_indexes, [[1, 0]], skip_start=True, skip_end=True)
px, _ = self._sample_curve(
@@ -469,10 +486,10 @@ def build(self):
self._set_coordinates_across(
[px], [[0]], corner_indexes, [[1, 0]], skip_start=True, skip_end=True)
# 1-3 curves
- start_indexes = [0, 0]
- corner_indexes = [0, self._box_count3]
- end_indexes = [0, self._element_count13]
- for i in range(1, self._box_count2):
+ start_indexes = [regular_count2, 0]
+ corner_indexes = [regular_count2, self._box_count3]
+ end_indexes = [regular_count2, self._element_count13]
+ for i in range(regular_count2 + 1, self._box_count2):
start_indexes[0] += 1
corner_indexes[0] += 1
end_indexes[0] += 1
@@ -492,7 +509,7 @@ def build(self):
# 2-3 curves
start_indexes = [self._element_count12, 0]
corner_indexes = [self._element_count12, self._element_count13]
- end_indexes = [0, self._element_count13]
+ end_indexes = [regular_count2, self._element_count13]
for i in range(1, self._box_count1):
start_indexes[0] -= 1
corner_indexes[0] -= 1
@@ -507,8 +524,8 @@ def build(self):
self._set_coordinates_across(
[px], [[0]], start_indexes, [[0, 1]], skip_start=True, skip_end=True, blend=True)
px, _ = self._sample_curve(
- corner[0], [-d for d in corner[2]], None, end[0], [-d for d in end[2]], None, self._box_count2,
- start_weight=1.0 + min_weight, end_weight=self._box_count2 + min_weight)
+ corner[0], [-d for d in corner[2]], None, end[0], [-d for d in end[2]], None, self._box_count2 - regular_count2,
+ start_weight=1.0 + min_weight, end_weight=self._box_count2 - regular_count2 + min_weight)
self._set_coordinates_across(
[px], [[0]], corner_indexes, [[-1, 0]], skip_start=True, skip_end=True, blend=True)
@@ -520,6 +537,13 @@ def build(self):
end_indexes[1] += 1
self._smooth_derivative_across(start_indexes, end_indexes, [[1, 0]], [1],
fix_start_direction=True, fix_end_direction=True)
+ if regular_count2:
+ # need to smooth d2 through regular elements to 3-way point
+ start_indexes[1] += 1
+ end_indexes[0] = regular_count2
+ end_indexes[1] += 1
+ self._smooth_derivative_across(start_indexes, end_indexes, [[1, 0]], [2],
+ fix_start_direction=True, fix_end_direction=True)
# smooth 1-3 curves
start_indexes = [0, 0]
end_indexes = [0, self._element_count13]
diff --git a/src/scaffoldmaker/utils/zinc_utils.py b/src/scaffoldmaker/utils/zinc_utils.py
index dcf2cd03..126d85aa 100644
--- a/src/scaffoldmaker/utils/zinc_utils.py
+++ b/src/scaffoldmaker/utils/zinc_utils.py
@@ -134,6 +134,27 @@ def mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers):
return
+def get_mesh_first_element_with_node(mesh, field, node):
+ """
+ Assumes all components of field have the same Elementfieldtemplate.
+ :param mesh: A Zinc Mesh or MeshGroup.
+ :param field: Field to query, since nodes are mapped by field.
+ :param node: A Zinc Node
+ :return: Zinc Element or None
+ """
+ elementiterator = mesh.createElementiterator()
+ element = elementiterator.next()
+ while element.isValid():
+ eft = element.getElementfieldtemplate(field, -1)
+ if eft.isValid():
+ node_count = eft.getNumberOfLocalNodes()
+ for n in range(1, node_count + 1):
+ if element.getNode(eft, n) == node:
+ return element
+ element = elementiterator.next()
+ return None
+
+
def get_nodeset_field_parameters(nodeset, field, only_value_labels=None):
"""
Returns parameters of field from nodes in nodeset in identifier order.
diff --git a/tests/test_ellipsoid.py b/tests/test_ellipsoid.py
index ccf1ece4..e312ae05 100644
--- a/tests/test_ellipsoid.py
+++ b/tests/test_ellipsoid.py
@@ -170,7 +170,7 @@ def test_ellipsoid_3D(self):
self.assertEqual(result, RESULT_OK)
self.assertAlmostEqual(surface_area, 27.86848567909992, delta=TOL)
# note exact ellipsoid volume is 4.0 / 3.0 * math.pi * a * b * c = 12.566370614359173
- self.assertAlmostEqual(volume, 12.557389634764395, delta=TOL)
+ self.assertAlmostEqual(volume, 12.557389634764352, delta=TOL)
for annotation_group in annotation_groups:
name = annotation_group.getName()
diff --git a/tests/test_lung.py b/tests/test_lung.py
index 384c0fa0..ddb0f4f9 100644
--- a/tests/test_lung.py
+++ b/tests/test_lung.py
@@ -4,6 +4,7 @@
from cmlibs.utils.zinc.finiteelement import evaluateFieldNodesetRange, findNodeWithName
from cmlibs.utils.zinc.general import ChangeManager
+from cmlibs.utils.zinc.group import mesh_group_to_identifier_ranges, nodeset_group_to_identifier_ranges
from cmlibs.zinc.context import Context
from cmlibs.zinc.field import Field
from cmlibs.zinc.result import RESULT_OK
@@ -12,6 +13,7 @@
from scaffoldmaker.meshtypes.meshtype_3d_lung1 import MeshType_3d_lung1
from scaffoldmaker.meshtypes.meshtype_3d_lung2 import MeshType_3d_lung2
from scaffoldmaker.meshtypes.meshtype_3d_lung3 import MeshType_3d_lung3
+from scaffoldmaker.meshtypes.meshtype_3d_lung4 import MeshType_3d_lung4
from scaffoldmaker.utils.meshrefinement import MeshRefinement
from testutils import assertAlmostEqualList, check_annotation_term_ids
@@ -930,10 +932,10 @@ def test_lung3_human(self):
fieldcache = fieldmodule.createFieldcache()
result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)
- self.assertAlmostEqual(surfaceArea, 4.887248605193398, delta=tol)
result, volume = volumeField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)
- self.assertAlmostEqual(volume, 0.301840881252553, delta=tol)
+ self.assertAlmostEqual(surfaceArea, 4.88745261466682, delta=tol)
+ self.assertAlmostEqual(volume, 0.301840862866859, delta=tol)
# check some annotationGroups:
expectedSizes3d = {
@@ -1000,23 +1002,6 @@ def test_lung3_human(self):
size = group.getMeshGroup(mesh1d).getSize()
self.assertEqual(expectedSizes1d[name], size, name)
- # # test finding a marker in scaffold
- # markerGroup = fieldmodule.findFieldByName("marker").castGroup()
- # markerNodes = markerGroup.getNodesetGroup(nodes)
- # self.assertEqual(7, markerNodes.getSize())
- # markerName = fieldmodule.findFieldByName("marker_name")
- # self.assertTrue(markerName.isValid())
- # markerLocation = fieldmodule.findFieldByName("marker_location")
- # self.assertTrue(markerLocation.isValid())
- # # test apex marker point
- # cache = fieldmodule.createFieldcache()
- # node = findNodeWithName(markerNodes, markerName, "apex of left lung")
- # self.assertTrue(node.isValid())
- # cache.setNode(node)
- # element, xi = markerLocation.evaluateMeshLocation(cache, 3)
- # self.assertEqual(40, element.getIdentifier())
- # assertAlmostEqualList(self, xi, [ 0.0, 1.0, 1.0 ], 1.0E-10)
-
# refine 2x2x2 and check result
# need to use original annotation groups to get temporaries
annotationGroups = originalAnnotationGroups
@@ -1062,23 +1047,283 @@ def test_lung3_human(self):
size = group.getMeshGroup(mesh2d).getSize()
self.assertEqual(expectedSizes2d[name] * (refineNumberOfElements ** 2), size, name)
- # # test finding a marker in refined scaffold
- # markerGroup = refineFieldmodule.findFieldByName("marker").castGroup()
- # refinedNodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
- # markerNodes = markerGroup.getNodesetGroup(refinedNodes)
- # self.assertEqual(7, markerNodes.getSize())
- # markerName = refineFieldmodule.findFieldByName("marker_name")
- # self.assertTrue(markerName.isValid())
- # markerLocation = refineFieldmodule.findFieldByName("marker_location")
- # self.assertTrue(markerLocation.isValid())
- # # test apex marker point
- # cache = refineFieldmodule.createFieldcache()
- # node = findNodeWithName(markerNodes, markerName, "apex of left lung")
- # self.assertTrue(node.isValid())
- # cache.setNode(node)
- # element, xi = markerLocation.evaluateMeshLocation(cache, 3)
- # self.assertEqual(319, element.getIdentifier())
- # assertAlmostEqualList(self, xi, [ 0.0, 1.0, 1.0 ], 1.0E-10)
+ def test_lung4_human(self):
+ """
+ Test creation of human lung scaffold lung4 with open fissures.
+ """
+ scaffold = MeshType_3d_lung4
+ parameter_set_names = scaffold.getParameterSetNames()
+ self.assertEqual(parameter_set_names,
+ ["Default", "Human 1 Coarse", "Human 1 Medium", "Human 1 Fine", "Ellipsoid"])
+ options = scaffold.getDefaultOptions("Human 1 Coarse")
+ self.assertEqual(16, len(options))
+ self.assertFalse(scaffold.checkOptions(options))
+
+ context = Context("Test")
+ region = context.getDefaultRegion()
+ fieldmodule = region.getFieldmodule()
+ self.assertTrue(region.isValid())
+ annotation_groups, _ = scaffold.generateBaseMesh(region, options)
+ base_annotation_groups = copy.copy(annotation_groups)
+ self.assertEqual(25, len(base_annotation_groups))
+ fieldmodule.defineAllFaces()
+ for annotation_group in annotation_groups:
+ annotation_group.addSubelements()
+ scaffold.defineFaceAnnotations(region, options, annotation_groups)
+ for annotation_group in annotation_groups:
+ if annotation_group not in base_annotation_groups:
+ annotation_group.addSubelements()
+ self.assertEqual(76, len(annotation_groups))
+
+ fieldcache = fieldmodule.createFieldcache()
+ coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+ self.assertTrue(coordinates.isValid())
+ nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+ self.assertEqual(416, nodes.getSize())
+ intersection_group_names = [
+ ("lower lobe of left lung", "upper lobe of left lung"),
+ ("lower lobe of right lung", "middle lobe of right lung"),
+ ("lower lobe of right lung", "upper lobe of right lung"),
+ ("middle lobe of right lung", "upper lobe of right lung")
+ ]
+ # check number of common nodes at hilum
+ with ChangeManager(fieldmodule):
+ for group_name1, group_name2 in intersection_group_names:
+ group = fieldmodule.createFieldGroup()
+ nodeset_group = group.createNodesetGroup(nodes)
+ nodeset_group.addNodesConditional(
+ fieldmodule.createFieldAnd(fieldmodule.findFieldByName(group_name1),
+ fieldmodule.findFieldByName(group_name2)))
+ self.assertEqual(3, nodeset_group.getSize())
+ del nodeset_group
+ del group
+
+ expected_mesh_sizes = {
+ 'mesh3d': (232, 0.23531518999948317),
+ 'mesh2d': (830, 9.322931842916137),
+ 'mesh1d': (1003, 110.4277660475001),
+ 'left lung.mesh3d': (116, 0.11765752113724519),
+ 'lower lobe of left lung.mesh3d': (44, 0.055984781964438096),
+ 'upper lobe of left lung.mesh3d': (72, 0.06167273917280671),
+ 'right lung.mesh3d': (116, 0.11765766886223888),
+ 'lower lobe of right lung.mesh3d': (44, 0.05598478196443805),
+ 'middle lobe of right lung.mesh3d': (24, 0.015908839863426272),
+ 'upper lobe of right lung.mesh3d': (48, 0.04576404703437447),
+ 'oblique fissure of lower lobe of left lung.mesh2d': (18, 0.23573085010104664),
+ 'oblique fissure of upper lobe of left lung.mesh2d': (20, 0.26116146713690147),
+ # base of middle lobe includes medial surface
+ 'base of middle lobe of right lung surface.mesh2d': (10, 0.11847510573774353),
+ 'horizontal fissure of middle lobe of right lung.mesh2d': (10, 0.08506039275120454),
+ 'horizontal fissure of upper lobe of right lung.mesh2d': (10, 0.08506039275120465),
+ 'oblique fissure of lower lobe of right lung.mesh2d': (18, 0.23573085010104675),
+ 'oblique fissure of middle lobe of right lung.mesh2d': (10, 0.11847510573774353),
+ 'oblique fissure of upper lobe of right lung.mesh2d': (10, 0.14268636139915797),
+ 'lateral edge of base of lower lobe of left lung.mesh1d': (3, 0.6163117802059481),
+ 'medial edge of base of lower lobe of left lung.mesh1d': (3, 0.5174206016915209),
+ 'posterior edge of lower lobe of left lung.mesh1d': (6, 0.7217563496699031),
+ 'lateral edge of base of lower lobe of right lung.mesh1d': (3, 0.6163117802059481),
+ 'medial edge of base of lower lobe of right lung.mesh1d': (3, 0.5174206016915209),
+ # base of middle lobe includes anterior edge
+ 'lateral edge of base of middle lobe of right lung.mesh1d': (7, 1.0751126241895612),
+ 'medial edge of base of middle lobe of right lung.mesh1d': (7, 0.8432549927571157),
+ 'posterior edge of lower lobe of right lung.mesh1d': (6, 0.7217563496699028)
+ }
+ TOL = 1.0E-6
+ with ChangeManager(fieldmodule):
+ one = fieldmodule.createFieldConstant(1.0)
+ for mesh_name, expected_sizes in expected_mesh_sizes.items():
+ expected_size, expected_val = expected_sizes
+ mesh_group = fieldmodule.findMeshByName(mesh_name)
+ size = mesh_group.getSize()
+ self.assertEqual(size, expected_size, msg=mesh_name)
+ # volume/area/length
+ val_integral = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh_group)
+ val_integral.setNumbersOfPoints(4)
+ result, val = val_integral.evaluateReal(fieldcache, 1)
+ self.assertEqual(result, RESULT_OK)
+ self.assertAlmostEqual(val, expected_val, msg=mesh_name, delta=TOL)
+
+ # check markers
+ marker_group = fieldmodule.findFieldByName("marker").castGroup()
+ marker_nodes = marker_group.getNodesetGroup(nodes)
+ self.assertEqual(9, marker_nodes.getSize())
+ for annotation_group in annotation_groups:
+ if annotation_group.isMarker():
+ group_field = annotation_group.getGroup()
+ self.assertEqual(group_field, marker_group)
+ null_group = fieldmodule.findFieldByName(annotation_group.getName()).castGroup()
+ self.assertFalse(null_group.isValid())
+
+ # refine 2x2x2 and check results
+ refine_region = region.createRegion()
+ refine_fieldmodule = refine_region.getFieldmodule()
+ options['Refine'] = True
+ options['Refine number of elements'] = 2
+ meshrefinement = MeshRefinement(region, refine_region, base_annotation_groups)
+ scaffold.refineMesh(meshrefinement, options)
+ refine_annotation_groups = meshrefinement.getAnnotationGroups()
+ self.assertEqual(25, len(refine_annotation_groups))
+ refine_fieldmodule.defineAllFaces()
+ for annotation_group in refine_annotation_groups:
+ annotation_group.addSubelements()
+ old_annotation_groups = copy.copy(refine_annotation_groups)
+ scaffold.defineFaceAnnotations(refine_region, options, refine_annotation_groups)
+ for annotation_group in refine_annotation_groups:
+ if annotation_group not in old_annotation_groups:
+ annotation_group.addSubelements()
+ self.assertEqual(76, len(refine_annotation_groups))
+
+ refine_fieldcache = refine_fieldmodule.createFieldcache()
+ refine_coordinates = refine_fieldmodule.findFieldByName("coordinates").castFiniteElement()
+ self.assertTrue(refine_coordinates.isValid())
+ refine_nodes = refine_fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+ self.assertEqual(2481, refine_nodes.getSize())
+ # check number of common nodes at hilum
+ with ChangeManager(refine_fieldmodule):
+ for group_name1, group_name2 in intersection_group_names:
+ group = refine_fieldmodule.createFieldGroup()
+ nodeset_group = group.createNodesetGroup(refine_nodes)
+ nodeset_group.addNodesConditional(
+ refine_fieldmodule.createFieldAnd(refine_fieldmodule.findFieldByName(group_name1),
+ refine_fieldmodule.findFieldByName(group_name2)))
+ self.assertEqual(5, nodeset_group.getSize())
+ del nodeset_group
+ del group
+
+ expected_refine_mesh_sizes = {
+ 'mesh3d': (1856, 0.230336480247037),
+ 'mesh2d': (6104, 16.332327494330084),
+ 'mesh1d': (6718, 358.30880646477294),
+ 'left lung.mesh3d': (928, 0.1151678888674793),
+ 'lower lobe of left lung.mesh3d': (352, 0.0548394588491132),
+ 'upper lobe of left lung.mesh3d': (576, 0.060328430018365395),
+ 'right lung.mesh3d': (928, 0.11516859137956122),
+ 'lower lobe of right lung.mesh3d': (352, 0.054839458849113155),
+ 'middle lobe of right lung.mesh3d': (192, 0.015593500770872093),
+ 'upper lobe of right lung.mesh3d': (384, 0.04473563175957532),
+ 'oblique fissure of lower lobe of left lung.mesh2d': (72, 0.23303251949030654),
+ 'oblique fissure of upper lobe of left lung.mesh2d': (80, 0.25693609525575456),
+ # base of middle lobe includes medial surface
+ 'base of middle lobe of right lung surface.mesh2d': (40, 0.11647413723834865),
+ 'horizontal fissure of middle lobe of right lung.mesh2d': (40, 0.08446588867418067),
+ 'horizontal fissure of upper lobe of right lung.mesh2d': (40, 0.08446588867418048),
+ 'oblique fissure of lower lobe of right lung.mesh2d': (72, 0.23303251949030673),
+ 'oblique fissure of middle lobe of right lung.mesh2d': (40, 0.11647413723834865),
+ 'oblique fissure of upper lobe of right lung.mesh2d': (40, 0.14046195801740563),
+ 'lateral edge of base of lower lobe of left lung.mesh1d': (6, 0.6113060310659696),
+ 'medial edge of base of lower lobe of left lung.mesh1d': (6, 0.5145764175197279),
+ 'posterior edge of lower lobe of left lung.mesh1d': (12, 0.7213630968439506),
+ 'lateral edge of base of lower lobe of right lung.mesh1d': (6, 0.6113060310659696),
+ 'medial edge of base of lower lobe of right lung.mesh1d': (6, 0.5145764175197279),
+ # base of middle lobe includes anterior edge
+ 'lateral edge of base of middle lobe of right lung.mesh1d': (14, 1.0700457132854544),
+ 'medial edge of base of middle lobe of right lung.mesh1d': (14, 0.8398626858403571),
+ 'posterior edge of lower lobe of right lung.mesh1d': (12, 0.7213630968439503)
+ }
+ with ChangeManager(refine_fieldmodule):
+ one = refine_fieldmodule.createFieldConstant(1.0)
+ for mesh_name, expected_sizes in expected_refine_mesh_sizes.items():
+ expected_size, expected_val = expected_sizes
+ mesh_group = refine_fieldmodule.findMeshByName(mesh_name)
+ size = mesh_group.getSize()
+ self.assertEqual(size, expected_size, msg=mesh_name)
+ # volume/area/length
+ val_integral = refine_fieldmodule.createFieldMeshIntegral(one, refine_coordinates, mesh_group)
+ val_integral.setNumbersOfPoints(4)
+ result, val = val_integral.evaluateReal(refine_fieldcache, 1)
+ self.assertEqual(result, RESULT_OK)
+ self.assertAlmostEqual(val, expected_val, msg=mesh_name, delta=TOL)
+
+ # check refine markers
+ refine_marker_group = refine_fieldmodule.findFieldByName("marker").castGroup()
+ refine_marker_nodes = refine_marker_group.getNodesetGroup(refine_nodes)
+ refine_lung_nodes = refine_fieldmodule.findFieldByName("lung").castGroup().getNodesetGroup(refine_nodes)
+ self.assertEqual(9, refine_marker_nodes.getSize())
+ for annotation_group in refine_annotation_groups:
+ if annotation_group.isMarker():
+ marker_node = annotation_group.getMarkerNode()
+ self.assertTrue(marker_node.isValid())
+ group_field = annotation_group.getGroup()
+ self.assertEqual(group_field, refine_marker_group)
+ self.assertTrue(refine_marker_nodes.containsNode(marker_node))
+ self.assertTrue(refine_lung_nodes.containsNode(marker_node))
+ null_group = refine_fieldmodule.findFieldByName(annotation_group.getName()).castGroup()
+ self.assertFalse(null_group.isValid())
+
+ def test_lung4_human_left(self):
+ """
+ Test creation of human lung scaffold lung4 left side only, that objects have expected identifiers.
+ """
+ scaffold = MeshType_3d_lung4
+ options = scaffold.getDefaultOptions("Human 1 Coarse")
+ options["Number of elements oblique"] = 4 # minimum, for speed
+ options["Right lung"] = False
+
+ context = Context("Test")
+ region = context.getDefaultRegion()
+ fieldmodule = region.getFieldmodule()
+ self.assertTrue(region.isValid())
+ annotation_groups, _ = scaffold.generateMesh(region, options)
+ self.assertEqual(31, len(annotation_groups))
+ # check no right annotations
+ for annotation_group in annotation_groups:
+ name = annotation_group.getName()
+ self.assertFalse("right" in name)
+ # get ranges of identifiers in meshes and nodes
+ expected_domain_count_ranges = [
+ ("mesh3d", 44, [[1, 44]]),
+ ("mesh2d", 164, [[1, 164]]),
+ ("mesh1d", 208, [[1, 208]]),
+ ("nodes", 93, [[1, 89], [187, 190]]),
+ ("marker.nodes", 4, [[187, 190]])
+ ]
+ for domain_name, expected_count, expected_ranges in expected_domain_count_ranges:
+ if "mesh" in domain_name:
+ domain = fieldmodule.findMeshByName(domain_name)
+ ranges = mesh_group_to_identifier_ranges(domain)
+ else:
+ domain = fieldmodule.findNodesetByName(domain_name)
+ ranges = nodeset_group_to_identifier_ranges(domain)
+ self.assertEqual(domain.getSize(), expected_count)
+ self.assertEqual(ranges, expected_ranges)
+
+ def test_lung4_human_right(self):
+ """
+ Test creation of human lung scaffold lung4 left side only, that objects have expected identifiers.
+ """
+ scaffold = MeshType_3d_lung4
+ options = scaffold.getDefaultOptions("Human 1 Coarse")
+ options["Number of elements oblique"] = 4 # minimum, for speed
+ options["Left lung"] = False
+
+ context = Context("Test")
+ region = context.getDefaultRegion()
+ fieldmodule = region.getFieldmodule()
+ self.assertTrue(region.isValid())
+ annotation_groups, _ = scaffold.generateMesh(region, options)
+ self.assertEqual(46, len(annotation_groups))
+ # check no right annotations
+ for annotation_group in annotation_groups:
+ name = annotation_group.getName()
+ self.assertFalse("left" in name)
+ # get ranges of identifiers in meshes and nodes
+ expected_domain_count_ranges = [
+ ("mesh3d", 44, [[45, 88]]),
+ ("mesh2d", 170, [[165, 334]]),
+ ("mesh1d", 222, [[209, 430]]),
+ ("nodes", 102, [[90, 186], [191, 195]]),
+ ("marker.nodes", 5, [[191, 195]])
+ ]
+ for domain_name, expected_count, expected_ranges in expected_domain_count_ranges:
+ if "mesh" in domain_name:
+ domain = fieldmodule.findMeshByName(domain_name)
+ ranges = mesh_group_to_identifier_ranges(domain)
+ else:
+ domain = fieldmodule.findNodesetByName(domain_name)
+ ranges = nodeset_group_to_identifier_ranges(domain)
+ self.assertEqual(domain.getSize(), expected_count, msg=domain_name)
+ self.assertEqual(ranges, expected_ranges, msg=domain_name)
+
if __name__ == "__main__":
unittest.main()
diff --git a/tests/test_uterus.py b/tests/test_uterus.py
index 756bb347..c666e111 100644
--- a/tests/test_uterus.py
+++ b/tests/test_uterus.py
@@ -121,9 +121,7 @@ def test_uterus1(self):
continue
annotationGroups.remove(annotationGroup)
self.assertEqual(22, len(annotationGroups))
- # also remove all faces and lines as not needed for refinement
- mesh2d.destroyAllElements()
- mesh1d.destroyAllElements()
+ # must keep all faces and lines as used for refinement
refineRegion = region.createRegion()
refineFieldmodule = refineRegion.getFieldmodule()