diff --git a/src/scaffoldmaker/annotation/bladder_terms.py b/src/scaffoldmaker/annotation/bladder_terms.py index a7c9febc..b937bc4c 100644 --- a/src/scaffoldmaker/annotation/bladder_terms.py +++ b/src/scaffoldmaker/annotation/bladder_terms.py @@ -16,10 +16,12 @@ ("dorsal part of serosa of neck of urinary bladder", "ILX:0739280"), ("dorsal part of serosa of urethra", "ILX:0739283"), ("dorsal part of serosa of urinary bladder", "ILX:0739248"), + ("left hemi-bladder", "ILX:0793661"), ("lumen of body of urinary bladder", "ILX:0739252"), ("lumen of neck of urinary bladder", "ILX:0739256"), ("lumen of urethra", "UBERON:0010390", "ILX:0736762"), ("neck of urinary bladder", "UBERON:0001258", "FMA:15912"), + ("right hemi-bladder", "ILX:0793662"), ("serosa of body of urinary bladder", "ILX:0739276"), ("serosa of neck of urinary bladder", "ILX:0739277"), ("serosa of urethra", "ILX:0739282"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index f1a3faaf..405add15 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -397,7 +397,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Neck diameter 2': 4.0, 'Wall thickness': 0.5, 'Neck angle degrees': 45, - 'Include urethra': True, + 'Include urethra': False, 'Number of elements along urethra': 8, 'Urethra diameter 1': 1.5, 'Urethra diameter 2': 1.0, @@ -925,9 +925,6 @@ def generateBaseMesh(cls, region, options): d2Final += d2List[(elementsCountThroughWall + 1) * elementsCountAround:] d3Final += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] - xFlat = d1Flat = d2Flat = [] - xOrgan = d1Organ = d2Organ = [] - # Obtain elements count along body and neck of the bladder for defining annotation groups elementsCountAlongBody = round(ureterPositionDown * elementsCountAlongBladder - 1) elementsCountAlongNeck = elementsCountAlongBladder - elementsCountAlongBody @@ -1028,6 +1025,17 @@ def generateBaseMesh(cls, region, options): urinaryBladderMeshGroup = bladderGroup.getMeshGroup(mesh) ureterMeshGroup = ureterGroup. getMeshGroup(mesh) + xOrgan = d1Organ = d2Organ = [] + + # Create ureters and urethra only for coordinates (not for flat coordinates) + if includeUrethra or includeUreter: + xFlat = d1Flat = d2Flat = [] + else: + # Obtain flat nodes coordinates + xFlat, d1Flat, d2Flat = obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, elementsCountThroughWall, + xFinal, d1Final, d2Final, bladderLength, neckDiameter1, xApexInner, d2ApexInner, + bladderWallThickness) + # Create nodes and elements nodeIdentifier, elementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xFinal, d1Final, d2Final, d3Final, xFlat, d1Flat, d2Flat, xOrgan, d1Organ, d2Organ, None, @@ -1042,8 +1050,12 @@ def generateBaseMesh(cls, region, options): # Create annotation groups for dorsal and ventral parts of the bladder and urethra bladderDorsalGroup = AnnotationGroup(region, get_bladder_term("dorsal part of bladder")) bladderVentralGroup = AnnotationGroup(region, get_bladder_term("ventral part of bladder")) + leftHemibladderGroup = AnnotationGroup(region, get_bladder_term("left hemi-bladder")) + rightHemibladderGroup = AnnotationGroup(region, get_bladder_term("right hemi-bladder")) dorsalBladderMeshGroup = bladderDorsalGroup.getMeshGroup(mesh) ventralBladderMeshGroup = bladderVentralGroup.getMeshGroup(mesh) + leftHemibladderMeshGroup = leftHemibladderGroup.getMeshGroup(mesh) + rightHemibladderMeshGroup = rightHemibladderGroup.getMeshGroup(mesh) if includeUrethra: urethraDorsalGroup = AnnotationGroup(region, get_bladder_term("dorsal part of urethra")) urethraVentralGroup = AnnotationGroup(region, get_bladder_term("ventral part of urethra")) @@ -1058,17 +1070,27 @@ def generateBaseMesh(cls, region, options): element = mesh.findElementByIdentifier(elementIdx) if e2 < elementsCountAlongBladder: ventralBladderMeshGroup.addElement(element) + if e1 < elementsCountAround // 2: + leftHemibladderMeshGroup.addElement(element) + else: + rightHemibladderMeshGroup.addElement(element) else: ventralUrethraMeshGroup.addElement(element) else: element = mesh.findElementByIdentifier(elementIdx) if e2 < elementsCountAlongBladder: dorsalBladderMeshGroup.addElement(element) + if e1 < elementsCountAround // 2: + leftHemibladderMeshGroup.addElement(element) + else: + rightHemibladderMeshGroup.addElement(element) else: dorsalUrethraMeshGroup.addElement(element) annotationGroups.append(bladderDorsalGroup) annotationGroups.append(bladderVentralGroup) + annotationGroups.append(leftHemibladderGroup) + annotationGroups.append(rightHemibladderGroup) if includeUrethra: annotationGroups.append(urethraDorsalGroup) annotationGroups.append(urethraVentralGroup) @@ -1525,3 +1547,181 @@ def generateUreterInlets(region, nodes, mesh, ureterDefaultOptions, elementsCoun mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers) return nodeIdentifier + +def obtainBladderFlatNodes(elementsCountAlongBladder, elementsCountAround, elementsCountThroughWall, + xFinal, d1Final, d2Final, bladderLength, neckDiameter1, xApexInner, d2ApexInner, + bladderWallThickness): + """ + Calculates flat coordinates for the bladder when it is opened into a flat preparation. + :param elementsCountAlongBladder: Number of elements along bladder. + :param elementsCountAround: Number of elements around bladder. + :param elementsCountThroughWall: Number of elements through wall. + :param xFinal, d1Final, d2Final: Coordinates and derivatives of all nodes of coordinates field. + :param bladderLength: Bladder length along the central path. + :param neckDiameter1: Major diameter of the bladder neck. + :param xApexInner, d2ApexInner: Coordinates and d2 of the apex located on the lumen of the bladder. + :param bladderWallThickness: Thickness of wall. + :return: Coordinates and derivatives of flat coordinates field. + """ + # Calculate ellipses circumference along the bladder inner layer(lumen) + circumferenceList = [] + minorNodeAlong_x = [] + minorNodeAlong_d2 = [] + for n1 in range(1, elementsCountAlongBladder + 1): + xNodesOnEllipse = [] + d1NodesOnEllipse = [] + idMinor = (n1 - 1) * elementsCountAround * (elementsCountThroughWall + 1) + elementsCountThroughWall + 1 + minorNodeAlong_x.append(xFinal[idMinor]) + minorNodeAlong_d2.append(d2Final[idMinor]) + for n2 in range(elementsCountAround): + id = n2 + idMinor + x = xFinal[id] + d1 = d1Final[id] + xNodesOnEllipse.append(x) + d1NodesOnEllipse.append(d1) + xNodesOnEllipse.append(xFinal[idMinor]) + d1NodesOnEllipse.append(d1Final[idMinor]) + circumference = interp.getCubicHermiteCurvesLength(xNodesOnEllipse, d1NodesOnEllipse) + circumferenceList.append(circumference) + maxCircumference = max(circumferenceList) + + # Find the angle at the bottom of the bladder neck + v1 = [0.0, 0.0, bladderLength] + v2 = [0.5 * neckDiameter1, 0.0, bladderLength] + alpha = vector.angleBetweenVectors(v1, v2) + + # Find apex to urethra arcLength in minor radius + minorNodeAlong_x.insert(0, xApexInner) + minorNodeAlong_d2.insert(0, d2ApexInner) + minorarcLength = interp.getCubicHermiteCurvesLength(minorNodeAlong_x, minorNodeAlong_d2) + + # calculate nodes coordinates based on Hammer projection formulation + xfnList = [] + angleAlongUnit = (math.pi - alpha) / elementsCountAlongBladder + angleAroundUnit = 2 * math.pi / elementsCountAround + for n2 in range(1, elementsCountAlongBladder + 1): + for n1 in range(elementsCountAround + 1): + phi = -math.pi / 2 - alpha + n2 * angleAlongUnit + if n1 < elementsCountAround // 2 + 1: + theta = n1 * angleAroundUnit + else: + theta = math.pi - n1 * angleAroundUnit + t = math.sqrt(1 + math.cos(phi) * math.cos(theta / 2)) + xScale = maxCircumference / 2 + yScale = bladderLength / 2 + x = [xScale * math.cos(phi) * math.sin(theta / 2) / t, + yScale * (math.sin(phi) / t + 1), + 0.0] + xfnList.append(x) + + # Rearrange the nodes + xfnListRearranged = [] + for n2 in range(elementsCountAlongBladder): + for n1 in range(elementsCountAround + 1): + if n1 < elementsCountAround // 2 + 1: + nodeIndex = n2 * (elementsCountAround + 1) + elementsCountAround // 2 - n1 + else: + nodeIndex = n2 * (elementsCountAround + 1) + n1 + xfnListRearranged.append(xfnList[nodeIndex]) + + # Define d1 for flat nodes + d1fnListRearranged = [] + for n2 in range(elementsCountAlongBladder): + for n1 in range(elementsCountAround + 1): + if n1 == elementsCountAround: + id = (elementsCountAround + 1) * n2 + d1 = [d1fnListRearranged[id][0], -d1fnListRearranged[id][1], d1fnListRearranged[id][2]] + else: + v1 = xfnListRearranged[(elementsCountAround + 1) * n2 + n1] + v2 = xfnListRearranged[(elementsCountAround + 1) * n2 + n1 + 1] + d1 = [v2[c] - v1[c] for c in range(3)] + d1fnListRearranged.append(d1) + # Smoothing the derivatives + smoothed_d1 = [] + for n2 in range(elementsCountAlongBladder): + xLineSmoothing = [] + d1LineSmoothing = [] + for n1 in range(elementsCountAround + 1): + xLineSmoothing.append(xfnListRearranged[n2 * (elementsCountAround + 1) + n1]) + d1LineSmoothing.append(d1fnListRearranged[n2 * (elementsCountAround + 1) + n1]) + smd1 = smoothCubicHermiteDerivativesLine(xLineSmoothing, d1LineSmoothing, fixAllDirections=False, + fixStartDerivative=True, fixEndDerivative=True, + fixStartDirection=False, fixEndDirection=False) + smoothed_d1 += smd1 + + # Define d2 for flat nodes + # Create lines from top go down the bladder + xNodesToDown = [] + d2NodesToDown = [] + for n2 in range(elementsCountAround + 1): + x = [xfnListRearranged[n1 * (elementsCountAround + 1) + n2] for n1 in range(0, elementsCountAlongBladder)] + xNodesToDown += x + for n1 in range(1, elementsCountAlongBladder + 1): + if n1 == elementsCountAlongBladder: + d2 = [0.0, minorarcLength / elementsCountAlongBladder, 0.0] + else: + v1 = xNodesToDown[elementsCountAlongBladder * n2 + n1 - 1] + v2 = xNodesToDown[elementsCountAlongBladder * n2 + n1] + d2 = [v2[c] - v1[c] for c in range(3)] + d2NodesToDown.append(d2) + # Smoothing the derivatives + smoothed_d2 = [] + for n1 in range(elementsCountAround + 1): + lineSmoothingNodes = [] + lineSmoothingNodes_d2 = [] + for n2 in range(elementsCountAlongBladder): + lineSmoothingNodes.append(xNodesToDown[n1 * elementsCountAlongBladder + n2]) + lineSmoothingNodes_d2.append(d2NodesToDown[n1 * elementsCountAlongBladder + n2]) + smd2 = smoothCubicHermiteDerivativesLine(lineSmoothingNodes, lineSmoothingNodes_d2, fixAllDirections=False, + fixStartDerivative=False, fixEndDerivative=True, + fixStartDirection=False, fixEndDirection=False) + smoothed_d2 += smd2 + # Re-arrange the derivatives order + d2fnListRearranged = [] + for n2 in range(elementsCountAlongBladder): + for n1 in range(elementsCountAround + 1): + rd2 = smoothed_d2[n1 * elementsCountAlongBladder + n2] + d2fnListRearranged.append(rd2) + + # Create the nodes through the wall + xflatList = [] + d1flatList = [] + d2flatList = [] + for n2 in range(elementsCountThroughWall + 1): + for n1 in range(len(xfnListRearranged)): + x = [xfnListRearranged[n1][0], xfnListRearranged[n1][1], n2 * bladderWallThickness / elementsCountThroughWall] + d1 = smoothed_d1[n1] + d2 = d2fnListRearranged[n1] + xflatList.append(x) + d1flatList.append(d1) + d2flatList.append(d2) + + # Apex derivatives + v1 = xApexInner + v2 = xfnListRearranged[0] + v3 = xfnListRearranged[elementsCountAround // 2] + v21 = [v2[c] - v1[c] for c in range(3)] + v31 = [v3[c] - v1[c] for c in range(3)] + d1Mag = vector.magnitude(v21) + d2Mag = vector.magnitude(v31) + + # Add apex nodes to the list + xFlat = [] + d1Flat = [] + d2Flat = [] + for n1 in range(elementsCountThroughWall + 1): + xApex = [xApexInner[0], xApexInner[1], xApexInner[2] + n1 * bladderWallThickness / elementsCountThroughWall] + xFlat.append(xApex) + d1Flat.append([-d1Mag, 0.0, 0.0]) + d2Flat.append([0.0, d2Mag, 0.0]) + + # Re-arrange the nodes + for n3 in range(1, elementsCountAlongBladder + 1): + for n2 in range(elementsCountThroughWall + 1): + for n1 in range(elementsCountAround + 1): + i = (n3 - 1) * (elementsCountAround + 1) + n2 * (elementsCountAround + 1) * elementsCountAlongBladder + n1 + xFlat.append(xflatList[i]) + d1Flat.append(d1flatList[i]) + d2Flat.append(d2flatList[i]) + + return xFlat, d1Flat, d2Flat diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 0cdbfaa0..34be1555 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -16,6 +16,7 @@ from scaffoldmaker.utils import vector from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite +from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabelsVersion from scaffoldmaker.utils.geometry import createCirclePoints @@ -536,6 +537,9 @@ def createNodesAndElements(region, elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) result = elementtemplate.defineField(coordinates, -1, eft) + elementtemplateApex = mesh.createElementtemplate() + elementtemplateApex.setElementShapeType(Element.SHAPE_TYPE_CUBE) + if xFlat: # Flat coordinates field bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) @@ -559,6 +563,14 @@ def createNodesAndElements(region, if useCrossDerivatives: flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) + flatNodetemplate3 = nodes.createNodetemplate() + flatNodetemplate3.defineField(flatCoordinates) + flatNodetemplate3.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) + flatNodetemplate3.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) + flatNodetemplate3.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + flatNodetemplate3.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + flatElementtemplate1 = mesh.createElementtemplate() flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) flatElementtemplate1.defineField(flatCoordinates, -1, eftFlat1) @@ -567,6 +579,13 @@ def createNodesAndElements(region, flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) flatElementtemplate2.defineField(flatCoordinates, -1, eftFlat2) + if closedProximalEnd: + flatElementtemplateApex1 = mesh.createElementtemplate() + flatElementtemplateApex1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + flatElementtemplateApex2 = mesh.createElementtemplate() + flatElementtemplateApex2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + if xOrgan: # Organ coordinates field bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) @@ -605,10 +624,26 @@ def createNodesAndElements(region, # Flat coordinates field if xFlat: nodeIdentifier = firstNodeIdentifier - for n2 in range(elementsCountAlong + 1): + if closedProximalEnd: + # Apexes + for n in range(elementsCountThroughWall + 1): + node = nodes.findNodeByIdentifier(nodeIdentifier) + node.merge(flatNodetemplate3) + cache.setNode(node) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xFlat[n]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Flat[n]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Flat[n]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Flat[n]) + if useCrossDerivatives: + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + for n2 in range(elementsCountAlong if closedProximalEnd else elementsCountAlong + 1): for n3 in range(elementsCountThroughWall + 1): for n1 in range(elementsCountAround): - i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 + if closedProximalEnd: + i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 + elementsCountThroughWall + 1 + else: + i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 node = nodes.findNodeByIdentifier(nodeIdentifier) node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) cache.setNode(node) @@ -642,26 +677,58 @@ def createNodesAndElements(region, nodeIdentifier = nodeIdentifier + 1 # create elements - elementtemplate3 = mesh.createElementtemplate() - elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) radiansPerElementAround = math.pi*2.0 / elementsCountAround + radiansPerElementAroundFlat = math.pi / elementsCountAround + allValueLabels = [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D2_DS1DS2] allAnnotationGroups = [] if closedProximalEnd: - # Create apex + # Create apex elements for e3 in range(elementsCountThroughWall): for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1) % elementsCountAround - eft1 = eftfactory.createEftShellPoleBottom(va * 100, vb * 100) - elementtemplate3.defineField(coordinates, -1, eft1) - element = mesh.createElement(elementIdentifier, elementtemplate3) + eftApex = eftfactory.createEftShellPoleBottom(va * 100, vb * 100) + elementtemplateApex.defineField(coordinates, -1, eftApex) + element = mesh.createElement(elementIdentifier, elementtemplateApex) bni1 = e3 + 1 bni2 = elementsCountThroughWall + 1 + elementsCountAround*e3 + e1 + 1 bni3 = elementsCountThroughWall + 1 + elementsCountAround*e3 + (e1 + 1) % elementsCountAround + 1 nodeIdentifiers = [bni1, bni2, bni3, bni1 + 1, bni2 + elementsCountAround, bni3 + elementsCountAround] - element.setNodesByIdentifier(eft1, nodeIdentifiers) + element.setNodesByIdentifier(eftApex, nodeIdentifiers) + onOpening = e1 > elementsCountAround - 2 + if xFlat: + vaf = e1 + elementsCountAround + vbf = vaf + 1 + eftApexFlat = eftfactory.createEftShellPoleBottom(vaf * 100, vbf * 100) + if e1 >= (elementsCountAround // 2): + remapEftNodeValueLabelsVersion(eftApexFlat, [1, 4], [Node.VALUE_LABEL_D_DS1], 2) + flatElementtemplateApex1.defineField(flatCoordinates, -1, eftApexFlat) + element.merge(flatElementtemplateApex1) + element.setNodesByIdentifier(eftApexFlat, nodeIdentifiers) + # set general linear map coefficients + e1b = e1 - (elementsCountAround // 2) + radiansAroundApexFlat = e1b * radiansPerElementAroundFlat + radiansAroundApexFlatNext = (e1b + 1) * radiansPerElementAroundFlat + scalefactorsFlat = [ + -1.0, + math.sin(radiansAroundApexFlat), math.cos(radiansAroundApexFlat), radiansPerElementAroundFlat, + math.sin(radiansAroundApexFlatNext), math.cos(radiansAroundApexFlatNext), radiansPerElementAroundFlat, + math.sin(radiansAroundApexFlat), math.cos(radiansAroundApexFlat), radiansPerElementAroundFlat, + math.sin(radiansAroundApexFlatNext), math.cos(radiansAroundApexFlatNext), radiansPerElementAroundFlat + ] + result = element.setScaleFactors(eftApexFlat, scalefactorsFlat) + if onOpening: + lnRemapV2 = [3, 6] + eftApexOpen = eftApexFlat + remapEftNodeValueLabelsVersion(eftApexOpen, lnRemapV2, allValueLabels, 2) + flatElementtemplateApex2.defineField(flatCoordinates, -1, eftApexOpen) + element.merge(flatElementtemplateApex2) + element.setNodesByIdentifier(eftApexOpen, nodeIdentifiers) + if eftApexOpen.getNumberOfLocalScaleFactors() > 0: + element.setScaleFactor(eftApexOpen, 1, -1.0) # set general linear map coefficients radiansAround = e1 * radiansPerElementAround radiansAroundNext = ((e1 + 1) % elementsCountAround) * radiansPerElementAround @@ -672,7 +739,7 @@ def createNodesAndElements(region, math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround ] - result = element.setScaleFactors(eft1, scalefactors) + result = element.setScaleFactors(eftApex, scalefactors) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[0] + \ annotationGroupsThroughWall[e3] diff --git a/tests/test_bladder.py b/tests/test_bladder.py index 9d51085d..84dd0d49 100644 --- a/tests/test_bladder.py +++ b/tests/test_bladder.py @@ -48,18 +48,18 @@ def test_bladderurethra1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = scaffold.generateBaseMesh(region, options) - self.assertEqual(13, len(annotationGroups)) + self.assertEqual(12, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(240, mesh3d.getSize()) + self.assertEqual(144, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(960, mesh2d.getSize()) + self.assertEqual(576, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(1201, mesh1d.getSize()) + self.assertEqual(721, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(487, nodes.getSize()) + self.assertEqual(295, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -67,7 +67,13 @@ def test_bladderurethra1(self): self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) assertAlmostEqualList(self, minimums, [-15.48570141314588, -12.992184072505665, -0.5], 1.0E-6) - assertAlmostEqualList(self, maximums, [15.485696373577879, 13.837536258199144, 127.68631532717487], 1.0E-6) + assertAlmostEqualList(self, maximums, [15.485696373577879, 12.992185118079435, 60.65303081523686], 1.0E-6) + + flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement() + self.assertTrue(flatCoordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(flatCoordinates, nodes) + assertAlmostEqualList(self, minimums, [-37.166408946939285, 0.0, 0.0], 1.0E-6) + assertAlmostEqualList(self, maximums, [43.162112471562814, 60.35206128046009, 0.5], 1.0E-6) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -76,19 +82,24 @@ def test_bladderurethra1(self): surfaceAreaField.setNumbersOfPoints(4) volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) volumeField.setNumbersOfPoints(3) + flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup) + flatSurfaceAreaField.setNumbersOfPoints(4) + fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 5334.9516480055845, delta=1.0E-6) + self.assertAlmostEqual(surfaceArea, 4564.8516929316065, delta=1.0E-6) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 2555.560965374249, delta=1.0E-6) + self.assertAlmostEqual(volume, 2221.200414557121, delta=1.0E-6) + result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(flatSurfaceArea, 4087.7670057497376, delta=1.0E-3) # check some annotationGroups: expectedSizes3d = { "dome of the bladder": 108, "neck of urinary bladder": 36, - "urethra": 96, "urinary bladder": 144 } for name in expectedSizes3d: @@ -125,7 +136,7 @@ def test_bladderurethra1(self): for annotationGroup in removeAnnotationGroups: annotationGroups.remove(annotationGroup) - self.assertEqual(13, len(annotationGroups)) + self.assertEqual(12, len(annotationGroups)) refineRegion = region.createRegion() refineFieldmodule = refineRegion.getFieldmodule() @@ -145,16 +156,16 @@ def test_bladderurethra1(self): for annotation in annotationGroups: if annotation not in oldAnnotationGroups: annotationGroup.addSubelements() - self.assertEqual(37, len(annotationGroups)) + self.assertEqual(30, len(annotationGroups)) mesh3d = refineFieldmodule.findMeshByDimension(3) - self.assertEqual(15360, mesh3d.getSize()) + self.assertEqual(9216, mesh3d.getSize()) mesh2d = refineFieldmodule.findMeshByDimension(2) - self.assertEqual(49920, mesh2d.getSize()) + self.assertEqual(29952, mesh2d.getSize()) mesh1d = refineFieldmodule.findMeshByDimension(1) - self.assertEqual(53764, mesh1d.getSize()) + self.assertEqual(32260, mesh1d.getSize()) nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(19210, nodes.getSize()) + self.assertEqual(11530, nodes.getSize()) datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize())