From 47b22e5d07ebc765134cbd79a40e14097b6fd9b5 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Mon, 28 Oct 2024 08:45:18 +1300 Subject: [PATCH] Add discontinuity on thorax, abdomen core boundaries --- .../meshtypes/meshtype_3d_wholebody2.py | 8 ++ src/scaffoldmaker/utils/eft_utils.py | 26 +++--- src/scaffoldmaker/utils/tubenetworkmesh.py | 81 ++++++++++++++++--- 3 files changed, 95 insertions(+), 20 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py index f965e7f8..d2137ce0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py +++ b/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 diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index 6e7366e9..d2d3e822 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/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,10 +1001,11 @@ 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 @@ -1012,11 +1013,14 @@ def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodePara 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 diff --git a/src/scaffoldmaker/utils/tubenetworkmesh.py b/src/scaffoldmaker/utils/tubenetworkmesh.py index f58aa670..f41e6770 100644 --- a/src/scaffoldmaker/utils/tubenetworkmesh.py +++ b/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): """