Merge pull request #266 from rchristie/main

Add discontinuity on thorax, abdomen core boundaries

src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py

@@ -775,9 +775,12 @@ def generateBaseMesh(cls, region, options):
 
         annotationAlongCounts = []
         annotationAroundCounts = []
+        defaultCoreBoundaryScalingMode = 1
+        annotationCoreBoundaryScalingMode = []
         for layoutAnnotationGroup in layoutAnnotationGroups:
             alongCount = 0
             aroundCount = 0
+            coreBoundaryScalingMode = 0
             name = layoutAnnotationGroup.getName()
             if "head" in name:
                 alongCount = elementsCountAlongHead
@@ -788,9 +791,11 @@ def generateBaseMesh(cls, region, options):
             elif "thorax" in name:
                 alongCount = elementsCountAlongThorax
                 aroundCount = elementsCountAroundTorso
+                coreBoundaryScalingMode = 2
             elif "abdomen" in name:
                 alongCount = elementsCountAlongAbdomen
                 aroundCount = elementsCountAroundTorso
+                coreBoundaryScalingMode = 2
             elif "arm to hand" in name:
                 alongCount = elementsCountAlongArmToHand
                 aroundCount = elementsCountAroundArm
@@ -805,6 +810,7 @@ def generateBaseMesh(cls, region, options):
                 aroundCount = elementsCountAroundLeg
             annotationAlongCounts.append(alongCount)
             annotationAroundCounts.append(aroundCount)
+            annotationCoreBoundaryScalingMode.append(coreBoundaryScalingMode)
 
         tubeNetworkMeshBuilder = BodyTubeNetworkMeshBuilder(
             networkMesh,
@@ -818,6 +824,8 @@ def generateBaseMesh(cls, region, options):
             elementsCountTransition=options['Number of elements across core transition'],
             defaultElementsCountCoreBoxMinor=options["Number of elements across core box minor"],
             annotationElementsCountsCoreBoxMinor=[],
+            defaultCoreBoundaryScalingMode=defaultCoreBoundaryScalingMode,
+            annotationCoreBoundaryScalingMode=annotationCoreBoundaryScalingMode,
             useOuterTrimSurfaces=True)
 
         meshDimension = 3

src/scaffoldmaker/utils/eft_utils.py

@@ -973,8 +973,8 @@ def getTricubicHermiteSerendipityElementNodeParameter(eft, scalefactors, nodePar
     Currently assumes extra versions are not used, but this is not checked.
     :param eft: Element field template.
     :param scalefactors: List of scale factors.
-    :param nodeParameters: List over 8 (3-D)
-    local nodes in Zinc ordering of 4 parameter vectors x, d1, d2, d3 each with 3 components.
+    :param nodeParameters: List over 8 (3-D) local nodes in Zinc ordering of 4 parameter vectors
+    x, d1, d2, d3 each with 3 components.
     :param localNodeIndex: 1-based local node index on original basis.
     :param valueLabel: A value label mapped by CubicHermiteSerendipityFunctionOffset.
     :return: Parameter.
@@ -1001,22 +1001,26 @@ def getTricubicHermiteSerendipityElementNodeParameter(eft, scalefactors, nodePar
     return parameter
 
 
-def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodeParameters, localNodeIndexes, valueLabel):
+def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodeParameters, localNodeIndexes, valueLabel,
+                                                     version=1):
     """
     Scale tricubic hermite serendipity element field template value label at local nodes to fit actual parameters.
-    Parameter must currently be an unscaled single term expression at that local node.
+    Calculation assumes value label parameter is an unscaled single term expression at the local nodes.
     :param eft: Element field template e.g. returned from determineCubicHermiteSerendipityEft.
     Note: be careful to pass a separate eft to this function from any cached which do not needing scaling!
     :param scalefactors: Existing scale factors for element
     :param nodeParameters: As passed to determineCubicHermiteSerendipityEft: list over 8 (3-D)
     local nodes in Zinc ordering of 4 parameter vectors x, d1, d2, d3 each with 3 components.
     :param localNodeIndexes: Local node indexes to scale value label at. Currently must be in [5, 6, 7, 8].
     :param valueLabel: Single value label to scale. Currently only implemened for D_DS3.
-    :return: Modified eft, scalefactors
+    :param version: Set to a number > 1 to use that derivative version instead of scale factors, and client is
+    responsible for assigning that derivative version in proportion to the returned scale factors.
+    :return: Modified eft, final scalefactors, add scalefactors (multiplying the affected derivatives; these are
+    included in final scalefactors if version == 1, otherwise client must use these to scale versioned derivatives).
     """
     assert len(nodeParameters) == 8
     assert valueLabel == Node.VALUE_LABEL_D_DS3
-    newScalefactors = copy.copy(scalefactors) if scalefactors else []
+    addScalefactors = []
     for ln in localNodeIndexes:
         assert ln in [5, 6, 7, 8]
         na = ln - 5
@@ -1027,6 +1031,10 @@ def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodePara
         db = nodeParameters[nb][3]
         dbScaled = computeCubicHermiteEndDerivative(xa, da, xb, db)
         scalefactor = magnitude(dbScaled) / magnitude(db)
-        newScalefactors.append(scalefactor)
-    addScaleEftNodesValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS3, 3)
-    return eft, newScalefactors
+        addScalefactors.append(scalefactor)
+    if version == 1:
+        newScalefactors = scalefactors + addScalefactors if scalefactors else addScalefactors
+        addScaleEftNodesValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS3, 3)
+        return eft, newScalefactors, addScalefactors
+    remapEftNodeValueLabelsVersion(eft, localNodeIndexes, [valueLabel], version)
+    return eft, scalefactors, addScalefactors

src/scaffoldmaker/utils/tubenetworkmesh.py

@@ -189,11 +189,50 @@ def getNewTrimAnnotationGroup(self):
         self._trimAnnotationGroupCount += 1
         return self.getOrCreateAnnotationGroup(("trim surface " + "{:03d}".format(self._trimAnnotationGroupCount), ""))
 
+    def resolveEftCoreBoundaryScaling(self, eft, scalefactors, nodeParameters, nodeIdentifiers, mode):
+        """
+        Resolve differences in d3 scaling across core-shell boundary by one of 2 modes.
+        Works for 8-noded tricubic Hermite serendipity basis only.
+        :param eft: Element field template to modify.
+        :param scalefactors: Existing scalefactors for use with eft.
+        :param nodeParameters: List over 8 (3-D) local nodes in Zinc ordering of 4 parameter vectors
+        x, d1, d2, d3 each with 3 components.
+        :param nodeIdentifiers: List over 8 3-D local nodes giving global node identifiers.
+        :param mode: 1 to set scale factors, 2 to add version 2 to d3 for the boundary nodes and assigning
+        values to that version equal to the scale factors x version 1.
+        :return: New eft, new scalefactors.
+        """
+        assert mode in (1, 2)
+        eft, scalefactors, addScalefactors = addTricubicHermiteSerendipityEftParameterScaling(
+            eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3, version=mode)
+        if mode == 2:
+            nodetemplate = self._nodes.createNodetemplate()
+            n = 4
+            for nodeIdentifier, scalefactor in zip(nodeIdentifiers[4:], addScalefactors):
+                node = self._nodes.findNodeByIdentifier(nodeIdentifier)
+                nodetemplate.defineFieldFromNode(self._coordinates, node)
+                versionsCount = nodetemplate.getValueNumberOfVersions(self._coordinates, -1, Node.VALUE_LABEL_D_DS3)
+                if versionsCount == 1:
+                    # make version 2 of d3 at the node and assign to it
+                    nodetemplate.setValueNumberOfVersions(self._coordinates, -1, Node.VALUE_LABEL_D_DS3, 2)
+                    # merge clears the current parameters so need to re-set
+                    node.merge(nodetemplate)
+                    self._fieldcache.setNode(node)
+                    for valueLabel, value in zip(
+                            (Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1,
+                             Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3), nodeParameters[n]):
+                        self._coordinates.setNodeParameters(self._fieldcache, -1, valueLabel, 1, value)
+                    d3_v2 = mult(nodeParameters[n][3], scalefactor)
+                    self._coordinates.setNodeParameters(self._fieldcache, -1, Node.VALUE_LABEL_D_DS3, 2, d3_v2)
+                n += 1
+        return eft, scalefactors
+
 
 class TubeNetworkMeshSegment(NetworkMeshSegment):
 
     def __init__(self, networkSegment, pathParametersList, elementsCountAround, elementsCountThroughShell,
-                 isCore=False, elementsCountCoreBoxMinor: int = 2, elementsCountTransition: int = 1):
+                 isCore=False, elementsCountCoreBoxMinor: int=2, elementsCountTransition: int=1,
+                 coreBoundaryScalingMode: int=1):
         """
         :param networkSegment: NetworkSegment this is built from.
         :param pathParametersList: [pathParameters] if 2-D or [outerPathParameters, innerPathParameters] if 3-D
@@ -210,6 +249,7 @@ def __init__(self, networkSegment, pathParametersList, elementsCountAround, elem
         self._elementsCountCoreBoxMajor = (elementsCountAround // 2) - elementsCountCoreBoxMinor
         self._elementsCountCoreBoxMinor = elementsCountCoreBoxMinor
         self._elementsCountTransition = elementsCountTransition
+        self._coreBoundaryScalingMode = coreBoundaryScalingMode
 
         assert elementsCountThroughShell > 0
         self._elementsCountThroughShell = elementsCountThroughShell
@@ -241,6 +281,9 @@ def __init__(self, networkSegment, pathParametersList, elementsCountAround, elem
         # [nAlong][nAcrossMajor][nAcrossMinor] format.
         self._boxElementIds = None  # [along][major][minor]
 
+    def getCoreBoundaryScalingMode(self):
+        return self._coreBoundaryScalingMode
+
     def getElementsCountAround(self):
         return self._elementsCountAround
 
@@ -1586,8 +1629,8 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None):
                             nodeLayouts.append(None)
                     eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts)
                     if self._elementsCountTransition == 1:
-                        eft, scalefactors = addTricubicHermiteSerendipityEftParameterScaling(
-                            eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3)
+                        eft, scalefactors = generateData.resolveEftCoreBoundaryScaling(
+                            eft, scalefactors, nodeParameters, nids, self._coreBoundaryScalingMode)
                     elementtemplate = mesh.createElementtemplate()
                     elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE)
                     elementtemplate.defineField(coordinates, -1, eft)
@@ -1625,8 +1668,8 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None):
                                 for n1 in (e1, n1p):
                                     nodeParameters.append(self.getRimCoordinates(n1, n2, n3))
                         eft = generateData.createElementfieldtemplate()
-                        eft, scalefactors = addTricubicHermiteSerendipityEftParameterScaling(
-                            eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3)
+                        eft, scalefactors = generateData.resolveEftCoreBoundaryScaling(
+                            eft, scalefactors, nodeParameters, nids, self._coreBoundaryScalingMode)
                         elementtemplateTransition = mesh.createElementtemplate()
                         elementtemplateTransition.setElementShapeType(Element.SHAPE_TYPE_CUBE)
                         elementtemplateTransition.defineField(coordinates, -1, eft)
@@ -2851,8 +2894,8 @@ def _generateTransitionElements(self, s, n2, mesh, elementtemplate, coordinates,
 
             eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts)
             if elementsCountTransition == 1:
-                eft, scalefactors = addTricubicHermiteSerendipityEftParameterScaling(
-                    eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3)
+                eft, scalefactors = generateData.resolveEftCoreBoundaryScaling(
+                    eft, scalefactors, nodeParameters, nids, segment.getCoreBoundaryScalingMode())
             elementtemplate.defineField(coordinates, -1, eft)
             element = mesh.createElement(elementIdentifier, elementtemplate)
             element.setNodesByIdentifier(eft, nids)
@@ -3037,8 +3080,8 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData):
                     if lastTransition:
                         # need to generate eft again otherwise modifying object in eftList, mucking up outer layers
                         eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts)
-                        eft, scalefactors = addTricubicHermiteSerendipityEftParameterScaling(
-                            eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3)
+                        eft, scalefactors = generateData.resolveEftCoreBoundaryScaling(
+                            eft, scalefactors, nodeParameters, nids, segment.getCoreBoundaryScalingMode())
 
                     elementtemplate.defineField(coordinates, -1, eft)
                     element = mesh.createElement(elementIdentifier, elementtemplate)
@@ -3060,6 +3103,7 @@ def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSeg
                  defaultElementsCountAround: int=8, annotationElementsCountsAround: list=[],
                  elementsCountThroughShell: int=1, isCore=False, elementsCountTransition: int=1,
                  defaultElementsCountCoreBoxMinor: int=2, annotationElementsCountsCoreBoxMinor: list=[],
+                 defaultCoreBoundaryScalingMode=1, annotationCoreBoundaryScalingMode=[],
                  useOuterTrimSurfaces=True):
         """
         Builds contiguous tube network meshes with smooth element size transitions at junctions, optionally with solid
@@ -3085,6 +3129,9 @@ def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSeg
         elements across core box minor/d3 direction for segments with any elements in the annotation group. Client must
         ensure exclusive map from segments. Groups with zero value or past end of this list use the
         defaultElementsCountCoreBoxMinor.
+        :param defaultCoreBoundaryScalingMode: Mode 1=scalefactor, 2=version.
+        :param annotationCoreBoundaryScalingMode: Boundary mode 1 or 2 in order of layoutAnnotationGroups,
+        or 0 to use default.
         :param useOuterTrimSurfaces: Set to False to use separate trim surfaces on inner and outer tubes. Ignored if
         no inner path.
         """
@@ -3097,6 +3144,8 @@ def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSeg
         self._elementsCountTransition = elementsCountTransition
         self._defaultElementsCountCoreBoxMinor = defaultElementsCountCoreBoxMinor
         self._annotationElementsCountsCoreBoxMinor = annotationElementsCountsCoreBoxMinor
+        self._defaultCoreBoundaryScalingMode = defaultCoreBoundaryScalingMode
+        self._annotationCoreBoundaryScalingMode = annotationCoreBoundaryScalingMode
         layoutFieldmodule = self._layoutRegion.getFieldmodule()
         self._layoutInnerCoordinates = layoutFieldmodule.findFieldByName("inner coordinates").castFiniteElement()
         if not self._layoutInnerCoordinates.isValid():
@@ -3114,6 +3163,7 @@ def createSegment(self, networkSegment):
         elementsCountAround = self._defaultElementsCountAround
         elementsCountCoreBoxMinor = self._defaultElementsCountCoreBoxMinor
 
+        coreBoundaryScalingMode = self._defaultCoreBoundaryScalingMode
         i = 0
         for layoutAnnotationGroup in self._layoutAnnotationGroups:
             if i >= len(self._annotationElementsCountsAround):
@@ -3135,10 +3185,19 @@ def createSegment(self, networkSegment):
                         elementsCountCoreBoxMinor = self._annotationElementsCountsCoreBoxMinor[i]
                         break
                 i += 1
-
+            i = 0
+            for layoutAnnotationGroup in self._layoutAnnotationGroups:
+                if i >= len(self._annotationCoreBoundaryScalingMode):
+                    break
+                if self._annotationCoreBoundaryScalingMode[i] > 0:
+                    if networkSegment.hasLayoutElementsInMeshGroup(
+                            layoutAnnotationGroup.getMeshGroup(self._layoutMesh)):
+                        coreBoundaryScalingMode = self._annotationCoreBoundaryScalingMode[i]
+                        break
+                i += 1
         return TubeNetworkMeshSegment(networkSegment, pathParametersList, elementsCountAround,
                                       self._elementsCountThroughShell, self._isCore, elementsCountCoreBoxMinor,
-                                      self._elementsCountTransition)
+                                      self._elementsCountTransition, coreBoundaryScalingMode)
 
     def createJunction(self, inSegments, outSegments):
         """