From f77713960d69fb2ae7db847c7e7394594e9c66cb Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 11 Oct 2021 17:13:58 +1300 Subject: [PATCH 01/21] Remove texture coordinates --- .../meshtypes/meshtype_3d_colon1.py | 14 +- .../meshtypes/meshtype_3d_colonsegment1.py | 202 ++++-------------- src/scaffoldmaker/utils/tubemesh.py | 99 ++------- 3 files changed, 62 insertions(+), 253 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 9f445a37..e45ba89b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -9,7 +9,7 @@ from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ - getTeniaColi, createFlatAndTextureCoordinatesTeniaColi, createNodesAndElementsTeniaColi + getTeniaColi, createFlatCoordinatesTeniaColi, createNodesAndElementsTeniaColi from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp @@ -628,29 +628,29 @@ def generateBaseMesh(cls, region, options): tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroupsAround, closedProximalEnd) - # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( + # Create flat coordinates + xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( xiList, relaxedLengthList, length, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) else: - # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( + # Create flat coordinates + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( xiList, relaxedLengthList, length, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 04c05ebe..c8a570e9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -321,28 +321,28 @@ def generateBaseMesh(cls, region, options): elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, tubeTCWidthList, tcThickness, sxRefList, annotationGroupsAround, closedProximalEnd) - # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( + # Create flat coordinates + xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( xiList, relaxedLengthList, segmentLength, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, transitElementList, closedProximalEnd) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) else: - # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( + # Create flat coordinates + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( xiList, relaxedLengthList, segmentLength, wallThickness, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, @@ -1177,7 +1177,7 @@ def getFullProfileFromHalfHaustrum(xHaustrumHalfSet, d1HaustrumHalfSet, def getXiListFromOuterLengthProfile(xInner, d1Inner, segmentAxis, wallThickness, transitElementList): """ - Gets a list of xi for flat and texture coordinates calculated + Gets a list of xi for flat coordinates calculated from outer arclength of elements around a segment (most relaxed state). :param xInner: Coordinates of points on inner surface around segment. :param d1Inner: Derivatives of points on inner surface around segment. @@ -1423,13 +1423,12 @@ def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, return x, d1, d2, d3 -def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, +def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, totalLengthAlong, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd): """ - Calculates flat coordinates and texture coordinates - for a colon scaffold with tenia coli when it is opened + Calculates flat coordinates for a colon scaffold with tenia coli when it is opened up into a flat preparation. :param xiList: List containing xi for each point around the outer surface of colon in its most relaxed state. @@ -1446,18 +1445,18 @@ def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, :param transitElementList: stores true if element around is an element that transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. :param closedProximalEnd: True when proximal end of tube is closed. - :return: coordinates and derivatives of flat and texture coordinates fields. + :return: coordinates and derivatives of flat coordinates fields. """ # Calculate flat coordinates factor = 3.0 if tcCount == 3 else 2.0 elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount - # Find flat and texture coordinates for colon - xFlatColon, d1FlatColon, d2FlatColon, xTextureColon, d1TextureColon, d2TextureColon \ - = tubemesh.createFlatAndTextureCoordinates(xiList, relaxedLengthList, - totalLengthAlong, wallThickness, elementsCountAround, elementsCountAlong, - elementsCountThroughWall, transitElementList) + # Find flat coordinates for colon + xFlatColon, d1FlatColon, d2FlatColon = tubemesh.createFlatCoordinates(xiList, relaxedLengthList, + totalLengthAlong, wallThickness, + elementsCountAround, elementsCountAlong, + elementsCountThroughWall, transitElementList) # Find flat coordinates for tenia coli xFlatListTC = [] @@ -1514,79 +1513,26 @@ def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, d2FlatListTC.append(d2Flat) d2FlatListTC = d2FlatListTC + d2FlatListTC[-((elementsCountAroundTC - 1)*tcCount + 1):] - # Find texture coordinates for tenia coli - wList = [] - dwList = [] - xiTexture = xiList[0] - xTextureListTC = [] - d1TextureListTC = [] - d2TextureListTC = [] - - for N in range(tcCount + 1): - idxTCMid = N*(elementsCountAroundTC + elementsCountAroundHaustrum) - TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) - TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) - dTC = xiTexture[idxTCMid] - xiTexture[TCStartIdx] if N > 0 else xiTexture[TCEndIdx] - xiTexture[idxTCMid] - v1 = [xiTexture[TCStartIdx], 0.0, 1.0] if N > 0 else [-dTC, 0.0, 1.0] - v2 = [ xiTexture[idxTCMid], 0.0, 2.0] - v3 = [ xiTexture[TCEndIdx], 0.0, 1.0] if N < tcCount else [ 1 + dTC , 0.0, 1.0] - d1 = d2 = d3 = [dTC, 0.0, 0.0] - nx = [v1, v2, v3] - nd1 = [d1, d2, d3] - sx, sd1, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) - if N > 0 and N < tcCount: - w = sx[1:-1] - dw = sd1[1:-1] if elementsCountAroundTC > 2 else nd1[1:-1] - elif N == 0: - w = sx[int(elementsCountAroundTC*0.5):-1] - dw = sd1[int(elementsCountAroundTC*0.5):-1] if elementsCountAroundTC > 2 else \ - nd1[int(elementsCountAroundTC*0.5):-1] - else: - w = sx[1:int(elementsCountAroundTC*0.5)+1] - dw = sd1[1:int(elementsCountAroundTC*0.5)+1] if elementsCountAroundTC > 2 else \ - nd1[1:int(elementsCountAroundTC*0.5)+1] - wList = wList + w - dwList = dwList + dw - - d2 = [0.0, 1.0 / elementsCountAlong, 0.0] - for n2 in range(elementsCountAlong + 1): - for n1 in range((elementsCountAroundTC-1) * tcCount + 1): - u = [ wList[n1][0], - 1.0 / elementsCountAlong * n2, - wList[n1][2]] - xTextureListTC.append(u) - d1TextureListTC.append(dwList[n1]) - d2TextureListTC.append(d2) - xFlat, d1Flat, d2Flat, _ = combineTeniaColiWithColon(xFlatColon, d1FlatColon, d2FlatColon, [], xFlatListTC, d1FlatListTC, d2FlatListTC, [], (elementsCountAroundTC - 1)*tcCount + 1, elementsCountAround + 1, elementsCountAlong, elementsCountThroughWall, closedProximalEnd) - xTexture, d1Texture, d2Texture, _ = combineTeniaColiWithColon(xTextureColon, d1TextureColon, - d2TextureColon, [], xTextureListTC, d1TextureListTC, d2TextureListTC, [], - (elementsCountAroundTC - 1)*tcCount + 1, elementsCountAround + 1, elementsCountAlong, - elementsCountThroughWall, closedProximalEnd) - - return xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture + return xFlat, d1Flat, d2Flat def createNodesAndElementsTeniaColi(region, x, d1, d2, d3, xFlat, d1Flat, d2Flat, - xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ - Create nodes and elements for the coordinates, flat coordinates, - and texture coordinates field. Note that flat and texture coordinates - not implemented for closedProximalEnd yet. + Create nodes and elements for the coordinates and flat coordinates fields. + Note that flat coordinates not implemented for closedProximalEnd yet. :param x, d1, d2, d3: coordinates and derivatives of coordinates field. :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of flat coordinates field. - :param xTexture, d1Texture, d2Texture, d3Texture: coordinates and derivatives - of texture coordinates field. :param elementsCountAroundTC: Number of elements around tenia coli. :param elementsCountAroundHaustrum: Number of elements around haustrum. :param elementsCountAlong: Number of elements along colon. @@ -1653,11 +1599,11 @@ def createNodesAndElementsTeniaColi(region, elementtemplate2.defineField(coordinates, -1, eft2) bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) - eftTexture3 = bicubichermitelinear.createEftBasic() - eftTexture4 = bicubichermitelinear.createEftOpenTube() - eftTexture5 = bicubichermitelinear.createEftWedgeXi1One() - eftTexture6 = bicubichermitelinear.createEftWedgeXi1Zero() - eftTexture7 = bicubichermitelinear.createEftWedgeXi1ZeroOpenTube() + eftFlat3 = bicubichermitelinear.createEftBasic() + eftFlat4 = bicubichermitelinear.createEftOpenTube() + eftFlat5 = bicubichermitelinear.createEftWedgeXi1One() + eftFlat6 = bicubichermitelinear.createEftWedgeXi1Zero() + eftFlat7 = bicubichermitelinear.createEftWedgeXi1ZeroOpenTube() if xFlat: # Create flat coordinates field @@ -1680,62 +1626,23 @@ def createNodesAndElementsTeniaColi(region, flatElementtemplate1 = mesh.createElementtemplate() flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate1.defineField(flatCoordinates, -1, eftTexture3) + flatElementtemplate1.defineField(flatCoordinates, -1, eftFlat3) flatElementtemplate2 = mesh.createElementtemplate() flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture4) + flatElementtemplate2.defineField(flatCoordinates, -1, eftFlat4) flatElementtemplate3 = mesh.createElementtemplate() flatElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate3.defineField(flatCoordinates, -1, eftTexture5) + flatElementtemplate3.defineField(flatCoordinates, -1, eftFlat5) flatElementtemplate4 = mesh.createElementtemplate() flatElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate4.defineField(flatCoordinates, -1, eftTexture6) + flatElementtemplate4.defineField(flatCoordinates, -1, eftFlat6) flatElementtemplate5 = mesh.createElementtemplate() flatElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate5.defineField(flatCoordinates, -1, eftTexture7) - - if xTexture: - # Create texture coordinates field - textureCoordinates = findOrCreateFieldTextureCoordinates(fm) - textureNodetemplate1 = nodes.createNodetemplate() - textureNodetemplate1.defineField(textureCoordinates) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - - textureNodetemplate2 = nodes.createNodetemplate() - textureNodetemplate2.defineField(textureCoordinates) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) - if useCrossDerivatives: - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) - - textureElementtemplate1 = mesh.createElementtemplate() - textureElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate1.defineField(textureCoordinates, -1, eftTexture3) - - textureElementtemplate2 = mesh.createElementtemplate() - textureElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate2.defineField(textureCoordinates, -1, eftTexture4) - - textureElementtemplate3 = mesh.createElementtemplate() - textureElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate3.defineField(textureCoordinates, -1, eftTexture5) - - textureElementtemplate4 = mesh.createElementtemplate() - textureElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate4.defineField(textureCoordinates, -1, eftTexture6) - - textureElementtemplate5 = mesh.createElementtemplate() - textureElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate5.defineField(textureCoordinates, -1, eftTexture7) + flatElementtemplate5.defineField(flatCoordinates, -1, eftFlat7) # create nodes for coordinates field for n in range(len(x)): @@ -1754,7 +1661,7 @@ def createNodesAndElementsTeniaColi(region, nodeIdentifier = nodeIdentifier + 1 # Create nodes for flat coordinates field - if xFlat and xTexture: + if xFlat: nodeIdentifier = firstNodeIdentifier for n2 in range(elementsCountAlong + 1): for n3 in range(elementsCountThroughWall + 1): @@ -1763,29 +1670,20 @@ def createNodesAndElementsTeniaColi(region, (elementsCountAround + 1)*n3 + n1 + n2*((elementsCountAroundTC - 1)*tcCount + 1) node = nodes.findNodeByIdentifier(nodeIdentifier) node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) - node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) cache.setNode(node) # print('NodeIdentifier', nodeIdentifier, 'version 1', xFlatList[i]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xFlat[i]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Flat[i]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Flat[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xTexture[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Texture[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Texture[i]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) if n1 == 0: # print('NodeIdentifier', nodeIdentifier, 'version 2', xFlatList[i+elementsCountAround]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xFlat[i+elementsCountAround]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Flat[i+elementsCountAround]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Flat[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xTexture[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Texture[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Texture[i+elementsCountAround]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 # Create flat coordinates nodes for tenia coli @@ -1793,27 +1691,18 @@ def createNodesAndElementsTeniaColi(region, j = i + 2 + nTC node = nodes.findNodeByIdentifier(nodeIdentifier) node.merge(flatNodetemplate2 if nTC == 0 else flatNodetemplate1) - node.merge(textureNodetemplate2 if nTC == 0 else textureNodetemplate1) cache.setNode(node) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xFlat[j]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Flat[j]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Flat[j]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xTexture[j]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Texture[j]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Texture[j]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) if nTC == 0: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xFlat[j+(elementsCountAroundTC-1)*tcCount]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Flat[j+(elementsCountAroundTC-1)*tcCount]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Flat[j+(elementsCountAroundTC-1)*tcCount]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xTexture[j+(elementsCountAroundTC-1)*tcCount]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Texture[j+(elementsCountAroundTC-1)*tcCount]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Texture[j+(elementsCountAroundTC-1)*tcCount]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 # Create elements @@ -2029,10 +1918,9 @@ def createNodesAndElementsTeniaColi(region, onOpening = e1 > elementsCountAround - 2 element = mesh.createElement(elementIdentifier, elementtemplate) element.setNodesByIdentifier(eft, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - element.merge(textureElementtemplate2 if onOpening else textureElementtemplate1) - element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture3, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat4 if onOpening else eftFlat3, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ annotationGroupsThroughWall[e3] @@ -2065,10 +1953,9 @@ def createNodesAndElementsTeniaColi(region, bni22 + now + tcOffset, bni31, bni31 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate if eTC < int(elementsCountAroundTC*0.5) - 1 else elementtemplate1) element.setNodesByIdentifier(eft if eTC < int(elementsCountAroundTC*0.5) - 1 else eft1, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate1 if eTC < int(elementsCountAroundTC*0.5) - 1 else flatElementtemplate3) - element.merge(textureElementtemplate1 if eTC < int(elementsCountAroundTC*0.5) - 1 else textureElementtemplate3) - element.setNodesByIdentifier(eftTexture3 if eTC < int(elementsCountAroundTC*0.5) - 1 else eftTexture5, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat3 if eTC < int(elementsCountAroundTC*0.5) - 1 else eftFlat5, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] if annotationGroups: @@ -2108,28 +1995,25 @@ def createNodesAndElementsTeniaColi(region, bni22 + now + tcOffset, bni32, bni32 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate2) element.setNodesByIdentifier(eft2, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate4) - element.merge(textureElementtemplate4) - element.setNodesByIdentifier(eftTexture6, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat6, nodeIdentifiers) elif eTC > 0 and eTC < elementsCountAroundTC - 1: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate) element.setNodesByIdentifier(eft, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate1) - element.merge(textureElementtemplate1) - element.setNodesByIdentifier(eftTexture3, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat3, nodeIdentifiers) else: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni31 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate1) element.setNodesByIdentifier(eft1, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate3) - element.merge(textureElementtemplate3) - element.setNodesByIdentifier(eftTexture5, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat5, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] @@ -2169,15 +2053,13 @@ def createNodesAndElementsTeniaColi(region, onOpening = (eTC == int(elementsCountAroundTC*0.5 - 1)) element = mesh.createElement(elementIdentifier, elementtemplate if eTC > 0 else elementtemplate2) element.setNodesByIdentifier(eft if eTC > 0 else eft2, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: if eTC > 0: element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - element.merge(textureElementtemplate2 if onOpening else textureElementtemplate1) - element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture3, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat4 if onOpening else eftFlat3, nodeIdentifiers) else: element.merge(flatElementtemplate5 if onOpening else flatElementtemplate4) - element.merge(textureElementtemplate5 if onOpening else textureElementtemplate4) - element.setNodesByIdentifier(eftTexture7 if onOpening else eftTexture6, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat7 if onOpening else eftFlat6, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[ elementsCountAround + int(elementsCountAroundTC * (tcCount - 0.5)) + eTC] + \ diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index cc552fd3..564ad4d0 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -335,12 +335,11 @@ def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, return xList, d1List, d2List, d3List, curvatureList -def createFlatAndTextureCoordinates(xiList, lengthAroundList, +def createFlatCoordinates(xiList, lengthAroundList, totalLengthAlong, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ - Calculates flat coordinates and texture coordinates - for a tube when it is opened into a flat preparation. + Calculates flat coordinates for a tube when it is opened into a flat preparation. :param xiList: List containing xi for each point around the outer surface of the tube. :param lengthAroundList: List of total arclength around the @@ -352,7 +351,7 @@ def createFlatAndTextureCoordinates(xiList, lengthAroundList, :param elementsCountThroughWall: Number of elements through wall. :param transitElementList: stores true if element around is a transition element between a big and small element. - :return: coordinates and derivatives of flat and texture coordinates fields. + :return: coordinates and derivatives of flat coordinates field. """ # Calculate flat coordinates and derivatives @@ -397,55 +396,20 @@ def createFlatAndTextureCoordinates(xiList, lengthAroundList, d2FlatList.append(d2Flat) d2FlatList = d2FlatList + d2FlatList[-(elementsCountAround+1)*(elementsCountThroughWall+1):] - # Calculate texture coordinates and derivatives - xTextureList = [] - d1Texture = [] - d1TextureList = [] - d2TextureList = [] - - d2 = [0.0, 1.0 / elementsCountAlong, 0.0] - xiTexture = xiList[0] - - for n1 in range(len(xiTexture)): - d1 = [xiTexture[n1] - xiTexture[n1-1] if n1 > 0 else xiTexture[n1+1] - xiTexture[n1], - 0.0, - 0.0] - d1Texture.append(d1) - - # To modify derivative along transition elements - for i in range(len(transitElementList)): - if transitElementList[i]: - d1Texture[i+1] = d1Texture[i+2] - - for n2 in range(elementsCountAlong + 1): - for n3 in range(elementsCountThroughWall + 1): - for n1 in range(elementsCountAround + 1): - u = [ xiTexture[n1], - 1.0 / elementsCountAlong * n2, - 1.0 / elementsCountThroughWall * n3] - d1 = d1Texture[n1] - xTextureList.append(u) - d1TextureList.append(d1) - d2TextureList.append(d2) - - return xFlatList, d1FlatList, d2FlatList, xTextureList, d1TextureList, d2TextureList + return xFlatList, d1FlatList, d2FlatList def createNodesAndElements(region, x, d1, d2, d3, xFlat, d1Flat, d2Flat, - xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ - Create nodes and elements for the coordinates, flat coordinates, - and texture coordinates field. + Create nodes and elements for the coordinates and flat coordinates fields. :param x, d1, d2, d3: coordinates and derivatives of coordinates field. :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of flat coordinates field. - :param xTexture, d1Texture, d2Texture, d3Texture: coordinates and derivatives - of texture coordinates field. :param elementsCountAround: Number of elements around tube. :param elementsCountAlong: Number of elements along tube. :param elementsCountThroughWall: Number of elements through wall. @@ -499,8 +463,8 @@ def createNodesAndElements(region, if xFlat: # Flat coordinates field bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) - eftTexture1 = bicubichermitelinear.createEftBasic() - eftTexture2 = bicubichermitelinear.createEftOpenTube() + eftFlat1 = bicubichermitelinear.createEftBasic() + eftFlat2 = bicubichermitelinear.createEftOpenTube() flatCoordinates = findOrCreateFieldCoordinates(fm, name="flat coordinates") flatNodetemplate1 = nodes.createNodetemplate() @@ -521,38 +485,11 @@ def createNodesAndElements(region, flatElementtemplate1 = mesh.createElementtemplate() flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate1.defineField(flatCoordinates, -1, eftTexture1) + flatElementtemplate1.defineField(flatCoordinates, -1, eftFlat1) flatElementtemplate2 = mesh.createElementtemplate() flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture2) - - if xTexture: - # Texture coordinates field - textureCoordinates = findOrCreateFieldTextureCoordinates(fm) - textureNodetemplate1 = nodes.createNodetemplate() - textureNodetemplate1.defineField(textureCoordinates) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - - textureNodetemplate2 = nodes.createNodetemplate() - textureNodetemplate2.defineField(textureCoordinates) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) - if useCrossDerivatives: - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) - - elementtemplate1 = mesh.createElementtemplate() - elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate1.defineField(textureCoordinates, -1, eftTexture1) - - elementtemplate2 = mesh.createElementtemplate() - elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate2.defineField(textureCoordinates, -1, eftTexture2) + flatElementtemplate2.defineField(flatCoordinates, -1, eftFlat2) # Create nodes # Coordinates field @@ -571,8 +508,8 @@ def createNodesAndElements(region, # print('NodeIdentifier = ', nodeIdentifier, x[n], d1[n], d2[n]) nodeIdentifier = nodeIdentifier + 1 - # Flat and texture coordinates fields - if xFlat and xTexture: + # Flat coordinates field + if xFlat: nodeIdentifier = firstNodeIdentifier for n2 in range(elementsCountAlong + 1): for n3 in range(elementsCountThroughWall + 1): @@ -580,29 +517,20 @@ def createNodesAndElements(region, i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 node = nodes.findNodeByIdentifier(nodeIdentifier) node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) - node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) cache.setNode(node) # print('NodeIdentifier', nodeIdentifier, 'version 1, xList Index =', i+1) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xFlat[i]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Flat[i]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Flat[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xTexture[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Texture[i]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Texture[i]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) if n1 == 0: # print('NodeIdentifier', nodeIdentifier, 'version 2, xList Index =', i+elementsCountAround+1) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xFlat[i+elementsCountAround]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Flat[i+elementsCountAround]) flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Flat[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xTexture[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Texture[i+elementsCountAround]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Texture[i+elementsCountAround]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 # create elements @@ -667,10 +595,9 @@ def createNodesAndElements(region, onOpening = e1 > elementsCountAround - 2 element = mesh.createElement(elementIdentifier, elementtemplate) element.setNodesByIdentifier(eft, nodeIdentifiers) - if xFlat and xTexture: + if xFlat: element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - element.merge(elementtemplate2 if onOpening else elementtemplate1) - element.setNodesByIdentifier(eftTexture2 if onOpening else eftTexture1, nodeIdentifiers) + element.setNodesByIdentifier(eftFlat2 if onOpening else eftFlat1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ From 00017017a4d762411d1f12f81f70dab6e899c1bb Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 12 Oct 2021 09:44:33 +1300 Subject: [PATCH 02/21] Remove texture coordinates in bladder, small intestine and cecum --- src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py | 3 +-- src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py | 5 ++--- src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py | 4 ++-- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 7ecbf4d0..ec1a0879 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -903,7 +903,6 @@ def generateBaseMesh(cls, region, options): d3Final += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] xFlat = d1Flat = d2Flat = d3Flat = [] - xTexture = d1Texture = d2Texture = d3Texture = [] # Obtain elements count along body and neck of the bladder for defining annotation groups elementsCountAlongBody = round(ureterPositionDown * elementsCountAlongBladder - 1) @@ -1004,7 +1003,7 @@ def generateBaseMesh(cls, region, options): # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xFinal, d1Final, d2Final, d3Final, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xFinal, d1Final, d2Final, d3Final, xFlat, d1Flat, d2Flat, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 78acad65..c7447b58 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -464,7 +464,6 @@ def generateBaseMesh(cls, region, options): d3Cecum += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] xFlat = d1Flat = d2Flat = d3Flat = [] - xTexture = d1Texture = d2Texture = d3Texture = [] # Create nodes and elements if tcThickness > 0: @@ -475,7 +474,7 @@ def generateBaseMesh(cls, region, options): tubeTCWidthList, tcThickness, sxRefList, annotationGroupsAround, closedProximalEnd) nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, @@ -483,7 +482,7 @@ def generateBaseMesh(cls, region, options): else: nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 13c4ca59..2b2d1f73 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -922,13 +922,13 @@ def generateBaseMesh(cls, region, options): flatWidthList, xiList = smallIntestineSegmentTubeMeshInnerPoints.getFlatWidthAndXiList() # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( xiList, flatWidthList, length, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, From 6aa7469953a79e750d4e2b64b0e0a63e87c1255b Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 12 Oct 2021 18:15:27 +1300 Subject: [PATCH 03/21] Add colon coordinates --- .../meshtypes/meshtype_3d_bladderurethra1.py | 5 +- .../meshtypes/meshtype_3d_cecum1.py | 7 +- .../meshtypes/meshtype_3d_colon1.py | 37 ++++- .../meshtypes/meshtype_3d_colonsegment1.py | 153 +++++++++++++++++- .../meshtypes/meshtype_3d_smallintestine1.py | 6 +- src/scaffoldmaker/utils/tubemesh.py | 85 ++++++++++ 6 files changed, 271 insertions(+), 22 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index ec1a0879..2021ce78 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -902,7 +902,8 @@ def generateBaseMesh(cls, region, options): d2Final += d2List[(elementsCountThroughWall + 1) * elementsCountAround:] d3Final += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] - xFlat = d1Flat = d2Flat = d3Flat = [] + 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) @@ -1003,7 +1004,7 @@ def generateBaseMesh(cls, region, options): # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xFinal, d1Final, d2Final, d3Final, xFlat, d1Flat, d2Flat, + region, xFinal, d1Final, d2Final, d3Final, xFlat, d1Flat, d2Flat, xOrgan, d1Organ, d2Organ, None, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index c7447b58..f411c032 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -463,7 +463,8 @@ def generateBaseMesh(cls, region, options): d2Cecum += d2List[(elementsCountThroughWall + 1) * elementsCountAround:] d3Cecum += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] - xFlat = d1Flat = d2Flat = d3Flat = [] + xFlat = d1Flat = d2Flat = [] + xOrgan = d1Organ = d2Organ = [] # Create nodes and elements if tcThickness > 0: @@ -474,7 +475,7 @@ def generateBaseMesh(cls, region, options): tubeTCWidthList, tcThickness, sxRefList, annotationGroupsAround, closedProximalEnd) nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xOrgan, d1Organ, d2Organ, None, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, @@ -482,7 +483,7 @@ def generateBaseMesh(cls, region, options): else: nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xOrgan, d1Organ, d2Organ, None, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index e45ba89b..5854e4ce 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -9,7 +9,7 @@ from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ - getTeniaColi, createFlatCoordinatesTeniaColi, createNodesAndElementsTeniaColi + getTeniaColi, createFlatCoordinatesTeniaColi, createColonCoordinatesTeniaColi, createNodesAndElementsTeniaColi from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils import interpolation as interp @@ -475,6 +475,11 @@ def generateBaseMesh(cls, region, options): useCubicHermiteThroughWall = not(segmentSettings['Use linear through wall']) elementsCountAlong = int(elementsCountAlongSegment*segmentCount) + # Colon coordinates + lengthToDiameterRatio = 24 + wallThicknessToDiameterRatio = 0.1 + teniaColiThicknessToDiameterRatio = 0.005 + firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -634,13 +639,22 @@ def generateBaseMesh(cls, region, options): elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd) + # Create colon coordinates + xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, + wallThicknessToDiameterRatio, + teniaColiThicknessToDiameterRatio, tcCount, + elementsCountAroundTC, + elementsCountAroundHaustrum, + elementsCountAlong, elementsCountThroughWall, + transitElementList, closedProximalEnd) + # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, - elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, - firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, - closedProximalEnd) + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xColon, d1Colon, d2Colon, + "colon coordinates", elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, + elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, + annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, + useCrossDerivatives, closedProximalEnd) else: # Create flat coordinates @@ -648,10 +662,17 @@ def generateBaseMesh(cls, region, options): xiList, relaxedLengthList, length, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) + # Create colon coordinates + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, lengthToDiameterRatio, + wallThicknessToDiameterRatio, + elementsCountAround, + elementsCountAlong, elementsCountThroughWall, + transitElementList) + # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, - elementsCountAround, elementsCountAlong, elementsCountThroughWall, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xColon, d1Colon, d2Colon, + "colon coordinates", elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index c8a570e9..4e474e58 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -312,6 +312,8 @@ def generateBaseMesh(cls, region, options): d2WarpedList, d3WarpedUnitList, contractedWallThicknessList, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) + xColonSegment = d1ColonSegment = d2ColonSegment = [] + relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() if tcThickness > 0: @@ -329,11 +331,11 @@ def generateBaseMesh(cls, region, options): # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, - elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, - tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, - firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, - closedProximalEnd) + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xColonSegment, d1ColonSegment, + d2ColonSegment, None, elementsCountAroundTC, elementsCountAroundHaustrum, + elementsCountAlongSegment, elementsCountThroughWall, tcCount, annotationGroupsAround, + annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) else: # Create flat coordinates xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( @@ -342,8 +344,8 @@ def generateBaseMesh(cls, region, options): # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, - elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xColonSegment, d1ColonSegment, + d2ColonSegment, None, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) @@ -1519,9 +1521,83 @@ def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, return xFlat, d1Flat, d2Flat +def createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, wallThicknessToDiameterRatio, + teniaColiThicknessToDiameterRatio, tcCount, elementsCountAroundTC, + elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, + transitElementList, closedProximalEnd): + """ + Calculates organ coordinates for a colon scaffold with tenia coli. Organ coordinate takes the form of a cylinder + with unit inner diameter, length of lengthToDiameterRatio, wall thickness of wallThicknessToDiameterRatio, with + tenia coli of teniaColiThicknessToDiameterRatio running along its length. + :param xiList: List containing xi for each point around the outer surface of colon in its most relaxed state. + :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ + :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. + :param teniaColiThicknessToDiameterRatio: Ratio of tenia coli thickness to inner diameter of organ. + :param tcCount: Number of tenia coli. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param elementsCountAlong: Number of elements along colon. + :param elementsCountThroughWall: Number of elements through wall. + :param transitElementList: stores true if element around is an element that transits from tenia coli to haustrum. + :param closedProximalEnd: True when proximal end of tube is closed. + :return: coordinates and derivatives of colon coordinates field. + """ + # Calculate organ coordinates + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount + + # Find organ coordinates for colon + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, lengthToDiameterRatio, + wallThicknessToDiameterRatio, elementsCountAround, + elementsCountAlong, elementsCountThroughWall, + transitElementList) + + # Find organ coordinates for tenia coli + xTC = [] + d1TC = [] + d2 = [0.0, lengthToDiameterRatio / elementsCountAlong, 0.0] + + for n2 in range(elementsCountAlong + 1): + faceStartIdx = elementsCountAround * (elementsCountThroughWall + 1) * n2 + \ + elementsCountAround * elementsCountThroughWall + xTCRaw = [] + d1TCRaw = [] + d2TCRaw = [] + for N in range(tcCount): + TCMidIdx = faceStartIdx + N*(elementsCountAroundTC + elementsCountAroundHaustrum) + TCStartIdx = TCMidIdx - int(elementsCountAroundTC * 0.5) if N > 0 else \ + TCMidIdx + tcCount * (elementsCountAroundTC + elementsCountAroundHaustrum) - int( + elementsCountAroundTC * 0.5) + TCEndIdx = TCMidIdx + int(elementsCountAroundTC*0.5) + v1 = xColon[TCStartIdx] + norm = vector.setMagnitude( + vector.crossproduct3(vector.normalise(d1Colon[TCMidIdx]), vector.normalise(d2Colon[TCMidIdx])), + teniaColiThicknessToDiameterRatio) + v2 = [xColon[TCMidIdx][c] + norm[c] for c in range(3)] + v3 = xColon[TCEndIdx] + nx = [v1, v2, v3] + nd1 = [d1Colon[TCStartIdx], + vector.setMagnitude(d1Colon[TCMidIdx], 2.0*vector.magnitude(d1Colon[TCMidIdx])), + d1Colon[TCEndIdx]] + sx, sd1 = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC)[0:2] + xTCRaw += sx[1:-1] + d1TCRaw += sd1[1:-1] + + xTC += xTCRaw[int(elementsCountAroundTC * 0.5 - 1):] + xTCRaw[:int(elementsCountAroundTC * 0.5 - 1)] + d1TC += d1TCRaw[int(elementsCountAroundTC * 0.5 - 1):] + d1TCRaw[:int(elementsCountAroundTC * 0.5 - 1)] + + d2TC = [d2 for n in range(len(xTC))] + + xColon, d1Colon, d2Colon, _ = combineTeniaColiWithColon(xColon, d1Colon, d2Colon, [], xTC, d1TC, d2TC, [], + (elementsCountAroundTC - 1) * tcCount, + elementsCountAround, elementsCountAlong, + elementsCountThroughWall, closedProximalEnd) + + return xColon, d1Colon, d2Colon + def createNodesAndElementsTeniaColi(region, x, d1, d2, d3, xFlat, d1Flat, d2Flat, + xOrgan, d1Organ, d2Organ, organCoordinateFieldName, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, @@ -1533,6 +1609,8 @@ def createNodesAndElementsTeniaColi(region, :param x, d1, d2, d3: coordinates and derivatives of coordinates field. :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of flat coordinates field. + :param xOrgan, d1Organ, d2Organ, d3Organ, organCoordinateFieldName: coordinates, + derivatives and name of organ coordinates field. :param elementsCountAroundTC: Number of elements around tenia coli. :param elementsCountAroundHaustrum: Number of elements around haustrum. :param elementsCountAlong: Number of elements along colon. @@ -1644,6 +1722,35 @@ def createNodesAndElementsTeniaColi(region, flatElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) flatElementtemplate5.defineField(flatCoordinates, -1, eftFlat7) + if xOrgan: + # Organ coordinates field + bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) + eftOrgan = bicubichermitelinear.createEftBasic() + + organCoordinates = findOrCreateFieldCoordinates(fm, name=organCoordinateFieldName) + organNodetemplate = nodes.createNodetemplate() + organNodetemplate.defineField(organCoordinates) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + + organElementtemplate = mesh.createElementtemplate() + organElementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + organElementtemplate.defineField(organCoordinates, -1, eftOrgan) + + # Tenia coli edge elements + organElementtemplate1 = mesh.createElementtemplate() + organElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + eftOrgan1 = bicubichermitelinear.createEftWedgeXi1One() + organElementtemplate1.defineField(organCoordinates, -1, eftOrgan1) + + organElementtemplate2 = mesh.createElementtemplate() + organElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + eftOrgan2 = bicubichermitelinear.createEftWedgeXi1Zero() + organElementtemplate2.defineField(organCoordinates, -1, eftOrgan2) + # create nodes for coordinates field for n in range(len(x)): node = nodes.createNode(nodeIdentifier, nodetemplate) @@ -1705,6 +1812,20 @@ def createNodesAndElementsTeniaColi(region, flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 + # Create nodes for organ coordinates field + if xOrgan: + nodeIdentifier = firstNodeIdentifier + for n in range(len(xOrgan)): + node = nodes.findNodeByIdentifier(nodeIdentifier) + node.merge(organNodetemplate) + cache.setNode(node) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xOrgan[n]) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Organ[n]) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Organ[n]) + if useCrossDerivatives: + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + # Create elements allAnnotationGroups = [] @@ -1921,6 +2042,9 @@ def createNodesAndElementsTeniaColi(region, if xFlat: element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) element.setNodesByIdentifier(eftFlat4 if onOpening else eftFlat3, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ annotationGroupsThroughWall[e3] @@ -1956,6 +2080,9 @@ def createNodesAndElementsTeniaColi(region, if xFlat: element.merge(flatElementtemplate1 if eTC < int(elementsCountAroundTC*0.5) - 1 else flatElementtemplate3) element.setNodesByIdentifier(eftFlat3 if eTC < int(elementsCountAroundTC*0.5) - 1 else eftFlat5, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate if eTC < int(elementsCountAroundTC*0.5) - 1 else organElementtemplate1) + element.setNodesByIdentifier(eftOrgan if eTC < int(elementsCountAroundTC*0.5) - 1 else eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] if annotationGroups: @@ -1998,6 +2125,9 @@ def createNodesAndElementsTeniaColi(region, if xFlat: element.merge(flatElementtemplate4) element.setNodesByIdentifier(eftFlat6, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate2) + element.setNodesByIdentifier(eftOrgan2, nodeIdentifiers) elif eTC > 0 and eTC < elementsCountAroundTC - 1: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] @@ -2006,6 +2136,9 @@ def createNodesAndElementsTeniaColi(region, if xFlat: element.merge(flatElementtemplate1) element.setNodesByIdentifier(eftFlat3, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, nodeIdentifiers) else: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni31 + now + tcOffset] @@ -2014,6 +2147,9 @@ def createNodesAndElementsTeniaColi(region, if xFlat: element.merge(flatElementtemplate3) element.setNodesByIdentifier(eftFlat5, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate1) + element.setNodesByIdentifier(eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] @@ -2060,6 +2196,9 @@ def createNodesAndElementsTeniaColi(region, else: element.merge(flatElementtemplate5 if onOpening else flatElementtemplate4) element.setNodesByIdentifier(eftFlat7 if onOpening else eftFlat6, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate if eTC > 0 else organElementtemplate2) + element.setNodesByIdentifier(eftOrgan if eTC > 0 else eftOrgan2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[ elementsCountAround + int(elementsCountAroundTC * (tcCount - 0.5)) + eTC] + \ diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 2b2d1f73..fd56b2de 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -921,14 +921,16 @@ def generateBaseMesh(cls, region, options): flatWidthList, xiList = smallIntestineSegmentTubeMeshInnerPoints.getFlatWidthAndXiList() - # Create flat and texture coordinates + # Create flat coordinates xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( xiList, flatWidthList, length, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) + xOrgan = d1Organ = d2Organ = [] + # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xOrgan, d1Organ, d2Organ, None, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 564ad4d0..165dfc5f 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -398,9 +398,57 @@ def createFlatCoordinates(xiList, lengthAroundList, return xFlatList, d1FlatList, d2FlatList +def createOrganCoordinates(xiList, lengthToDiameterRatio, wallThicknessToDiameterRatio, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): + """ + Calculates organ coordinates and derivatives represented by a cylindrical tube with a unit inner diameter, + length equivalent to lengthToDiameterRatio and wall thickness of wallThicknessToDiameterRatio. + :param xiList: List containing xi for each point around the outer surface of the tube. + :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ + :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. + :param elementsCountAround: Number of elements around tube. + :param elementsCountAlong: Number of elements along tube. + :param elementsCountThroughWall: Number of elements through wall. + :param transitElementList: stores true if element around is a transition element between a big and small element. + :return: coordinates and derivatives of organ coordinates field. + """ + # Calculate organ coordinates and derivatives + xOrganList = [] + d1OrganList = [] + d2OrganList = [] + d2 = [0.0, lengthToDiameterRatio / elementsCountAlong, 0.0] + + for n2 in range(elementsCountAlong + 1): + cx = [0.0, lengthToDiameterRatio / elementsCountAlong * n2, 0.0] + xiFace = xiList[n2] + for n3 in range(elementsCountThroughWall + 1): + radius = 0.5 + wallThicknessToDiameterRatio/elementsCountThroughWall * n3 + axis1 = [0.0, 0.0, radius] + axis2 = [radius, 0.0, 0.0] + d1List = [] + for n1 in range(len(xiFace) - 1): + dTheta = (xiFace[n1 + 1 if n1 < len(xiFace) - 1 else 0] - xiFace[n1]) * 2.0 * math.pi + radiansAround = 2.0 * math.pi * xiFace[n1] + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + xOrganList.append([(cx[c] + cosRadiansAround * axis1[c] + sinRadiansAround * axis2[c]) for c in range(3)]) + d1List.append([ dTheta*(-sinRadiansAround*axis1[c] + cosRadiansAround*axis2[c]) for c in range(3) ]) + d2OrganList.append(d2) + + # To modify derivative along transition elements + for i in range(len(transitElementList)): + if transitElementList[i]: + d1List[i] = vector.setMagnitude(d1List[i], vector.magnitude(d1List[i - 1])) + d1List[i + 1] = vector.setMagnitude(d1List[i+ 1], vector.magnitude(d1List[(i + 2) % elementsCountAround])) + + d1OrganList += d1List + + return xOrganList, d1OrganList, d2OrganList + def createNodesAndElements(region, x, d1, d2, d3, xFlat, d1Flat, d2Flat, + xOrgan, d1Organ, d2Organ, organCoordinateFieldName, elementsCountAround, elementsCountAlong, elementsCountThroughWall, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, @@ -410,6 +458,8 @@ def createNodesAndElements(region, :param x, d1, d2, d3: coordinates and derivatives of coordinates field. :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of flat coordinates field. + :param xOrgan, d1Organ, d2Organ, d3Organ, organCoordinateFieldName: coordinates, derivatives and name of organ + coordinates field. :param elementsCountAround: Number of elements around tube. :param elementsCountAlong: Number of elements along tube. :param elementsCountThroughWall: Number of elements through wall. @@ -491,6 +541,24 @@ def createNodesAndElements(region, flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) flatElementtemplate2.defineField(flatCoordinates, -1, eftFlat2) + if xOrgan: + # Organ coordinates field + bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) + eftOrgan = bicubichermitelinear.createEftBasic() + + organCoordinates = findOrCreateFieldCoordinates(fm, name=organCoordinateFieldName) + organNodetemplate = nodes.createNodetemplate() + organNodetemplate.defineField(organCoordinates) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + + organElementtemplate = mesh.createElementtemplate() + organElementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + organElementtemplate.defineField(organCoordinates, -1, eftOrgan) + # Create nodes # Coordinates field for n in range(len(x)): @@ -533,6 +601,20 @@ def createNodesAndElements(region, flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 + # Organ coordinates field + if xOrgan: + nodeIdentifier = firstNodeIdentifier + for n in range(len(xOrgan)): + node = nodes.findNodeByIdentifier(nodeIdentifier) + node.merge(organNodetemplate) + cache.setNode(node) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xOrgan[n]) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Organ[n]) + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Organ[n]) + if useCrossDerivatives: + organCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + # create elements elementtemplate3 = mesh.createElementtemplate() elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) @@ -598,6 +680,9 @@ def createNodesAndElements(region, if xFlat: element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) element.setNodesByIdentifier(eftFlat2 if onOpening else eftFlat1, nodeIdentifiers) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ From 83258f382cea67b7a058aa97e65fb171d0ca01df Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 12 Oct 2021 18:18:47 +1300 Subject: [PATCH 04/21] Update tenia coli width for mouse colon --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 5854e4ce..3834b02a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -313,11 +313,11 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Proximal inner radius'] = 1.0 options['Proximal tenia coli width'] = 0.5 options['Proximal-transverse inner radius'] = 0.9 - options['Proximal-transverse tenia coli width'] = 0.7 + options['Proximal-transverse tenia coli width'] = 0.5 options['Transverse-distal inner radius'] = 0.7 - options['Transverse-distal tenia coli width'] = 1.0 + options['Transverse-distal tenia coli width'] = 0.5 options['Distal inner radius'] = 0.7 - options['Distal tenia coli width'] = 1.0 + options['Distal tenia coli width'] = 0.5 elif 'Pig 1' in parameterSetName: options['Number of segments'] = 120 options['Proximal length'] = 3000.0 From 0d97f4bf94d851339b8af2ceef1366579c31176a Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 13 Oct 2021 11:42:22 +1300 Subject: [PATCH 05/21] Allow varying thickness through wall elements --- .../meshtypes/meshtype_3d_bladderurethra1.py | 2 + .../meshtypes/meshtype_3d_cecum1.py | 3 +- .../meshtypes/meshtype_3d_colon1.py | 19 +++++-- .../meshtypes/meshtype_3d_colonsegment1.py | 50 +++++++++++++++---- .../meshtypes/meshtype_3d_smallintestine1.py | 5 +- src/scaffoldmaker/utils/tubemesh.py | 49 +++++++++++++----- 6 files changed, 97 insertions(+), 31 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 2021ce78..cfe8f1b1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -860,9 +860,11 @@ def generateBaseMesh(cls, region, options): wallThicknessList = [bladderWallThickness] * (elementsCountAlong + 1) transitElementList = [0] * elementsCountAround + relativeThicknessList = [] xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList, wallThicknessList, + relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index f411c032..18d19591 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -430,8 +430,9 @@ def generateBaseMesh(cls, region, options): # Create coordinates and derivatives wallThicknessList = [wallThickness] * (elementsCountAlong + 1) + relativeThicknessList = [] xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, - d2WarpedList, d3WarpedUnitList, wallThicknessList, + d2WarpedList, d3WarpedUnitList, wallThicknessList, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) # Deal with multiple nodes at end point for closed proximal end diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 3834b02a..c7d07294 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -471,6 +471,10 @@ def generateBaseMesh(cls, region, options): elementsCountAlongSegment = segmentSettings['Number of elements along segment'] elementsCountThroughWall = segmentSettings['Number of elements through wall'] wallThickness = segmentSettings['Wall thickness'] + mucosaRelThickness = segmentSettings['Mucosa relative thickness'] + submucosaRelThickness = segmentSettings['Submucosa relative thickness'] + circularRelThickness = segmentSettings['Circular muscle layer relative thickness'] + longitudinalRelThickness = segmentSettings['Longitudinal muscle layer relative thickness'] useCrossDerivatives = segmentSettings['Use cross derivatives'] useCubicHermiteThroughWall = not(segmentSettings['Use linear through wall']) elementsCountAlong = int(elementsCountAlongSegment*segmentCount) @@ -580,6 +584,9 @@ def generateBaseMesh(cls, region, options): d3UnitExtrude = [] sxRefExtrudeList = [] + relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, circularRelThickness, + longitudinalRelThickness] + # Create object colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, @@ -618,7 +625,7 @@ def generateBaseMesh(cls, region, options): # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xExtrude, d1Extrude, - d2Extrude, d3UnitExtrude, contractedWallThicknessList, + d2Extrude, d3UnitExtrude, contractedWallThicknessList, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() @@ -635,12 +642,13 @@ def generateBaseMesh(cls, region, options): # Create flat coordinates xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( - xiList, relaxedLengthList, length, wallThickness, tcCount, tcThickness, + xiList, relaxedLengthList, length, wallThickness, relativeThicknessList, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd) # Create colon coordinates - xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, + xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, relativeThicknessList, + lengthToDiameterRatio, wallThicknessToDiameterRatio, teniaColiThicknessToDiameterRatio, tcCount, elementsCountAroundTC, @@ -659,11 +667,12 @@ def generateBaseMesh(cls, region, options): else: # Create flat coordinates xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( - xiList, relaxedLengthList, length, wallThickness, elementsCountAround, + xiList, relaxedLengthList, length, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) # Create colon coordinates - xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, lengthToDiameterRatio, + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessList, + lengthToDiameterRatio, wallThicknessToDiameterRatio, elementsCountAround, elementsCountAlong, elementsCountThroughWall, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 4e474e58..59eb5a7b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -54,7 +54,7 @@ def getDefaultOptions(parameterSetName='Default'): 'Number of elements around tenia coli' : 2, 'Number of elements around haustrum' : 8, 'Number of elements along segment' : 4, - 'Number of elements through wall' : 1, + 'Number of elements through wall' : 4, 'Start phase': 0.0, 'Start inner radius': 43.5, 'Start inner radius derivative': 0.0, @@ -72,6 +72,10 @@ def getDefaultOptions(parameterSetName='Default'): 'End tenia coli width derivative': 0.0, 'Tenia coli thickness': 1.6, 'Wall thickness': 1.6, + 'Mucosa relative thickness': 0.3, + 'Submucosa relative thickness': 0.3, + 'Circular muscle layer relative thickness': 0.3, + 'Longitudinal muscle layer relative thickness': 0.1, 'Use cross derivatives' : False, 'Use linear through wall' : True, 'Refine' : False, @@ -103,6 +107,11 @@ def getDefaultOptions(parameterSetName='Default'): options['End tenia coli width'] = 0.8 options['Tenia coli thickness'] = 0.0 options['Wall thickness'] = 0.55 + options['Mucosa relative thickness'] = 0.4 + options['Submucosa relative thickness'] = 0.1 + options['Circular muscle layer relative thickness'] = 0.3 + options['Longitudinal muscle layer relative thickness'] = 0.2 + elif 'Pig' in parameterSetName: options['Start inner radius'] = 20.0 options['End inner radius'] = 20.0 @@ -116,6 +125,11 @@ def getDefaultOptions(parameterSetName='Default'): options['End tenia coli width'] = 5.0 options['Tenia coli thickness'] = 0.5 options['Wall thickness'] = 2.0 + options['Mucosa relative thickness'] = 0.34 + options['Submucosa relative thickness'] = 0.25 + options['Circular muscle layer relative thickness'] = 0.25 + options['Longitudinal muscle layer relative thickness'] = 0.16 + return options @staticmethod @@ -142,6 +156,10 @@ def getOrderedOptionNames(): 'End tenia coli width derivative', 'Tenia coli thickness', 'Wall thickness', + 'Mucosa relative thickness', + 'Submucosa relative thickness', + 'Circular muscle layer relative thickness', + 'Longitudinal muscle layer relative thickness', 'Use cross derivatives', 'Use linear through wall', 'Refine', @@ -179,7 +197,11 @@ def checkOptions(options): 'Segment length mid derivative factor', 'Segment length', 'Tenia coli thickness', - 'Wall thickness']: + 'Wall thickness', + 'Mucosa relative thickness', + 'Submucosa relative thickness', + 'Circular muscle layer relative thickness', + 'Longitudinal muscle layer relative thickness' ]: if options[key] < 0.0: options[key] = 0.0 if options['Corner inner radius factor'] < 0.1: @@ -232,6 +254,10 @@ def generateBaseMesh(cls, region, options): endTCWidthDerivative = options['End tenia coli width derivative'] tcThickness = options['Tenia coli thickness'] wallThickness = options['Wall thickness'] + mucosaRelThickness = options['Mucosa relative thickness'] + submucosaRelThickness = options['Submucosa relative thickness'] + circularRelThickness = options['Circular muscle layer relative thickness'] + longitudinalRelThickness = options['Longitudinal muscle layer relative thickness'] useCrossDerivatives = options['Use cross derivatives'] useCubicHermiteThroughWall = not(options['Use linear through wall']) elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount @@ -306,10 +332,12 @@ def generateBaseMesh(cls, region, options): closedProximalEnd) contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() + relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, + circularRelThickness, longitudinalRelThickness] # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, - d2WarpedList, d3WarpedUnitList, contractedWallThicknessList, + d2WarpedList, d3WarpedUnitList, contractedWallThicknessList, relativeThicknessList, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) xColonSegment = d1ColonSegment = d2ColonSegment = [] @@ -325,7 +353,7 @@ def generateBaseMesh(cls, region, options): # Create flat coordinates xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( - xiList, relaxedLengthList, segmentLength, wallThickness, tcCount, tcThickness, + xiList, relaxedLengthList, segmentLength, wallThickness, relativeThicknessList, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, transitElementList, closedProximalEnd) @@ -339,7 +367,7 @@ def generateBaseMesh(cls, region, options): else: # Create flat coordinates xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( - xiList, relaxedLengthList, segmentLength, wallThickness, elementsCountAround, + xiList, relaxedLengthList, segmentLength, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) # Create nodes and elements @@ -1426,7 +1454,7 @@ def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, return x, d1, d2, d3 def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, - totalLengthAlong, wallThickness, tcCount, tcThickness, + totalLengthAlong, wallThickness, relativeThicknessList, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd): """ @@ -1438,6 +1466,7 @@ def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, its most relaxed state for each element along. :param totalLengthAlong: Total length along colon. :param wallThickness: Thickness of wall. + :param relativeThicknessList: Relative thickness of each element through wall. :param tcCount: Number of tenia coli. :param tcThickness: Thickness of tenia coli at its thickest region. :param elementsCountAroundTC: Number of elements around tenia coli. @@ -1455,8 +1484,8 @@ def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount # Find flat coordinates for colon - xFlatColon, d1FlatColon, d2FlatColon = tubemesh.createFlatCoordinates(xiList, relaxedLengthList, - totalLengthAlong, wallThickness, + xFlatColon, d1FlatColon, d2FlatColon = tubemesh.createFlatCoordinates(xiList, relaxedLengthList, totalLengthAlong, + wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) @@ -1521,7 +1550,7 @@ def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, return xFlat, d1Flat, d2Flat -def createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, wallThicknessToDiameterRatio, +def createColonCoordinatesTeniaColi(xiList, relativeThicknessList, lengthToDiameterRatio, wallThicknessToDiameterRatio, teniaColiThicknessToDiameterRatio, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd): @@ -1530,6 +1559,7 @@ def createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, wallThickness with unit inner diameter, length of lengthToDiameterRatio, wall thickness of wallThicknessToDiameterRatio, with tenia coli of teniaColiThicknessToDiameterRatio running along its length. :param xiList: List containing xi for each point around the outer surface of colon in its most relaxed state. + :param relativeThicknessList: Relative thickness for each element through wall. :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. :param teniaColiThicknessToDiameterRatio: Ratio of tenia coli thickness to inner diameter of organ. @@ -1546,7 +1576,7 @@ def createColonCoordinatesTeniaColi(xiList, lengthToDiameterRatio, wallThickness elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount # Find organ coordinates for colon - xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, lengthToDiameterRatio, + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessList, lengthToDiameterRatio, wallThicknessToDiameterRatio, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index fd56b2de..bf1eae7d 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -859,6 +859,7 @@ def generateBaseMesh(cls, region, options): d1Extrude = [] d2Extrude = [] d3UnitExtrude = [] + relativeThicknessList = [] # Create object smallIntestineSegmentTubeMeshInnerPoints = CylindricalSegmentTubeMeshInnerPoints( @@ -916,14 +917,14 @@ def generateBaseMesh(cls, region, options): # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xExtrude, d1Extrude, - d2Extrude, d3UnitExtrude, [wallThickness]*(elementsCountAlong+1), + d2Extrude, d3UnitExtrude, [wallThickness]*(elementsCountAlong+1), relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) flatWidthList, xiList = smallIntestineSegmentTubeMeshInnerPoints.getFlatWidthAndXiList() # Create flat coordinates xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( - xiList, flatWidthList, length, wallThickness, elementsCountAround, + xiList, flatWidthList, length, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) xOrgan = d1Organ = d2Organ = [] diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 165dfc5f..c1550de3 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -238,7 +238,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, return xWarpedList, d1WarpedList, d2WarpedListFinal, d3WarpedUnitList def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, - wallThicknessList, elementsCountAround, + wallThicknessList, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ Generates coordinates from inner to outer surface using coordinates @@ -248,6 +248,7 @@ def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, :param d2Inner: Derivatives on inner surface along tube :param d3Inner: Derivatives on inner surface through wall :param wallThicknessList: Wall thickness for each element along tube + :param relativeThicknessList: Relative wall thickness for each element through wall :param elementsCountAround: Number of elements around tube :param elementsCountAlong: Number of elements along tube :param elementsCountThroughWall: Number of elements through tube wall @@ -265,6 +266,14 @@ def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, d2List = [] d3List = [] + if relativeThicknessList: + xi3 = 0.0 + xi3List = [0.0] + for n3 in range(elementsCountThroughWall): + xi3 += relativeThicknessList[n3] + xi3List.append(xi3) + relativeThicknessList.append(relativeThicknessList[-1]) + for n2 in range(elementsCountAlong + 1): wallThickness = wallThicknessList[n2] for n1 in range(elementsCountAround): @@ -305,7 +314,7 @@ def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, curvatureAlong.append(curvature) for n3 in range(elementsCountThroughWall + 1): - xi3 = 1/elementsCountThroughWall * n3 + xi3 = xi3List[n3] if relativeThicknessList else 1.0/elementsCountThroughWall * n3 for n1 in range(elementsCountAround): n = n2*elementsCountAround + n1 norm = d3Inner[n] @@ -330,20 +339,17 @@ def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, curvatureList.append(curvature) #dx_ds3 - d3 = [c * wallThickness/elementsCountThroughWall for c in norm] + d3 = [c * wallThickness * (relativeThicknessList[n3] if relativeThicknessList else 1.0/elementsCountThroughWall) for c in norm] d3List.append(d3) return xList, d1List, d2List, d3List, curvatureList -def createFlatCoordinates(xiList, lengthAroundList, - totalLengthAlong, wallThickness, elementsCountAround, - elementsCountAlong, elementsCountThroughWall, transitElementList): +def createFlatCoordinates(xiList, lengthAroundList, totalLengthAlong, wallThickness, relativeThicknessList, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ Calculates flat coordinates for a tube when it is opened into a flat preparation. - :param xiList: List containing xi for each point around - the outer surface of the tube. - :param lengthAroundList: List of total arclength around the - outer surface for each element along. + :param xiList: List containing xi for each point around the outer surface of the tube. + :param lengthAroundList: List of total arclength around the outer surface for each element along. :param totalLengthAlong: Total length along tube. :param wallThickness: Thickness of wall. :param elementsCountAround: Number of elements around tube. @@ -353,6 +359,13 @@ def createFlatCoordinates(xiList, lengthAroundList, transition element between a big and small element. :return: coordinates and derivatives of flat coordinates field. """ + if relativeThicknessList: + xi3 = 0.0 + xi3List = [0.0] + for n3 in range(elementsCountThroughWall): + xi3 += relativeThicknessList[n3] + xi3List.append(xi3) + relativeThicknessList.append(relativeThicknessList[-1]) # Calculate flat coordinates and derivatives xFlatList = [] @@ -373,7 +386,7 @@ def createFlatCoordinates(xiList, lengthAroundList, xPad = (lengthAroundList[0] - lengthAround)*0.5 for n3 in range(elementsCountThroughWall + 1): - z = wallThickness / elementsCountThroughWall * n3 + z = wallThickness * (xi3List[n3] if relativeThicknessList else 1.0 / elementsCountThroughWall * n3) for n1 in range(elementsCountAround + 1): xFlat = [xPad + xiFace[n1] * lengthAround, totalLengthAlong / elementsCountAlong * n2, @@ -398,12 +411,13 @@ def createFlatCoordinates(xiList, lengthAroundList, return xFlatList, d1FlatList, d2FlatList -def createOrganCoordinates(xiList, lengthToDiameterRatio, wallThicknessToDiameterRatio, +def createOrganCoordinates(xiList, relativeThicknessList, lengthToDiameterRatio, wallThicknessToDiameterRatio, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ Calculates organ coordinates and derivatives represented by a cylindrical tube with a unit inner diameter, length equivalent to lengthToDiameterRatio and wall thickness of wallThicknessToDiameterRatio. :param xiList: List containing xi for each point around the outer surface of the tube. + :param relativeThicknessList: Relative thickness of each element through wall. :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. :param elementsCountAround: Number of elements around tube. @@ -412,6 +426,14 @@ def createOrganCoordinates(xiList, lengthToDiameterRatio, wallThicknessToDiamete :param transitElementList: stores true if element around is a transition element between a big and small element. :return: coordinates and derivatives of organ coordinates field. """ + if relativeThicknessList: + xi3 = 0.0 + xi3List = [0.0] + for n3 in range(elementsCountThroughWall): + xi3 += relativeThicknessList[n3] + xi3List.append(xi3) + relativeThicknessList.append(relativeThicknessList[-1]) + # Calculate organ coordinates and derivatives xOrganList = [] d1OrganList = [] @@ -422,7 +444,8 @@ def createOrganCoordinates(xiList, lengthToDiameterRatio, wallThicknessToDiamete cx = [0.0, lengthToDiameterRatio / elementsCountAlong * n2, 0.0] xiFace = xiList[n2] for n3 in range(elementsCountThroughWall + 1): - radius = 0.5 + wallThicknessToDiameterRatio/elementsCountThroughWall * n3 + xi3 = xi3List[n3] if relativeThicknessList else 1.0/elementsCountThroughWall * n3 + radius = 0.5 + wallThicknessToDiameterRatio * xi3 axis1 = [0.0, 0.0, radius] axis2 = [radius, 0.0, 0.0] d1List = [] From e235d770c045e599a8e3168cbf75e7e6a795758e Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 13 Oct 2021 11:50:40 +1300 Subject: [PATCH 06/21] Restrict number of elements through wall to 4 in colon --- src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 59eb5a7b..1c86130b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -171,7 +171,6 @@ def getOrderedOptionNames(): @staticmethod def checkOptions(options): for key in [ - 'Number of elements through wall', 'Refine number of elements around', 'Refine number of elements along segment', 'Refine number of elements through wall']: @@ -184,6 +183,8 @@ def checkOptions(options): options[key] = 2 if options['Number of elements around haustrum'] < 4: options['Number of elements around haustrum'] = 4 + if options['Number of elements through wall'] != 4: + options['Number of elements through wall'] = 4 for key in [ 'Number of elements around tenia coli', 'Number of elements around haustrum']: From a1923f99eed511b259ae543e5582774623823922 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 13 Oct 2021 12:42:59 +1300 Subject: [PATCH 07/21] Annotate colon wall layers --- src/scaffoldmaker/annotation/colon_terms.py | 3 + .../meshtypes/meshtype_3d_colon1.py | 43 ++++++++++++-- .../meshtypes/meshtype_3d_colonsegment1.py | 56 ++++++++++++++----- 3 files changed, 85 insertions(+), 17 deletions(-) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index 0c4fc010..cc27bccf 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -7,10 +7,13 @@ ( "ascending colon", "UBERON:0001156", "FMA:14545", "ILX:0734393"), ( "descending colon", "UBERON:0001158", "FMA:14547", "ILX:0724444"), ( "caecum", "UBERON:0001153", "FMA:14541", "ILX:0732270"), + ( "Circular muscle layer of colon", "ILX:0772428"), ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"), ( "colonic mucosa", "UBERON:0000317", "FMA:14984", "ILX:0731046"), ( "distal colon", "UBERON:0008971", "ILX:0727523"), + ( "Longitudinal muscle layer of colon", "ILX:0775554"), ( "mesenteric zone", "None"), + ( "myenteric nerve plexus", "UBERON:0002439", "FMA:63252", "ILX:0725342"), ( "non-mesenteric zone", "None"), ( "proximal colon", "UBERON:0008972", "ILX:0733240"), ( "serosa of colon", "UBERON:0003335", "FMA:14990", "ILX:0736932"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index c7d07294..31cf2110 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -5,7 +5,8 @@ """ import copy -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from opencmiss.zinc.element import Element +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ @@ -574,9 +575,11 @@ def generateBaseMesh(cls, region, options): for n in range(elementsCount): annotationGroupsAlong.append(annotationGroupAlong[i]) - annotationGroupsThroughWall = [] - for i in range(elementsCountThroughWall): - annotationGroupsThroughWall.append([ ]) + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], [circularMuscleGroup], [longitudinalMuscleGroup]] xExtrude = [] d1Extrude = [] @@ -702,3 +705,35 @@ def refineMesh(cls, meshrefinement, options): meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) return + + @classmethod + def defineFaceAnnotations(cls, region, options, annotationGroups): + ''' + 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 annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. + ''' + # Create 2d surface mesh groups + fm = region.getFieldmodule() + longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + get_colon_term("Longitudinal muscle layer of colon")) + circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + get_colon_term("Circular muscle layer of colon")) + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + + is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) + is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) + is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) + + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + + myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, + get_colon_term("myenteric nerve plexus")) + myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 1c86130b..c1a19845 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -7,11 +7,11 @@ import copy import math -from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldTextureCoordinates +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear @@ -305,14 +305,11 @@ def generateBaseMesh(cls, region, options): for i in range(elementsCountAlongSegment): annotationGroupsAlong.append([ ]) - # serosaGroup = AnnotationGroup(region, get_colon_term("serosa of colon")) - # submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) - # mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) - # annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], [serosaGroup]] - - annotationGroupsThroughWall = [] - for i in range(elementsCountThroughWall): - annotationGroupsThroughWall.append([ ]) + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], [circularMuscleGroup], [longitudinalMuscleGroup]] # Create inner points nSegment = 0 @@ -396,6 +393,37 @@ def refineMesh(cls, meshrefinement, options): refineElementsCountThroughWall) return + @classmethod + def defineFaceAnnotations(cls, region, options, annotationGroups): + ''' + 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 annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. + ''' + # Create 2d surface mesh groups + fm = region.getFieldmodule() + longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + get_colon_term("Longitudinal muscle layer of colon")) + circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + get_colon_term("Circular muscle layer of colon")) + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + + is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) + is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) + is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) + + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + + myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) + myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) + class ColonSegmentTubeMeshInnerPoints: """ Generates inner profile of a colon segment for use by tubemesh. @@ -2115,7 +2143,8 @@ def createNodesAndElementsTeniaColi(region, element.merge(organElementtemplate if eTC < int(elementsCountAroundTC*0.5) - 1 else organElementtemplate1) element.setNodesByIdentifier(eftOrgan if eTC < int(elementsCountAroundTC*0.5) - 1 else eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 - annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] + annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] + \ + annotationGroupsThroughWall[-1] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2183,7 +2212,8 @@ def createNodesAndElementsTeniaColi(region, element.setNodesByIdentifier(eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( - elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] + elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] + \ + annotationGroupsThroughWall[-1] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2233,7 +2263,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[ elementsCountAround + int(elementsCountAroundTC * (tcCount - 0.5)) + eTC] + \ - annotationGroupsAlong[e2] + annotationGroupsAlong[e2] + annotationGroupsThroughWall[-1] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: From a5bb3ceef81420296ee19562b8c460acf1c78d26 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 13 Oct 2021 13:11:33 +1300 Subject: [PATCH 08/21] Update unit tests --- tests/test_colon.py | 46 ++++++++++++++++++------------------ tests/test_colonsegment.py | 36 ++++++++-------------------- tests/test_smallintestine.py | 16 +------------ 3 files changed, 34 insertions(+), 64 deletions(-) diff --git a/tests/test_colon.py b/tests/test_colon.py index cafa05f7..5020d32c 100644 --- a/tests/test_colon.py +++ b/tests/test_colon.py @@ -91,7 +91,7 @@ def test_colon1(self): del tmpRegion annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options) - self.assertEqual(7, len(annotationGroups)) + self.assertEqual(11, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) @@ -99,13 +99,13 @@ def test_colon1(self): for annotationGroup in annotationGroups: annotationGroup.addSubelements() mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(432, mesh3d.getSize()) + self.assertEqual(1512, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(1656, mesh2d.getSize()) + self.assertEqual(4986, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(2043, mesh1d.getSize()) + self.assertEqual(5463, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(819, nodes.getSize()) + self.assertEqual(1989, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -121,10 +121,10 @@ def test_colon1(self): assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) assertAlmostEqualList(self, maximums, [ 186.72988844629867, 77.41781871321301, 3.2000000000000006 ], 1.0E-6) - textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() - minimums, maximums = evaluateFieldNodesetRange(textureCoordinates, nodes) - assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) - assertAlmostEqualList(self, maximums, [ 0.9812471574796385, 1.0, 2.0 ], 1.0E-6) + colonCoordinates = fieldmodule.findFieldByName("colon coordinates").castFiniteElement() + minimums, maximums = evaluateFieldNodesetRange(colonCoordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.6, 0.0, -0.6], 1.0E-4) + assertAlmostEqualList(self, maximums, [ 0.6, 24.0, 0.605], 1.0E-4) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -139,7 +139,7 @@ def test_colon1(self): self.assertAlmostEqual(surfaceArea, 14612.416789097502, delta=1.0E-6) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 26826.06954301569, delta=1.0E-6) + self.assertAlmostEqual(volume, 26825.42839677291, delta=1.0E-6) def test_mousecolon1(self): """ @@ -150,38 +150,38 @@ def test_mousecolon1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options) - self.assertEqual(6, len(annotationGroups)) + self.assertEqual(10, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(400, mesh3d.getSize()) + self.assertEqual(1600, mesh3d.getSize()) coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() self.assertTrue(coordinates.isValid()) flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement() self.assertTrue(flatCoordinates.isValid()) - textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() - self.assertTrue(textureCoordinates.isValid()) + colonCoordinates = fieldmodule.findFieldByName("colon coordinates").castFiniteElement() + self.assertTrue(colonCoordinates.isValid()) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1) flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup) flatSurfaceAreaField.setNumbersOfPoints(4) - textureSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, faceMeshGroup) - textureSurfaceAreaField.setNumbersOfPoints(4) - textureVolumeField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, mesh3d) - textureVolumeField.setNumbersOfPoints(3) + colonSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, colonCoordinates, faceMeshGroup) + colonSurfaceAreaField.setNumbersOfPoints(4) + colonVolumeField = fieldmodule.createFieldMeshIntegral(one, colonCoordinates, mesh3d) + colonVolumeField.setNumbersOfPoints(3) fieldcache = fieldmodule.createFieldcache() result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(flatSurfaceArea, 629.440514706522, delta=1.0E-6) - result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1) + self.assertAlmostEqual(flatSurfaceArea, 629.4883774904393, delta=1.0E-6) + result, colonSurfaceArea = colonSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6) - result, textureVolume = textureVolumeField.evaluateReal(fieldcache, 1) + self.assertAlmostEqual(colonSurfaceArea, 90.4578820802557, delta=1.0E-6) + result, colonVolume = colonVolumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureVolume, 1.0, delta=1.0E-6) + self.assertAlmostEqual(colonVolume, 8.290058800222006, delta=1.0E-6) if __name__ == "__main__": unittest.main() diff --git a/tests/test_colonsegment.py b/tests/test_colonsegment.py index 95cbb77a..725c0772 100644 --- a/tests/test_colonsegment.py +++ b/tests/test_colonsegment.py @@ -18,11 +18,11 @@ def test_humancolonsegment1(self): parameterSetNames = MeshType_3d_colonsegment1.getParameterSetNames() self.assertEqual(parameterSetNames, [ "Default", "Cattle 1", "Human 1", "Mouse 1", "Pig 1" ]) options = MeshType_3d_colonsegment1.getDefaultOptions("Human 1") - self.assertEqual(27, len(options)) + self.assertEqual(31, len(options)) self.assertEqual(0.0, options.get("Start phase")) self.assertEqual(2, options.get("Number of elements around tenia coli")) self.assertEqual(4, options.get("Number of elements along segment")) - self.assertEqual(1, options.get("Number of elements through wall")) + self.assertEqual(4, options.get("Number of elements through wall")) self.assertEqual(43.5, options.get("Start inner radius")) self.assertEqual(0.0, options.get("Start inner radius derivative")) self.assertEqual(33.0, options.get("End inner radius")) @@ -38,18 +38,18 @@ def test_humancolonsegment1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options) - self.assertEqual(3, len(annotationGroups)) + self.assertEqual(7, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(144, mesh3d.getSize()) + self.assertEqual(504, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(576, mesh2d.getSize()) + self.assertEqual(1746, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(747, mesh1d.getSize()) + self.assertEqual(2007, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(315, nodes.getSize()) + self.assertEqual(765, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -65,11 +65,6 @@ def test_humancolonsegment1(self): assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) assertAlmostEqualList(self, maximums, [397.2736607240895, 50.0, 3.2000000000000006], 1.0E-6) - textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() - minimums, maximums = evaluateFieldNodesetRange(textureCoordinates, nodes) - assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) - assertAlmostEqualList(self, maximums, [ 0.9887737064410717, 1.0, 2.0 ], 1.0E-6) - with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1) @@ -94,18 +89,16 @@ def test_mousecolonsegment1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options) - self.assertEqual(2, len(annotationGroups)) + self.assertEqual(6, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(40, mesh3d.getSize()) + self.assertEqual(160, mesh3d.getSize()) coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() self.assertTrue(coordinates.isValid()) flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement() self.assertTrue(flatCoordinates.isValid()) - textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() - self.assertTrue(textureCoordinates.isValid()) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -114,10 +107,7 @@ def test_mousecolonsegment1(self): surfaceAreaField.setNumbersOfPoints(4) flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup) flatSurfaceAreaField.setNumbersOfPoints(4) - textureSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, faceMeshGroup) - textureSurfaceAreaField.setNumbersOfPoints(4) - textureVolumeField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, mesh3d) - textureVolumeField.setNumbersOfPoints(3) + fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) @@ -125,12 +115,6 @@ def test_mousecolonsegment1(self): result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) self.assertAlmostEqual(flatSurfaceArea, surfaceArea, delta=1.0E-3) - result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1) - self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6) - result, textureVolume = textureVolumeField.evaluateReal(fieldcache, 1) - self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureVolume, 1.0, delta=1.0E-6) if __name__ == "__main__": unittest.main() diff --git a/tests/test_smallintestine.py b/tests/test_smallintestine.py index dad7a4d1..425d6418 100644 --- a/tests/test_smallintestine.py +++ b/tests/test_smallintestine.py @@ -97,11 +97,6 @@ def test_smallintestine1(self): assertAlmostEqualList(self, minimums, [ -1.39038154442654, 0.0, 0.0 ], 1.0E-6) assertAlmostEqualList(self, maximums, [ 4.891237158967401, 25.293706698841913, 0.1 ], 1.0E-6) - textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() - minimums, maximums = evaluateFieldNodesetRange(textureCoordinates, nodes) - assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) - assertAlmostEqualList(self, maximums, [ 0.875, 1.0, 1.0 ], 1.0E-6) - with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1) @@ -111,10 +106,7 @@ def test_smallintestine1(self): volumeField.setNumbersOfPoints(3) flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup) flatSurfaceAreaField.setNumbersOfPoints(4) - textureSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, faceMeshGroup) - textureSurfaceAreaField.setNumbersOfPoints(4) - textureVolumeField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, mesh3d) - textureVolumeField.setNumbersOfPoints(3) + fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) @@ -125,12 +117,6 @@ def test_smallintestine1(self): result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) self.assertAlmostEqual(flatSurfaceArea, 171.37026123844635, delta=1.0E-3) - result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1) - self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6) - result, textureVolume = textureVolumeField.evaluateReal(fieldcache, 1) - self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(textureVolume, 1.0, delta=1.0E-6) if __name__ == "__main__": unittest.main() From 15731d7dcba03dcfa9dc5b4fbfc53283597cc058 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 12:34:39 +1300 Subject: [PATCH 09/21] Match tenai coli thickness to 0.25 of wall thickness --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 31cf2110..93254926 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -483,7 +483,7 @@ def generateBaseMesh(cls, region, options): # Colon coordinates lengthToDiameterRatio = 24 wallThicknessToDiameterRatio = 0.1 - teniaColiThicknessToDiameterRatio = 0.005 + teniaColiThicknessToDiameterRatio = 0.025 firstNodeIdentifier = 1 firstElementIdentifier = 1 From 5569558edae798d7569a0d17f60c6a036cac85d7 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 12:50:17 +1300 Subject: [PATCH 10/21] Use 0.25 as thickness ratio for all layers for colon coordinates --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 5 +++-- src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py | 2 +- src/scaffoldmaker/utils/tubemesh.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 93254926..1c66bbf7 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -484,6 +484,7 @@ def generateBaseMesh(cls, region, options): lengthToDiameterRatio = 24 wallThicknessToDiameterRatio = 0.1 teniaColiThicknessToDiameterRatio = 0.025 + relativeThicknessListColonCoordinates = [0.25, 0.25, 0.25, 0.25] firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -650,7 +651,7 @@ def generateBaseMesh(cls, region, options): elementsCountThroughWall, transitElementList, closedProximalEnd) # Create colon coordinates - xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, relativeThicknessList, + xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, relativeThicknessListColonCoordinates, lengthToDiameterRatio, wallThicknessToDiameterRatio, teniaColiThicknessToDiameterRatio, tcCount, @@ -674,7 +675,7 @@ def generateBaseMesh(cls, region, options): elementsCountAlong, elementsCountThroughWall, transitElementList) # Create colon coordinates - xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessList, + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessListColonCoordinates, lengthToDiameterRatio, wallThicknessToDiameterRatio, elementsCountAround, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index c1a19845..b1d63f00 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -1588,7 +1588,7 @@ def createColonCoordinatesTeniaColi(xiList, relativeThicknessList, lengthToDiame with unit inner diameter, length of lengthToDiameterRatio, wall thickness of wallThicknessToDiameterRatio, with tenia coli of teniaColiThicknessToDiameterRatio running along its length. :param xiList: List containing xi for each point around the outer surface of colon in its most relaxed state. - :param relativeThicknessList: Relative thickness for each element through wall. + :param relativeThicknessList: Relative thickness for each element through wall for colon coordinates. :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. :param teniaColiThicknessToDiameterRatio: Ratio of tenia coli thickness to inner diameter of organ. diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index c1550de3..8a3caf1f 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -417,7 +417,7 @@ def createOrganCoordinates(xiList, relativeThicknessList, lengthToDiameterRatio, Calculates organ coordinates and derivatives represented by a cylindrical tube with a unit inner diameter, length equivalent to lengthToDiameterRatio and wall thickness of wallThicknessToDiameterRatio. :param xiList: List containing xi for each point around the outer surface of the tube. - :param relativeThicknessList: Relative thickness of each element through wall. + :param relativeThicknessList: Relative thickness of each element through wall for organ coordinates. :param lengthToDiameterRatio: Ratio of total length along organ to inner diameter of organ :param wallThicknessToDiameterRatio: Ratio of wall thickness to inner diameter of organ. :param elementsCountAround: Number of elements around tube. From 873a18de517256bc568c9f29ac1192c70d163360 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 15:10:05 +1300 Subject: [PATCH 11/21] Increase mesenteric zone width in cattle colon --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 1c66bbf7..c5fc4853 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -295,9 +295,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Transverse length'] = 3500.0 options['Distal length'] = 1650.0 options['Proximal inner radius'] = 35.0 - options['Proximal tenia coli width'] = 3.0 + options['Proximal tenia coli width'] = 12.0 options['Proximal-transverse inner radius'] = 15.0 - options['Proximal-transverse tenia coli width'] = 3.0 + options['Proximal-transverse tenia coli width'] = 6.0 options['Transverse-distal inner radius'] = 9.0 options['Transverse-distal tenia coli width'] = 3.0 options['Distal inner radius'] = 10.5 From a2e5d91e54302ec652eccde1c714ff2f2511860c Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 18:05:48 +1300 Subject: [PATCH 12/21] Allow 1 to 4 elements through wall --- .../meshtypes/meshtype_3d_colon1.py | 74 +++++++++------- .../meshtypes/meshtype_3d_colonsegment1.py | 86 ++++++++++++------- 2 files changed, 99 insertions(+), 61 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index c5fc4853..865b8bee 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -6,7 +6,7 @@ import copy from opencmiss.zinc.element import Element -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ @@ -483,8 +483,8 @@ def generateBaseMesh(cls, region, options): # Colon coordinates lengthToDiameterRatio = 24 wallThicknessToDiameterRatio = 0.1 - teniaColiThicknessToDiameterRatio = 0.025 - relativeThicknessListColonCoordinates = [0.25, 0.25, 0.25, 0.25] + teniaColiThicknessToDiameterRatio = wallThicknessToDiameterRatio/elementsCountThroughWall + relativeThicknessListColonCoordinates = [1.0/elementsCountThroughWall for n3 in range(elementsCountThroughWall)] firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -576,20 +576,34 @@ def generateBaseMesh(cls, region, options): for n in range(elementsCount): annotationGroupsAlong.append(annotationGroupAlong[i]) - mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) - submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) - annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], [circularMuscleGroup], [longitudinalMuscleGroup]] - xExtrude = [] d1Extrude = [] d2Extrude = [] d3UnitExtrude = [] sxRefExtrudeList = [] - relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, circularRelThickness, - longitudinalRelThickness] + if elementsCountThroughWall > 1: + relativeThicknessList = [] + annotationGroupsThroughWall = [] + if mucosaRelThickness > 0.0: + relativeThicknessList.append(mucosaRelThickness) + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + annotationGroupsThroughWall.append([mucosaGroup]) + if submucosaRelThickness > 0.0: + relativeThicknessList.append(submucosaRelThickness) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + annotationGroupsThroughWall.append([submucosaGroup]) + if circularRelThickness > 0.0: + relativeThicknessList.append(circularRelThickness) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + annotationGroupsThroughWall.append([circularMuscleGroup]) + if longitudinalRelThickness > 0.0: + relativeThicknessList.append(longitudinalRelThickness) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + annotationGroupsThroughWall.append([longitudinalMuscleGroup]) + else: + relativeThicknessList = [1.0] + annotationGroupsThroughWall = [[]] # Create object colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( @@ -719,22 +733,22 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): ''' # Create 2d surface mesh groups fm = region.getFieldmodule() - longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - get_colon_term("Longitudinal muscle layer of colon")) - circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - get_colon_term("Circular muscle layer of colon")) - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) - is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) - is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - - serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) - serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_colon_term("myenteric nerve plexus")) - myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) + longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") + circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") + + if longitudinalMuscleGroup and circularMuscleGroup: + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + + is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) + is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) + is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) + + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + + myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, + get_colon_term("myenteric nerve plexus")) + myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index b1d63f00..0c085c93 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -11,7 +11,7 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear @@ -183,7 +183,7 @@ def checkOptions(options): options[key] = 2 if options['Number of elements around haustrum'] < 4: options['Number of elements around haustrum'] = 4 - if options['Number of elements through wall'] != 4: + if options['Number of elements through wall'] > 4: options['Number of elements through wall'] = 4 for key in [ 'Number of elements around tenia coli', @@ -305,12 +305,6 @@ def generateBaseMesh(cls, region, options): for i in range(elementsCountAlongSegment): annotationGroupsAlong.append([ ]) - mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) - submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) - annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], [circularMuscleGroup], [longitudinalMuscleGroup]] - # Create inner points nSegment = 0 closedProximalEnd = False @@ -330,8 +324,29 @@ def generateBaseMesh(cls, region, options): closedProximalEnd) contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() - relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, - circularRelThickness, longitudinalRelThickness] + + if elementsCountThroughWall > 1: + relativeThicknessList = [] + annotationGroupsThroughWall = [] + if mucosaRelThickness > 0.0: + relativeThicknessList.append(mucosaRelThickness) + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + annotationGroupsThroughWall.append([mucosaGroup]) + if submucosaRelThickness > 0.0: + relativeThicknessList.append(submucosaRelThickness) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + annotationGroupsThroughWall.append([submucosaGroup]) + if circularRelThickness > 0.0: + relativeThicknessList.append(circularRelThickness) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + annotationGroupsThroughWall.append([circularMuscleGroup]) + if longitudinalRelThickness > 0.0: + relativeThicknessList.append(longitudinalRelThickness) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + annotationGroupsThroughWall.append([longitudinalMuscleGroup]) + else: + relativeThicknessList = [1.0] + annotationGroupsThroughWall = [[]] # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, @@ -405,24 +420,29 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): ''' # Create 2d surface mesh groups fm = region.getFieldmodule() - longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - get_colon_term("Longitudinal muscle layer of colon")) - circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - get_colon_term("Circular muscle layer of colon")) - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) - is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) - is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) + longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") + circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") + + if longitudinalMuscleGroup and circularMuscleGroup: + # longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + # get_colon_term("Longitudinal muscle layer of colon")) + # circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, + # get_colon_term("Circular muscle layer of colon")) + mesh2d = fm.findMeshByDimension(2) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + + is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) + is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) + is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) - serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) - myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) + myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) + myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) class ColonSegmentTubeMeshInnerPoints: """ @@ -1887,6 +1907,7 @@ def createNodesAndElementsTeniaColi(region, # Create elements allAnnotationGroups = [] + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) if closedProximalEnd: elementtemplate3 = mesh.createElementtemplate() @@ -1965,7 +1986,8 @@ def createNodesAndElementsTeniaColi(region, element.setNodesByIdentifier(eft4 if tetrahedralElement else eft6, nodeIdentifiers) element.setScaleFactors(eft4 if tetrahedralElement else eft6, scalefactors) elementIdentifier = elementIdentifier + 1 - annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + \ + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2018,7 +2040,8 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( - elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[0] + elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[0] + \ + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2065,7 +2088,8 @@ def createNodesAndElementsTeniaColi(region, element.setNodesByIdentifier(eft5 if tetrahedralElement else eft6, nodeIdentifiers) element.setScaleFactors(eft5 if tetrahedralElement else eft6, scalefactors) elementIdentifier = elementIdentifier + 1 - annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + \ + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2144,7 +2168,7 @@ def createNodesAndElementsTeniaColi(region, element.setNodesByIdentifier(eftOrgan if eTC < int(elementsCountAroundTC*0.5) - 1 else eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] + \ - annotationGroupsThroughWall[-1] + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2213,7 +2237,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] + \ - annotationGroupsThroughWall[-1] + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2263,7 +2287,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[ elementsCountAround + int(elementsCountAroundTC * (tcCount - 0.5)) + eTC] + \ - annotationGroupsAlong[e2] + annotationGroupsThroughWall[-1] + annotationGroupsAlong[e2] + [longitudinalMuscleGroup] if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: From 71351f3b51395b6bc741cc71c1d13931e61f7f43 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 18:22:14 +1300 Subject: [PATCH 13/21] Annotate serosa using 2D outside surface of colon --- .../meshtypes/meshtype_3d_colon1.py | 20 +++++++------ .../meshtypes/meshtype_3d_colonsegment1.py | 29 +++++++++---------- 2 files changed, 24 insertions(+), 25 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 865b8bee..d5451ed0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -6,7 +6,7 @@ import copy from opencmiss.zinc.element import Element -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName, getAnnotationGroupForTerm from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ @@ -733,22 +733,24 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): ''' # Create 2d surface mesh groups fm = region.getFieldmodule() + mesh2d = fm.findMeshByDimension(2) + + colonGroup = getAnnotationGroupForTerm(annotationGroups, get_colon_term("colon")) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_colon = colonGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_colon, is_exterior_face_xi3_1) + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") if longitudinalMuscleGroup and circularMuscleGroup: - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) - serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 0c085c93..4b0be096 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -11,7 +11,7 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, mergeAnnotationGroups, findOrCreateAnnotationGroupForTerm, findAnnotationGroupByName, getAnnotationGroupForTerm from scaffoldmaker.annotation.colon_terms import get_colon_term from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear @@ -301,9 +301,10 @@ def generateBaseMesh(cls, region, options): radiusAlongSegment, dRadiusAlongSegment, tcWidthAlongSegment, startPhase) # Create annotation - annotationGroupsAlong = [] + colonGroup = AnnotationGroup(region, get_colon_term("colon")) + annotationGroupsAlong = [ ] for i in range(elementsCountAlongSegment): - annotationGroupsAlong.append([ ]) + annotationGroupsAlong.append([colonGroup]) # Create inner points nSegment = 0 @@ -420,27 +421,23 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): ''' # Create 2d surface mesh groups fm = region.getFieldmodule() + mesh2d = fm.findMeshByDimension(2) + + colonGroup = getAnnotationGroupForTerm(annotationGroups, get_colon_term("colon")) + is_exterior = fm.createFieldIsExterior() + is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_colon = colonGroup.getFieldElementGroup(mesh2d) + is_serosa = fm.createFieldAnd(is_colon, is_exterior_face_xi3_1) + serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) + serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") if longitudinalMuscleGroup and circularMuscleGroup: - # longitudinalMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - # get_colon_term("Longitudinal muscle layer of colon")) - # circularMuscleGroup = getAnnotationGroupForTerm(annotationGroups, - # get_colon_term("Circular muscle layer of colon")) - mesh2d = fm.findMeshByDimension(2) - is_exterior = fm.createFieldIsExterior() - is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_serosa = fm.createFieldAnd(is_longitudinalMuscle, is_exterior_face_xi3_1) is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - - serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) - serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) From 6146e3ac500cb1f47fff5dcd9a3df98c3f916288 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 18 Oct 2021 18:28:04 +1300 Subject: [PATCH 14/21] Remove annotation group for myenteric nerve plexus --- src/scaffoldmaker/annotation/colon_terms.py | 1 - src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 12 ------------ .../meshtypes/meshtype_3d_colonsegment1.py | 10 ---------- 3 files changed, 23 deletions(-) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index cc27bccf..8bab6ab4 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -13,7 +13,6 @@ ( "distal colon", "UBERON:0008971", "ILX:0727523"), ( "Longitudinal muscle layer of colon", "ILX:0775554"), ( "mesenteric zone", "None"), - ( "myenteric nerve plexus", "UBERON:0002439", "FMA:63252", "ILX:0725342"), ( "non-mesenteric zone", "None"), ( "proximal colon", "UBERON:0008972", "ILX:0733240"), ( "serosa of colon", "UBERON:0003335", "FMA:14990", "ILX:0736932"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index d5451ed0..b360d82f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -742,15 +742,3 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): is_serosa = fm.createFieldAnd(is_colon, is_exterior_face_xi3_1) serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - - longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") - circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") - - if longitudinalMuscleGroup and circularMuscleGroup: - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) - is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, - get_colon_term("myenteric nerve plexus")) - myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 4b0be096..d8581d72 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -431,16 +431,6 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) - longitudinalMuscleGroup = findAnnotationGroupByName(annotationGroups, "Longitudinal muscle layer of colon") - circularMuscleGroup = findAnnotationGroupByName(annotationGroups, "Circular muscle layer of colon") - - if longitudinalMuscleGroup and circularMuscleGroup: - is_longitudinalMuscle = longitudinalMuscleGroup.getFieldElementGroup(mesh2d) - is_circularMuscle = circularMuscleGroup.getFieldElementGroup(mesh2d) - is_myentericPlexus = fm.createFieldAnd(is_longitudinalMuscle, is_circularMuscle) - myentericPlexus = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("myenteric nerve plexus")) - myentericPlexus.getMeshGroup(mesh2d).addElementsConditional(is_myentericPlexus) - class ColonSegmentTubeMeshInnerPoints: """ Generates inner profile of a colon segment for use by tubemesh. From a99581c66c6751e0bdd35e92c12a569c4393f9fe Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 10:17:18 +1300 Subject: [PATCH 15/21] Make tenia coli thickness 0.25 of wall thickness in colon coordinates --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index b360d82f..25b9c398 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -483,7 +483,7 @@ def generateBaseMesh(cls, region, options): # Colon coordinates lengthToDiameterRatio = 24 wallThicknessToDiameterRatio = 0.1 - teniaColiThicknessToDiameterRatio = wallThicknessToDiameterRatio/elementsCountThroughWall + teniaColiThicknessToDiameterRatio = 0.25*wallThicknessToDiameterRatio relativeThicknessListColonCoordinates = [1.0/elementsCountThroughWall for n3 in range(elementsCountThroughWall)] firstNodeIdentifier = 1 From 27a77fc7c231abbf59db3dd89597da536cf84c5a Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 10:32:30 +1300 Subject: [PATCH 16/21] Recalculate thickness proportions for less than 4 elements through wall --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 3 +++ src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 25b9c398..d732b88e 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -601,6 +601,9 @@ def generateBaseMesh(cls, region, options): relativeThicknessList.append(longitudinalRelThickness) longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) annotationGroupsThroughWall.append([longitudinalMuscleGroup]) + totalProportions = sum(relativeThicknessList) + for i in range(len(relativeThicknessList)): + relativeThicknessList[i] = relativeThicknessList[i] / totalProportions else: relativeThicknessList = [1.0] annotationGroupsThroughWall = [[]] diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index d8581d72..946b9660 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -345,6 +345,9 @@ def generateBaseMesh(cls, region, options): relativeThicknessList.append(longitudinalRelThickness) longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) annotationGroupsThroughWall.append([longitudinalMuscleGroup]) + totalProportions = sum(relativeThicknessList) + for i in range(len(relativeThicknessList)): + relativeThicknessList[i] = relativeThicknessList[i] / totalProportions else: relativeThicknessList = [1.0] annotationGroupsThroughWall = [[]] From 8ffbd3c453609ae632cd2d75aeb960e5e6759f47 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 10:56:33 +1300 Subject: [PATCH 17/21] Add annotation for inner surface of colonic mucosa --- src/scaffoldmaker/annotation/colon_terms.py | 1 + src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 5 +++++ src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py | 4 ++++ 3 files changed, 10 insertions(+) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index 8bab6ab4..e887a5fb 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -11,6 +11,7 @@ ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"), ( "colonic mucosa", "UBERON:0000317", "FMA:14984", "ILX:0731046"), ( "distal colon", "UBERON:0008971", "ILX:0727523"), + ( "inner surface of colonic mucosa", "None"), ( "Longitudinal muscle layer of colon", "ILX:0775554"), ( "mesenteric zone", "None"), ( "non-mesenteric zone", "None"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index d732b88e..34024407 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -741,7 +741,12 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): colonGroup = getAnnotationGroupForTerm(annotationGroups, get_colon_term("colon")) is_exterior = fm.createFieldIsExterior() is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) is_colon = colonGroup.getFieldElementGroup(mesh2d) is_serosa = fm.createFieldAnd(is_colon, is_exterior_face_xi3_1) serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + is_mucosaInnerSurface = fm.createFieldAnd(is_colon, is_exterior_face_xi3_0) + mucosaInnerSurface = findOrCreateAnnotationGroupForTerm(annotationGroups, region, + get_colon_term("inner surface of colonic mucosa")) + mucosaInnerSurface.getMeshGroup(mesh2d).addElementsConditional(is_mucosaInnerSurface) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 946b9660..88393167 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -429,10 +429,14 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): colonGroup = getAnnotationGroupForTerm(annotationGroups, get_colon_term("colon")) is_exterior = fm.createFieldIsExterior() is_exterior_face_xi3_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + is_exterior_face_xi3_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) is_colon = colonGroup.getFieldElementGroup(mesh2d) is_serosa = fm.createFieldAnd(is_colon, is_exterior_face_xi3_1) serosa = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("serosa of colon")) serosa.getMeshGroup(mesh2d).addElementsConditional(is_serosa) + is_mucosaInnerSurface = fm.createFieldAnd(is_colon, is_exterior_face_xi3_0) + mucosaInnerSurface = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_colon_term("inner surface of colonic mucosa")) + mucosaInnerSurface.getMeshGroup(mesh2d).addElementsConditional(is_mucosaInnerSurface) class ColonSegmentTubeMeshInnerPoints: """ From b0505407353f08f51dd1b3cb66f787fee14e2fce Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 11:20:08 +1300 Subject: [PATCH 18/21] Remove longitudinal muscle annotation for tenia coli for 1 element through wall --- .../meshtypes/meshtype_3d_colonsegment1.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 88393167..6d038f92 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -1901,7 +1901,12 @@ def createNodesAndElementsTeniaColi(region, # Create elements allAnnotationGroups = [] - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + + for group in annotationGroupsThroughWall: + longitudinalMuscle = findAnnotationGroupByName(group, "Longitudinal muscle layer of colon") + + if longitudinalMuscle: + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) if closedProximalEnd: elementtemplate3 = mesh.createElementtemplate() @@ -1981,7 +1986,7 @@ def createNodesAndElementsTeniaColi(region, element.setScaleFactors(eft4 if tetrahedralElement else eft6, scalefactors) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + \ - [longitudinalMuscleGroup] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2035,7 +2040,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[0] + \ - [longitudinalMuscleGroup] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2083,7 +2088,7 @@ def createNodesAndElementsTeniaColi(region, element.setScaleFactors(eft5 if tetrahedralElement else eft6, scalefactors) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[0] + \ - [longitudinalMuscleGroup] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2162,7 +2167,7 @@ def createNodesAndElementsTeniaColi(region, element.setNodesByIdentifier(eftOrgan if eTC < int(elementsCountAroundTC*0.5) - 1 else eftOrgan1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] + \ - [longitudinalMuscleGroup] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2231,7 +2236,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[elementsCountAround + int( elementsCountAroundTC * 0.5) + N * elementsCountAroundTC + eTC] + annotationGroupsAlong[e2] + \ - [longitudinalMuscleGroup] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2281,7 +2286,7 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[ elementsCountAround + int(elementsCountAroundTC * (tcCount - 0.5)) + eTC] + \ - annotationGroupsAlong[e2] + [longitudinalMuscleGroup] + annotationGroupsAlong[e2] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: From b2cc273cc0ac0ee289e9ce17ee0adf40a734c9cf Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 11:37:00 +1300 Subject: [PATCH 19/21] Remove upper case --- src/scaffoldmaker/annotation/colon_terms.py | 4 ++-- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 4 ++-- src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index e887a5fb..8680a69b 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -7,12 +7,12 @@ ( "ascending colon", "UBERON:0001156", "FMA:14545", "ILX:0734393"), ( "descending colon", "UBERON:0001158", "FMA:14547", "ILX:0724444"), ( "caecum", "UBERON:0001153", "FMA:14541", "ILX:0732270"), - ( "Circular muscle layer of colon", "ILX:0772428"), + ( "circular muscle layer of colon", "ILX:0772428"), ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"), ( "colonic mucosa", "UBERON:0000317", "FMA:14984", "ILX:0731046"), ( "distal colon", "UBERON:0008971", "ILX:0727523"), ( "inner surface of colonic mucosa", "None"), - ( "Longitudinal muscle layer of colon", "ILX:0775554"), + ( "longitudinal muscle layer of colon", "ILX:0775554"), ( "mesenteric zone", "None"), ( "non-mesenteric zone", "None"), ( "proximal colon", "UBERON:0008972", "ILX:0733240"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 34024407..8373c2cd 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -595,11 +595,11 @@ def generateBaseMesh(cls, region, options): annotationGroupsThroughWall.append([submucosaGroup]) if circularRelThickness > 0.0: relativeThicknessList.append(circularRelThickness) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) annotationGroupsThroughWall.append([circularMuscleGroup]) if longitudinalRelThickness > 0.0: relativeThicknessList.append(longitudinalRelThickness) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) annotationGroupsThroughWall.append([longitudinalMuscleGroup]) totalProportions = sum(relativeThicknessList) for i in range(len(relativeThicknessList)): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 6d038f92..ead73548 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -339,11 +339,11 @@ def generateBaseMesh(cls, region, options): annotationGroupsThroughWall.append([submucosaGroup]) if circularRelThickness > 0.0: relativeThicknessList.append(circularRelThickness) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("Circular muscle layer of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) annotationGroupsThroughWall.append([circularMuscleGroup]) if longitudinalRelThickness > 0.0: relativeThicknessList.append(longitudinalRelThickness) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) annotationGroupsThroughWall.append([longitudinalMuscleGroup]) totalProportions = sum(relativeThicknessList) for i in range(len(relativeThicknessList)): @@ -1903,10 +1903,10 @@ def createNodesAndElementsTeniaColi(region, allAnnotationGroups = [] for group in annotationGroupsThroughWall: - longitudinalMuscle = findAnnotationGroupByName(group, "Longitudinal muscle layer of colon") + longitudinalMuscle = findAnnotationGroupByName(group, "longitudinal muscle layer of colon") if longitudinalMuscle: - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("Longitudinal muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) if closedProximalEnd: elementtemplate3 = mesh.createElementtemplate() From f5eec0ab8ab91b9eb5b38cac6bff9fc6605947be Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 11:40:21 +1300 Subject: [PATCH 20/21] Update unit tests --- tests/test_colon.py | 2 +- tests/test_colonsegment.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_colon.py b/tests/test_colon.py index 5020d32c..d73a87a9 100644 --- a/tests/test_colon.py +++ b/tests/test_colon.py @@ -124,7 +124,7 @@ def test_colon1(self): colonCoordinates = fieldmodule.findFieldByName("colon coordinates").castFiniteElement() minimums, maximums = evaluateFieldNodesetRange(colonCoordinates, nodes) assertAlmostEqualList(self, minimums, [-0.6, 0.0, -0.6], 1.0E-4) - assertAlmostEqualList(self, maximums, [ 0.6, 24.0, 0.605], 1.0E-4) + assertAlmostEqualList(self, maximums, [ 0.6, 24.0, 0.625], 1.0E-4) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) diff --git a/tests/test_colonsegment.py b/tests/test_colonsegment.py index 725c0772..2f010ba6 100644 --- a/tests/test_colonsegment.py +++ b/tests/test_colonsegment.py @@ -38,7 +38,7 @@ def test_humancolonsegment1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options) - self.assertEqual(7, len(annotationGroups)) + self.assertEqual(8, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) @@ -89,7 +89,7 @@ def test_mousecolonsegment1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options) - self.assertEqual(6, len(annotationGroups)) + self.assertEqual(7, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) From c20a504817fe05eb171b4230d6573c4d3d243f3a Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Oct 2021 14:28:04 +1300 Subject: [PATCH 21/21] Limit number of elements through wall to 1 and 4 --- .../meshtypes/meshtype_3d_colon1.py | 33 ++++++----------- .../meshtypes/meshtype_3d_colonsegment1.py | 35 ++++++------------- 2 files changed, 21 insertions(+), 47 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 8373c2cd..e952207b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -582,31 +582,18 @@ def generateBaseMesh(cls, region, options): d3UnitExtrude = [] sxRefExtrudeList = [] - if elementsCountThroughWall > 1: - relativeThicknessList = [] - annotationGroupsThroughWall = [] - if mucosaRelThickness > 0.0: - relativeThicknessList.append(mucosaRelThickness) - mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) - annotationGroupsThroughWall.append([mucosaGroup]) - if submucosaRelThickness > 0.0: - relativeThicknessList.append(submucosaRelThickness) - submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) - annotationGroupsThroughWall.append([submucosaGroup]) - if circularRelThickness > 0.0: - relativeThicknessList.append(circularRelThickness) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) - annotationGroupsThroughWall.append([circularMuscleGroup]) - if longitudinalRelThickness > 0.0: - relativeThicknessList.append(longitudinalRelThickness) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) - annotationGroupsThroughWall.append([longitudinalMuscleGroup]) - totalProportions = sum(relativeThicknessList) - for i in range(len(relativeThicknessList)): - relativeThicknessList[i] = relativeThicknessList[i] / totalProportions - else: + if elementsCountThroughWall == 1: relativeThicknessList = [1.0] annotationGroupsThroughWall = [[]] + else: + relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, + circularRelThickness, longitudinalRelThickness] + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) + annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], + [circularMuscleGroup], [longitudinalMuscleGroup]] # Create object colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index ead73548..9469f7a2 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -183,7 +183,7 @@ def checkOptions(options): options[key] = 2 if options['Number of elements around haustrum'] < 4: options['Number of elements around haustrum'] = 4 - if options['Number of elements through wall'] > 4: + if options['Number of elements through wall'] != (1 or 4): options['Number of elements through wall'] = 4 for key in [ 'Number of elements around tenia coli', @@ -326,31 +326,18 @@ def generateBaseMesh(cls, region, options): contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() - if elementsCountThroughWall > 1: - relativeThicknessList = [] - annotationGroupsThroughWall = [] - if mucosaRelThickness > 0.0: - relativeThicknessList.append(mucosaRelThickness) - mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) - annotationGroupsThroughWall.append([mucosaGroup]) - if submucosaRelThickness > 0.0: - relativeThicknessList.append(submucosaRelThickness) - submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) - annotationGroupsThroughWall.append([submucosaGroup]) - if circularRelThickness > 0.0: - relativeThicknessList.append(circularRelThickness) - circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) - annotationGroupsThroughWall.append([circularMuscleGroup]) - if longitudinalRelThickness > 0.0: - relativeThicknessList.append(longitudinalRelThickness) - longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) - annotationGroupsThroughWall.append([longitudinalMuscleGroup]) - totalProportions = sum(relativeThicknessList) - for i in range(len(relativeThicknessList)): - relativeThicknessList[i] = relativeThicknessList[i] / totalProportions - else: + if elementsCountThroughWall == 1: relativeThicknessList = [1.0] annotationGroupsThroughWall = [[]] + else: + relativeThicknessList = [mucosaRelThickness, submucosaRelThickness, + circularRelThickness, longitudinalRelThickness] + mucosaGroup = AnnotationGroup(region, get_colon_term("colonic mucosa")) + submucosaGroup = AnnotationGroup(region, get_colon_term("submucosa of colon")) + circularMuscleGroup = AnnotationGroup(region, get_colon_term("circular muscle layer of colon")) + longitudinalMuscleGroup = AnnotationGroup(region, get_colon_term("longitudinal muscle layer of colon")) + annotationGroupsThroughWall = [[mucosaGroup], [submucosaGroup], + [circularMuscleGroup], [longitudinalMuscleGroup]] # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList,