diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index 0c4fc010..8680a69b 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -7,9 +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"), ( "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"), ( "proximal colon", "UBERON:0008972", "ILX:0733240"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 7ecbf4d0..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, @@ -902,8 +904,8 @@ def generateBaseMesh(cls, region, options): d2Final += d2List[(elementsCountThroughWall + 1) * elementsCountAround:] d3Final += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] - xFlat = d1Flat = d2Flat = d3Flat = [] - xTexture = d1Texture = d2Texture = d3Texture = [] + 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) @@ -1004,7 +1006,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, 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 78acad65..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 @@ -463,8 +464,8 @@ def generateBaseMesh(cls, region, options): d2Cecum += d2List[(elementsCountThroughWall + 1) * elementsCountAround:] d3Cecum += d3List[(elementsCountThroughWall + 1) * elementsCountAround:] - xFlat = d1Flat = d2Flat = d3Flat = [] - xTexture = d1Texture = d2Texture = d3Texture = [] + xFlat = d1Flat = d2Flat = [] + xOrgan = d1Organ = d2Organ = [] # Create nodes and elements if tcThickness > 0: @@ -475,7 +476,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, xOrgan, d1Organ, d2Organ, None, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroupsAround, annotationGroupsAlong, annotationGroupsThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, @@ -483,7 +484,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, 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 9f445a37..e952207b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -5,11 +5,12 @@ """ import copy -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from opencmiss.zinc.element import Element +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,\ - getTeniaColi, createFlatAndTextureCoordinatesTeniaColi, 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 @@ -294,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 @@ -313,11 +314,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 @@ -471,10 +472,20 @@ 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) + # Colon coordinates + lengthToDiameterRatio = 24 + wallThicknessToDiameterRatio = 0.1 + teniaColiThicknessToDiameterRatio = 0.25*wallThicknessToDiameterRatio + relativeThicknessListColonCoordinates = [1.0/elementsCountThroughWall for n3 in range(elementsCountThroughWall)] + firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -565,16 +576,25 @@ def generateBaseMesh(cls, region, options): for n in range(elementsCount): annotationGroupsAlong.append(annotationGroupAlong[i]) - annotationGroupsThroughWall = [] - for i in range(elementsCountThroughWall): - annotationGroupsThroughWall.append([ ]) - xExtrude = [] d1Extrude = [] d2Extrude = [] d3UnitExtrude = [] sxRefExtrudeList = [] + 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( region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, @@ -613,7 +633,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() @@ -628,30 +648,48 @@ def generateBaseMesh(cls, region, options): tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroupsAround, closedProximalEnd) - # Create flat and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( - xiList, relaxedLengthList, length, wallThickness, tcCount, tcThickness, + # Create flat coordinates + xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( + xiList, relaxedLengthList, length, wallThickness, relativeThicknessList, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd) + # Create colon coordinates + xColon, d1Colon, d2Colon = createColonCoordinatesTeniaColi(xiList, relativeThicknessListColonCoordinates, + 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, xTexture, d1Texture, d2Texture, - 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 and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( - xiList, relaxedLengthList, length, wallThickness, elementsCountAround, + # Create flat coordinates + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( + xiList, relaxedLengthList, length, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) + # Create colon coordinates + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessListColonCoordinates, + lengthToDiameterRatio, + wallThicknessToDiameterRatio, + elementsCountAround, + elementsCountAlong, elementsCountThroughWall, + transitElementList) + # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, - 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) @@ -672,3 +710,30 @@ 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() + 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_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 04c05ebe..9469f7a2 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, 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 @@ -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', @@ -153,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']: @@ -166,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'] != (1 or 4): + options['Number of elements through wall'] = 4 for key in [ 'Number of elements around tenia coli', 'Number of elements around haustrum']: @@ -179,7 +198,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 +255,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 @@ -274,18 +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([ ]) - - # 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([ ]) + annotationGroupsAlong.append([colonGroup]) # Create inner points nSegment = 0 @@ -307,11 +326,26 @@ def generateBaseMesh(cls, region, options): contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() + 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, - d2WarpedList, d3WarpedUnitList, contractedWallThicknessList, + d2WarpedList, d3WarpedUnitList, contractedWallThicknessList, relativeThicknessList, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) + xColonSegment = d1ColonSegment = d2ColonSegment = [] + relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() if tcThickness > 0: @@ -321,29 +355,29 @@ 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( - xiList, relaxedLengthList, segmentLength, wallThickness, tcCount, tcThickness, + # Create flat coordinates + xFlat, d1Flat, d2Flat = createFlatCoordinatesTeniaColi( + xiList, relaxedLengthList, segmentLength, wallThickness, relativeThicknessList, 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, - 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 and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( - xiList, relaxedLengthList, segmentLength, wallThickness, elementsCountAround, + # Create flat coordinates + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( + xiList, relaxedLengthList, segmentLength, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, - 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) @@ -365,6 +399,32 @@ 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() + 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_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: """ Generates inner profile of a colon segment for use by tubemesh. @@ -1177,7 +1237,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 +1483,12 @@ def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, return x, d1, d2, d3 -def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, - totalLengthAlong, wallThickness, tcCount, tcThickness, +def createFlatCoordinatesTeniaColi(xiList, relaxedLengthList, + totalLengthAlong, wallThickness, relativeThicknessList, 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. @@ -1437,6 +1496,7 @@ def createFlatAndTextureCoordinatesTeniaColi(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. @@ -1446,18 +1506,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, relativeThicknessList, + elementsCountAround, elementsCountAlong, + elementsCountThroughWall, transitElementList) # Find flat coordinates for tenia coli xFlatListTC = [] @@ -1514,79 +1574,103 @@ 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 + +def createColonCoordinatesTeniaColi(xiList, relativeThicknessList, 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 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. + :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 - return xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture + # Find organ coordinates for colon + xColon, d1Colon, d2Colon = tubemesh.createOrganCoordinates(xiList, relativeThicknessList, 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, - xTexture, d1Texture, d2Texture, + xOrgan, d1Organ, d2Organ, organCoordinateFieldName, 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 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. @@ -1653,11 +1737,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,63 +1764,53 @@ 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) + 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: - 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) - + 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) @@ -1754,7 +1828,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 +1837,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,32 +1858,43 @@ 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 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 = [] + 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() elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) @@ -1896,7 +1972,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 longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -1949,7 +2026,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 longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -1996,7 +2074,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 longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2029,10 +2108,12 @@ 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) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ annotationGroupsThroughWall[e3] @@ -2065,12 +2146,15 @@ 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) + 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] + annotationGroups = annotationGroupsAround[elementsCountAround + eTC] + annotationGroupsAlong[e2] + \ + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2108,31 +2192,38 @@ 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) + 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] 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) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, 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) + if xOrgan: + element.merge(organElementtemplate1) + 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] + \ + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: @@ -2169,19 +2260,20 @@ 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) + 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] + \ - annotationGroupsAlong[e2] + annotationGroupsAlong[e2] + ([longitudinalMuscleGroup] if longitudinalMuscle else []) if annotationGroups: allAnnotationGroups = mergeAnnotationGroups(allAnnotationGroups, annotationGroups) for annotationGroup in annotationGroups: diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 13c4ca59..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,19 +917,21 @@ 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 and texture coordinates - xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( - xiList, flatWidthList, length, wallThickness, elementsCountAround, + # Create flat coordinates + xFlat, d1Flat, d2Flat = tubemesh.createFlatCoordinates( + xiList, flatWidthList, length, wallThickness, relativeThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) + xOrgan = d1Organ = d2Organ = [] + # 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, 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 cc552fd3..8a3caf1f 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,21 +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 createFlatAndTextureCoordinates(xiList, lengthAroundList, - totalLengthAlong, wallThickness, elementsCountAround, - elementsCountAlong, elementsCountThroughWall, transitElementList): +def createFlatCoordinates(xiList, lengthAroundList, totalLengthAlong, wallThickness, relativeThicknessList, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList): """ - Calculates flat coordinates and texture 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. + 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 totalLengthAlong: Total length along tube. :param wallThickness: Thickness of wall. :param elementsCountAround: Number of elements around tube. @@ -352,8 +357,15 @@ 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. """ + 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 = [] @@ -374,7 +386,7 @@ def createFlatAndTextureCoordinates(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, @@ -397,55 +409,80 @@ 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] + return xFlatList, d1FlatList, d2FlatList - 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] +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 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. + :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. + """ + 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 = [] + 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): - 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 + 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 = [] + 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, - xTexture, d1Texture, d2Texture, + xOrgan, d1Organ, d2Organ, organCoordinateFieldName, 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 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. @@ -499,8 +536,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 +558,29 @@ 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) + flatElementtemplate2.defineField(flatCoordinates, -1, eftFlat2) - 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 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: - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) + organNodetemplate.setValueNumberOfVersions(organCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - 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) + organElementtemplate = mesh.createElementtemplate() + organElementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + organElementtemplate.defineField(organCoordinates, -1, eftOrgan) # Create nodes # Coordinates field @@ -571,8 +599,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,31 +608,36 @@ 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 + # 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) @@ -667,10 +700,12 @@ 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) + if xOrgan: + element.merge(organElementtemplate) + element.setNodesByIdentifier(eftOrgan, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 annotationGroups = annotationGroupsAround[e1] + annotationGroupsAlong[e2] + \ diff --git a/tests/test_colon.py b/tests/test_colon.py index cafa05f7..d73a87a9 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.625], 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..2f010ba6 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(8, 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(7, 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()