From 2433b2ac50f5b0ef24c580f23e2e6c75f41d9c29 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 16 Jun 2020 12:44:32 +1200 Subject: [PATCH 01/12] Edit parameters to better represent pig colon --- .../meshtypes/meshtype_3d_colon1.py | 84 +++++++++---------- .../meshtypes/meshtype_3d_colonsegment1.py | 2 +- 2 files changed, 43 insertions(+), 43 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 88166050..e7a2e3bc 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -100,43 +100,43 @@ class MeshType_3d_colon1(Scaffold_base): }, 'meshEdits' : exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ 163.7, -25.2, 12.2 ], [ -21.7, 50.1, -18.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 117.2, 32.8, -2.6 ], [ -64.3, 34.4, -3.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 38.8, 68.0, 2.6 ], [ -83.9, 24.9, 14.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -42.8, 65.5, 15.5 ], [ -69.9, -46.7, 10.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -93.1, -11.4, 30.6 ], [ 5.7, -104.4, 9.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -51.1, -83.7, 32.0 ], [ 69.7, -44.1, 1.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 43.0, -88.5, 28.5 ], [ 71.0, 41.0, -1.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 94.9, -5.8, 33.7 ], [ -14.3, 60.1, 0.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 56.6, 72.9, 45.0 ], [ -108.7, 54.1, 22.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -37.0, 64.6, 59.3 ], [ -72.5, -33.9, 6.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -94.1, -15.6, 78.2 ], [ 8.5, -82.1, 13.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -46.1, -90.1, 82.9 ], [ 76.6, -39.9, 11.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 45.8, -87.7, 92.3 ], [ 99.1, 55.4, 12.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 110.7, 2.2, 102.2 ], [ 0.4, 61.7, 12.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 73.9, 58.3, 110.7 ], [ -86.3, 19.7, 9.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -21.8, 25.4, 118.0 ], [ -78.1, -83.9, 23.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -41.7, -49.8, 134.9 ], [ 10.8, -41.9, 18.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 4.0, -46.5, 133.5 ], [ 8.3, 43.8, -7.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 12.0, 2.6, 130.0 ], [ 27.8, 34.6, 1.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 67.4, 25.7, 125.9 ], [ 66.6, -38.7, -10.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 75.1, -21.4, 122.5 ], [ -17.5, -31.0, 4.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 42.7, -66.4, 124.2 ], [ -47.7, -41.1, -2.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -2.1, -88.6, 122.8 ], [ -29.6, -5.3, -1.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -52.1, -78.6, 117.3 ], [ -42.2, 33.8, -6.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -72.2, -6.5, 104.4 ], [ 16.8, 80.4, -10.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -28.6, 45.7, 90.9 ], [ 47.4, 37.5, -9.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 67.1, 63.6, 78.6 ], [ 77.1, -44.5, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 89.1, -11.8, 67.7 ], [ 0.0, -73.3, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 40.5, -79.2, 62.3 ], [ -63.5, -36.7, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -36.4, -80.2, 59.6 ], [ -60.1, 26.4, -8.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -69.2, -19.3, 47.0 ], [ 0.0, 73.3, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -22.4, 39.1, 34.8 ], [ 63.5, 36.7, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 40.6, 46.0, 32.9 ], [ 25.4, -23.1, -7.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 66.2, 11.0, 24.6 ], [ 13.7, -51.4, -1.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 32.0, -49.3, 24.4 ], [ -51.1, -50.1, -1.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -51.1, -45.3, 15.7 ], [ -38.0, 52.2, -16.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -122.9, 124.2, -36.0 ], [ -21.3, 64.5, -18.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ] ] ) + [ [ -7.2, 83.3, -20.7 ], [ -65.2, -8.1, 7.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -68.5, 52.8, -9.6 ], [ -40.1, -36.1, 10.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -97.4, -26.3, 5.7 ], [ 18.0, -93.2, 13.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -56.8, -90.5, 14.1 ], [ 65.5, -41.4, 7.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 48.9, -100.8, 24.0 ], [ 112.2, 40.1, 19.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 114.8, -12.6, 38.7 ], [ 8.2, 96.1, 14.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 60.3, 83.5, 43.7 ], [ -108.7, 54.1, 22.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -41.2, 90.7, 56.3 ], [ -89.0, -32.4, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -107.9, -9.7, 76.6 ], [ 11.1, -94.4, 11.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -57.3, -91.9, 81.3 ], [ 71.2, -31.2, 5.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 51.2, -89.4, 97.2 ], [ 99.1, 55.4, 12.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 91.6, 9.3, 103.6 ], [ 4.7, 51.2, 3.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 61.6, 111.8, 109.6 ], [ -85.2, 46.1, 2.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -54.6, 91.9, 129.4 ], [ -92.7, -55.0, 14.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -109.0, 5.6, 156.9 ], [ 23.6, -108.2, 27.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -59.1, -62.5, 170.8 ], [ 74.0, -20.1, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 23.5, -53.2, 179.7 ], [ 84.6, 47.0, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 62.3, 30.1, 187.5 ], [ -12.8, 58.0, 0.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 22.4, 45.2, 181.1 ], [ -23.6, -34.5, -7.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -1.9, 4.9, 180.5 ], [ -41.3, -30.9, 7.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -45.1, -12.6, 194.4 ], [ -40.5, -4.6, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -71.7, -2.2, 197.2 ], [ -25.2, 35.8, -6.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -65.8, 42.1, 182.3 ], [ 26.6, 37.6, -15.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -14.1, 81.2, 163.5 ], [ 41.0, 10.3, -9.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 61.7, 86.1, 156.4 ], [ 77.9, -40.7, 8.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 92.9, 20.5, 150.3 ], [ 0.0, -73.3, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 48.9, -65.0, 142.8 ], [ -82.8, -80.0, -1.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -54.3, -90.8, 134.0 ], [ -60.1, 26.4, -8.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -89.9, 11.2, 115.0 ], [ 34.9, 125.1, -27.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -17.4, 74.2, 91.1 ], [ 78.8, 19.1, -15.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 43.4, 50.2, 73.3 ], [ 30.2, -36.0, -9.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 62.4, -5.1, 63.5 ], [ 10.9, -54.2, -2.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 32.7, -51.7, 56.1 ], [ -38.6, -29.8, -8.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -38.1, -28.6, 46.8 ], [ -66.0, 83.3, -12.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -5.4, 32.0, 26.0 ], [ 48.1, 17.6, -21.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 146.6, 101.2, -41.2 ], [ 63.3, 35.3, -31.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 312.3, 199.5, -123.4 ], [ 39.7, 24.3, -20.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ] ] ) } ), 'Pig 2': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { @@ -228,14 +228,14 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Distal tenia coli width'] = 1.0 elif 'Pig 1' in parameterSetName: options['Number of segments'] = 120 - options['Proximal length'] = 2610.0 + options['Proximal length'] = 3000.0 options['Transverse length'] = 200.0 options['Distal length'] = 200.0 - options['Proximal inner radius'] = 20.0 + options['Proximal inner radius'] = 38.0 options['Proximal tenia coli width'] = 5.0 - options['Proximal-transverse inner radius'] = 10.5 + options['Proximal-transverse inner radius'] = 14.0 options['Proximal-transverse tenia coli width'] = 5.0 - options['Transverse-distal inner radius'] = 8.0 + options['Transverse-distal inner radius'] = 10.5 options['Transverse-distal tenia coli width'] = 5.0 options['Distal inner radius'] = 8.0 options['Distal tenia coli width'] = 5.0 @@ -391,7 +391,7 @@ def generateBaseMesh(cls, region, options): centralPath.generate(tmpRegion) cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion) # for i in range(len(cx)): - # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') + # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') del tmpRegion # find arclength of colon diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 8a6c96e9..2f375488 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -95,7 +95,7 @@ def getDefaultOptions(parameterSetName='Default'): options['Start inner radius'] = 20.0 options['End inner radius'] = 20.0 options['Corner inner radius factor'] = 0.0 - options['Haustrum inner radius factor'] = 0.3 + options['Haustrum inner radius factor'] = 0.2 options['Segment length end derivative factor'] = 0.8 options['Segment length mid derivative factor'] = 2.0 options['Segment length'] = 25.0 From 79679920d366b843e9eb36f4d09eb00a898d3197 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 17 Jun 2020 14:46:22 +1200 Subject: [PATCH 02/12] Create annotation groups for proximal, transverse and distal colon --- src/scaffoldmaker/annotation/colon_terms.py | 20 +++-- .../meshtypes/meshtype_3d_cecum1.py | 7 +- .../meshtypes/meshtype_3d_colon1.py | 54 ++++++++++-- .../meshtypes/meshtype_3d_colonsegment1.py | 85 ++++++++++--------- .../meshtypes/meshtype_3d_smallintestine1.py | 6 +- src/scaffoldmaker/utils/tubemesh.py | 8 +- tests/test_colon.py | 4 +- 7 files changed, 118 insertions(+), 66 deletions(-) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index f59e6e3e..9f172a35 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -4,13 +4,19 @@ # convention: preferred name, preferred id, followed by any other ids and alternative names colon_terms = [ - ( "colon", "UBERON:0001155", "FMA:14543" ), - ( "mesenteric zone", None ), - ( "non-mesenteric zone", None ), - ( "taenia coli", "UBERON:0012419", "FMA:15041" ), - ( "tenia libera", None ), - ( "tenia mesocolica", None ), - ( "tenia omentalis", None ) + ( "ascending colon", "UBERON:0001156", "FMA:14545", "ILX:0734393"), + ( "descending colon", "UBERON:0001158", "FMA:14547", "ILX:0724444"), + ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"), + ( "distal colon", "UBERON:0008971", "ILX:0727523"), + ( "mesenteric zone", None), + ( "non-mesenteric zone", None), + ( "proximal colon", "UBERON:0008972", "ILX:0733240"), + ( "spiral colon", "UBERON:0010239", "ILX:0735018"), + ( "taenia coli", "UBERON:0012419", "FMA:15041", "ILX:0731555"), + ( "taenia libera", "ILX:0739285"), + ( "taenia mesocolica", "ILX:0739284"), + ( "taenia omentalis", "ILX:0739286"), + ( "transverse colon", "UBERON:0001157", "FMA:14546", "ILX:0728767") ] def get_colon_term(name : str): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 417a240c..e3cb2087 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -14,7 +14,6 @@ from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import matrix -from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils.tracksurface import TrackSurface, TrackSurfacePosition from scaffoldmaker.utils import tubemesh from scaffoldmaker.utils import vector @@ -359,9 +358,11 @@ def generateBaseMesh(cls, region, options): segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, innerRadiusAlongCecum, dInnerRadiusAlongCecum, tcWidthAlongCecum, startPhase) + annotationArrayAlong = [''] * elementsCountAlongSegment + for nSegment in range(segmentCount): # Make regular segments - xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray \ + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround \ = colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) # Replace first half of first segment with apex and sample along apex and second half of segment @@ -446,7 +447,7 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=True) # Add ostium on track surface between two tenia on the last segment diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index e7a2e3bc..21deae24 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -5,11 +5,13 @@ """ import copy +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +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 +from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1, ColonSegmentTubeMeshInnerPoints,\ + getTeniaColi, createFlatAndTextureCoordinatesTeniaColi, createNodesAndElementsTeniaColi from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage -from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh from scaffoldmaker.utils.zinc_utils import exnodeStringFromNodeValues @@ -405,6 +407,7 @@ def generateBaseMesh(cls, region, options): length += arcLength segmentLength = length / segmentCount # print('Length = ', length) + elementAlongLength = length / elementsCountAlong # Sample central path sx, sd1, se, sxi, ssf = interp.sampleCubicHermiteCurves(cx, cd1, elementsCountAlong) @@ -423,6 +426,37 @@ def generateBaseMesh(cls, region, options): tcWidthList, elementsCountAlong) + # Create annotation groups for colon sections + elementsAlongInProximal = round(proximalLength/elementAlongLength) + elementsAlongInTransverse = round(transverseLength/elementAlongLength) + elementsAlongInDistal = elementsCountAlong - elementsAlongInProximal - elementsAlongInTransverse + + if tcCount == 1: + proximalGroup = AnnotationGroup(region, get_colon_term("proximal colon")) + transverseGroup = AnnotationGroup(region, get_colon_term("transverse colon")) + distalGroup = AnnotationGroup(region, get_colon_term("distal colon")) + annotationGroups = [proximalGroup, transverseGroup, distalGroup] + annotationArrayAlong = (['proximal colon'] * elementsAlongInProximal + + ['transverse colon'] * elementsAlongInTransverse + + ['distal colon'] * elementsAlongInDistal) + + elif tcCount == 2: + spiralGroup = AnnotationGroup(region, get_colon_term("spiral colon")) + transverseGroup = AnnotationGroup(region, get_colon_term("transverse colon")) + distalGroup = AnnotationGroup(region, get_colon_term("distal colon")) + annotationGroups = [spiralGroup, transverseGroup, distalGroup] + annotationArrayAlong = (['spiral colon'] * elementsAlongInProximal + + ['transverse colon'] * elementsAlongInTransverse + + ['distal colon'] * elementsAlongInDistal) + elif tcCount == 3: + ascendingGroup = AnnotationGroup(region, get_colon_term("ascending colon")) + transverseGroup = AnnotationGroup(region, get_colon_term("transverse colon")) + descendingGroup = AnnotationGroup(region, get_colon_term("descending colon")) + annotationGroups = [ascendingGroup, transverseGroup, descendingGroup] + annotationArrayAlong = (['ascending colon'] * elementsAlongInProximal + + ['transverse colon'] * elementsAlongInTransverse + + ['descending colon'] * elementsAlongInDistal) + xExtrude = [] d1Extrude = [] d2Extrude = [] @@ -439,7 +473,7 @@ def generateBaseMesh(cls, region, options): for nSegment in range(segmentCount): # Create inner points - xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray\ + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroupsAround, annotationArrayAround\ = colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) # Project reference point for warping onto central path @@ -465,6 +499,8 @@ def generateBaseMesh(cls, region, options): contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() + annotationGroups += annotationGroupsAround + # Create coordinates and derivatives xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xExtrude, d1Extrude, d2Extrude, d3UnitExtrude, contractedWallThicknessList, @@ -474,10 +510,10 @@ def generateBaseMesh(cls, region, options): if tcThickness > 0: tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() - xList, d1List, d2List, d3List, annotationGroups, annotationArray = getTeniaColi( + xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi( region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroups, annotationArray) + tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroups, annotationArrayAround) # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( @@ -489,8 +525,8 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tcCount, annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives) + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives) else: # Create flat and texture coordinates @@ -502,8 +538,8 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 2f375488..bff51eab 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -18,7 +18,6 @@ from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.geometry import createCirclePoints from scaffoldmaker.utils import matrix -from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh from scaffoldmaker.utils import vector @@ -259,10 +258,12 @@ def generateBaseMesh(cls, region, options): segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, radiusAlongSegment, dRadiusAlongSegment, tcWidthAlongSegment, startPhase) + annotationArrayAlong = [''] * elementsCountAlongSegment + # Create inner points nSegment = 0 - xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray = \ + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround = \ colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) # Project reference point for warping onto central path @@ -287,10 +288,10 @@ def generateBaseMesh(cls, region, options): if tcThickness > 0: tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() - xList, d1List, d2List, d3List, annotationGroups, annotationArray = getTeniaColi( + xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi( region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, - tubeTCWidthList, tcThickness, sxRefList, annotationGroups, annotationArray) + tubeTCWidthList, tcThickness, sxRefList, annotationGroups, annotationArrayAround) # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( @@ -302,8 +303,8 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, - tcCount, annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives) + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives) else: # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( @@ -314,8 +315,8 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, - annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) return annotationGroups @@ -1200,9 +1201,8 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tubeTCWidthList, tcThickness, sx, - annotationGroups, annotationArray): + annotationGroups, annotationArrayAround): """ - EDIT Create equally spaced points for tenia coli over the outer surface of the colon. Points are sampled from a cubic hermite curve running through the left edge of tenia coli @@ -1224,7 +1224,7 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, :param tubeTCWidthList: List of tenia coli width along tube length. :param tcThickness: Thickness of tenia coli at its thickest part. :param sx: Coordinates of central path. - :param annotationGroups, annotationArray: annotations for tube + :param annotationGroups, annotationArrayAround: annotations for tube :return: coordinates, derivatives, annotationGroups for colon with tenia coli. """ @@ -1298,19 +1298,23 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, elementsCountAlong, elementsCountThroughWall) # Update annotation groups + if tcCount == 2: + tcGroup = AnnotationGroup(region, get_colon_term("taenia coli")) + annotationGroups += [tcGroup] + annotationArrayAround += (['taenia coli']*elementsCountAroundTC*tcCount) if tcCount == 3: - tlGroup = AnnotationGroup(region, get_colon_term("tenia libera")) - tmGroup = AnnotationGroup(region, get_colon_term("tenia mesocolica")) - toGroup = AnnotationGroup(region, get_colon_term("tenia omentalis")) + tlGroup = AnnotationGroup(region, get_colon_term("taenia libera")) + tmGroup = AnnotationGroup(region, get_colon_term("taenia mesocolica")) + toGroup = AnnotationGroup(region, get_colon_term("taenia omentalis")) annotationGroupsTC = [tlGroup, tmGroup, toGroup] - annotationArrayTC = (['tenia omentalis']*(elementsCountAroundTC//2) + - ['tenia libera']*elementsCountAroundTC + - ['tenia mesocolica']*elementsCountAroundTC + - ['tenia omentalis']*(elementsCountAroundTC//2)) + annotationArrayTC = (['taenia omentalis']*(elementsCountAroundTC//2) + + ['taenia libera']*elementsCountAroundTC + + ['taenia mesocolica']*elementsCountAroundTC + + ['taenia omentalis']*(elementsCountAroundTC//2)) annotationGroups += annotationGroupsTC - annotationArray += annotationArrayTC + annotationArrayAround += annotationArrayTC - return x, d1, d2, d3, annotationGroups, annotationArray + return x, d1, d2, d3, annotationGroups, annotationArrayAround def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, d3TC, nodesCountAroundTC, elementsCountAround, elementsCountAlong, @@ -1506,7 +1510,7 @@ def createNodesAndElementsTeniaColi(region, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, - annotationGroups, annotationArray, + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives): """ @@ -1523,7 +1527,8 @@ def createNodesAndElementsTeniaColi(region, :param elementsCountThroughWall: Number of elements through wall. :param tcCount: Number of tenia coli. :param annotationGroups: stores information about annotation groups. - :param annotationArray: stores annotation names of elements around. + :param annotationArrayAround: stores annotation names of elements around. + :param annotationArrayAlong: stores annotation names of elements along. :param firstNodeIdentifier, firstElementIdentifier: first node and element identifier to use. :param useCubicHermiteThroughWall: use linear when false @@ -1763,7 +1768,8 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 if annotationGroups: for annotationGroup in annotationGroups: - if annotationArray[e1] == annotationGroup._name: + if annotationArrayAround[e1] == annotationGroup._name or \ + annotationArrayAlong[e2] == annotationGroup._name: meshGroup = annotationGroup.getMeshGroup(mesh) meshGroup.addElement(element) @@ -1785,11 +1791,11 @@ def createNodesAndElementsTeniaColi(region, 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) elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - for annotationGroup in annotationGroups: - if annotationArray[elementsCountAround + eTC] == annotationGroup._name: - meshGroup = annotationGroup.getMeshGroup(mesh) - meshGroup.addElement(element) + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + eTC] == annotationGroup._name or \ + annotationArrayAlong[e2] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) for N in range(tcCount - 1): for eTC in range(elementsCountAroundTC): @@ -1826,12 +1832,12 @@ def createNodesAndElementsTeniaColi(region, element.merge(textureElementtemplate3) element.setNodesByIdentifier(eftTexture5, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - for annotationGroup in annotationGroups: - if annotationArray[elementsCountAround + int(elementsCountAroundTC*0.5) + - N*elementsCountAroundTC + eTC] == annotationGroup._name: - meshGroup = annotationGroup.getMeshGroup(mesh) - meshGroup.addElement(element) + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + int(elementsCountAroundTC*0.5) + + N*elementsCountAroundTC + eTC] == annotationGroup._name or \ + annotationArrayAlong[e2] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) for eTC in range(int(elementsCountAroundTC*0.5)): bni21 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 1 + tcOffset1 + \ @@ -1859,11 +1865,12 @@ def createNodesAndElementsTeniaColi(region, element.merge(textureElementtemplate5 if onOpening else textureElementtemplate4) element.setNodesByIdentifier(eftTexture7 if onOpening else eftTexture6, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - for annotationGroup in annotationGroups: - if annotationArray[elementsCountAround + int(elementsCountAroundTC*2.5) + eTC] == annotationGroup._name: - meshGroup = annotationGroup.getMeshGroup(mesh) - meshGroup.addElement(element) + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + int(elementsCountAroundTC*(tcCount - 0.5)) + eTC] \ + == annotationGroup._name \ + or annotationArrayAlong[e2] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) fm.endChange() diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index f0dcb491..7cd236fe 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -8,7 +8,6 @@ from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage -from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import tubemesh from scaffoldmaker.utils.tubemesh import CylindricalSegmentTubeMeshInnerPoints @@ -330,13 +329,14 @@ def generateBaseMesh(cls, region, options): # Create annotation groups annotationGroups = [] - annotationArray = [''] * (elementsCountAround) + annotationArrayAround = [''] * (elementsCountAround) + annotationArrayAlong = [''] * (elementsCountAlong) # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) return annotationGroups diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 5fbce2f7..372d2839 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -471,7 +471,7 @@ def createNodesAndElements(region, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArray, + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ @@ -486,7 +486,8 @@ def createNodesAndElements(region, :param elementsCountAlong: Number of elements along tube. :param elementsCountThroughWall: Number of elements through wall. :param annotationGroups: stores information about annotation groups. - :param annotationArray: stores annotation names of elements around. + :param annotationArrayAround: stores annotation names of elements around. + :param annotationArrayAlong: stores annotation names of elements along. :param firstNodeIdentifier, firstElementIdentifier: first node and element identifier to use. :param useCubicHermiteThroughWall: use linear when false @@ -700,7 +701,8 @@ def createNodesAndElements(region, elementIdentifier = elementIdentifier + 1 if annotationGroups: for annotationGroup in annotationGroups: - if annotationArray[e1] == annotationGroup._name: + if annotationArrayAround[e1] == annotationGroup._name or \ + annotationArrayAlong[e2] == annotationGroup._name: meshGroup = annotationGroup.getMeshGroup(mesh) meshGroup.addElement(element) diff --git a/tests/test_colon.py b/tests/test_colon.py index 5549f840..94c531ce 100644 --- a/tests/test_colon.py +++ b/tests/test_colon.py @@ -90,7 +90,7 @@ def test_colon1(self): del tmpRegion annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options) - self.assertEqual(3, len(annotationGroups)) + self.assertEqual(6, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) @@ -149,7 +149,7 @@ def test_mousecolon1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options) - self.assertEqual(2, len(annotationGroups)) + self.assertEqual(5, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) From 67dc6c14b5b8ca9f88dde7573147a5a49ee6b74c Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 17 Jun 2020 15:03:58 +1200 Subject: [PATCH 03/12] Add annotation groups for small intestine --- .../annotation/smallintestine_terms.py | 22 +++++++++++++++++++ .../meshtypes/meshtype_3d_smallintestine1.py | 22 ++++++++++++++----- 2 files changed, 39 insertions(+), 5 deletions(-) create mode 100644 src/scaffoldmaker/annotation/smallintestine_terms.py diff --git a/src/scaffoldmaker/annotation/smallintestine_terms.py b/src/scaffoldmaker/annotation/smallintestine_terms.py new file mode 100644 index 00000000..63e1d00b --- /dev/null +++ b/src/scaffoldmaker/annotation/smallintestine_terms.py @@ -0,0 +1,22 @@ +""" +Common resource for small intestine annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +smallintestine_terms = [ + ( "duodenum", "UBERON:0002114", "FMA:7206", "ILX:0726125"), + ( "ileum", "UBERON:0002116", "FMA:7208", "ILX:0728151"), + ( "jejunum", "UBERON:0002115", "FMA:7207", "ILX:0724224"), + ( "small intestine", "UBERON:0002108", "FMA:7200", "ILX:0726770") + ] + +def get_smallintestine_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in smallintestine_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Small intestine annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 7cd236fe..1a969158 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -5,6 +5,8 @@ """ import copy +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.smallintestine_terms import get_smallintestine_term from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -244,6 +246,7 @@ def generateBaseMesh(cls, region, options): # print(e+1, arcLength) length += arcLength segmentLength = length / segmentCount + elementAlongLength = length / elementsCountAlong # print('Length = ', length) # Sample central path @@ -256,6 +259,20 @@ def generateBaseMesh(cls, region, options): innerRadiusSegmentList, dInnerRadiusSegmentList = interp.sampleParameterAlongLine(lengthList, innerRadiusList, segmentCount) + # Create annotation groups for small intestine sections + elementsAlongDuodenum = round(duodenumLength / elementAlongLength) + elementsAlongJejunum = round(jejunumLength / elementAlongLength) + elementsAlongIleum = elementsCountAlong - elementsAlongDuodenum - elementsAlongJejunum + + duodenumGroup = AnnotationGroup(region, get_smallintestine_term("duodenum")) + jejunumGroup = AnnotationGroup(region, get_smallintestine_term("jejunum")) + ileumGroup = AnnotationGroup(region, get_smallintestine_term("ileum")) + annotationGroups = [duodenumGroup, jejunumGroup, ileumGroup] + annotationArrayAlong = (['duodenum'] * elementsAlongDuodenum + + ['jejunum'] * elementsAlongJejunum + + ['ileum'] * elementsAlongIleum) + annotationArrayAround = [''] * (elementsCountAround) + xExtrude = [] d1Extrude = [] d2Extrude = [] @@ -327,11 +344,6 @@ def generateBaseMesh(cls, region, options): xiList, flatWidthList, length, wallThickness, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) - # Create annotation groups - annotationGroups = [] - annotationArrayAround = [''] * (elementsCountAround) - annotationArrayAlong = [''] * (elementsCountAlong) - # Create nodes and elements nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, From 76cca022ce8f8494331893b24d81703e465d1d84 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 17 Jun 2020 17:16:31 +1200 Subject: [PATCH 04/12] Allow reduced haustrum appearance in transverse and distal pig colon --- .../meshtypes/meshtype_3d_cecum1.py | 4 +++- .../meshtypes/meshtype_3d_colon1.py | 16 ++++++++++++---- .../meshtypes/meshtype_3d_colonsegment1.py | 18 ++++++++++++------ 3 files changed, 27 insertions(+), 11 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index e3cb2087..0d9f49d9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -345,6 +345,8 @@ def generateBaseMesh(cls, region, options): [endTCWidth], [endTCWidthDerivative], xi)[0] tcWidthAlongCecum.append(tcWidth) + haustrumInnerRadiusFactorAlongCecum = [haustrumInnerRadiusFactor] * (elementsCountAlong + 1) + xToSample = [] d1ToSample = [] d2ToSample = [] @@ -355,7 +357,7 @@ def generateBaseMesh(cls, region, options): colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, - segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongCecum, innerRadiusAlongCecum, dInnerRadiusAlongCecum, tcWidthAlongCecum, startPhase) annotationArrayAlong = [''] * elementsCountAlongSegment diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 21deae24..3fc8ef6f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -236,11 +236,11 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Proximal inner radius'] = 38.0 options['Proximal tenia coli width'] = 5.0 options['Proximal-transverse inner radius'] = 14.0 - options['Proximal-transverse tenia coli width'] = 5.0 + options['Proximal-transverse tenia coli width'] = 4.0 options['Transverse-distal inner radius'] = 10.5 - options['Transverse-distal tenia coli width'] = 5.0 + options['Transverse-distal tenia coli width'] = 3.0 options['Distal inner radius'] = 8.0 - options['Distal tenia coli width'] = 5.0 + options['Distal tenia coli width'] = 1.5 elif 'Pig 2' in parameterSetName: options['Number of segments'] = 3 options['Proximal length'] = 30.0 @@ -426,6 +426,14 @@ def generateBaseMesh(cls, region, options): tcWidthList, elementsCountAlong) + # Account for reduced haustrum appearance in transverse and distal pig colon + if tcCount == 2: + haustrumInnerRadiusFactorList = [haustrumInnerRadiusFactor, haustrumInnerRadiusFactor*0.75, + haustrumInnerRadiusFactor*0.5, haustrumInnerRadiusFactor*0.2] + haustrumInnerRadiusFactorAlongElementList = \ + interp.sampleParameterAlongLine(lengthList, haustrumInnerRadiusFactorList, elementsCountAlong)[0] + else: + haustrumInnerRadiusFactorAlongElementList = [haustrumInnerRadiusFactor]*(elementsCountAlong+1) # Create annotation groups for colon sections elementsAlongInProximal = round(proximalLength/elementAlongLength) elementsAlongInTransverse = round(transverseLength/elementAlongLength) @@ -467,7 +475,7 @@ def generateBaseMesh(cls, region, options): colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, - segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongElementList, innerRadiusAlongElementList, dInnerRadiusAlongElementList, tcWidthAlongElementList, startPhase) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index bff51eab..3ebb272d 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -252,10 +252,12 @@ def generateBaseMesh(cls, region, options): [endTCWidth], [endTCWidthDerivative], xi)[0] tcWidthAlongSegment.append(tcWidth) + haustrumInnerRadiusFactorAlongSegment = [haustrumInnerRadiusFactor]*(elementsCountAlongSegment + 1) + colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, - segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongSegment, radiusAlongSegment, dRadiusAlongSegment, tcWidthAlongSegment, startPhase) annotationArrayAlong = [''] * elementsCountAlongSegment @@ -343,7 +345,7 @@ class ColonSegmentTubeMeshInnerPoints: def __init__(self, region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, segmentLength, wallThickness, - cornerInnerRadiusFactor, haustrumInnerRadiusFactor, innerRadiusAlongElementList, + cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongElementList, innerRadiusAlongElementList, dInnerRadiusAlongElementList, tcWidthAlongElementList, startPhase): self._region = region @@ -356,7 +358,7 @@ def __init__(self, region, elementsCountAroundTC, elementsCountAroundHaustrum, self._segmentLength = segmentLength self._wallThickness = wallThickness self._cornerInnerRadiusFactor = cornerInnerRadiusFactor - self._haustrumInnerRadiusFactor = haustrumInnerRadiusFactor + self._haustrumInnerRadiusFactorAlongElementList = haustrumInnerRadiusFactorAlongElementList self._innerRadiusAlongElementList = innerRadiusAlongElementList self._dInnerRadiusAlongElementList = dInnerRadiusAlongElementList self._tcWidthAlongElementList = tcWidthAlongElementList @@ -376,13 +378,16 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment): tcWidthSegmentList = self._tcWidthAlongElementList[nSegment*self._elementsCountAlongSegment: (nSegment + 1)*self._elementsCountAlongSegment + 1] + haustrumInnerRadiusFactorSegmentList = self._haustrumInnerRadiusFactorAlongElementList[nSegment * self._elementsCountAlongSegment: + (nSegment + 1) * self._elementsCountAlongSegment + 1] + xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, \ contractedWallThicknessSegment, segmentAxis, annotationGroups, annotationArray\ = getColonSegmentInnerPoints(self._region, self._elementsCountAroundTC, self._elementsCountAroundHaustrum, self._elementsCountAlongSegment, self._tcCount, self._segmentLengthEndDerivativeFactor, self._segmentLengthMidDerivativeFactor, self._segmentLength, self._wallThickness, - self._cornerInnerRadiusFactor, self._haustrumInnerRadiusFactor, + self._cornerInnerRadiusFactor, haustrumInnerRadiusFactorSegmentList, radiusSegmentList, dRadiusSegmentList, tcWidthSegmentList, self._startPhase) @@ -414,7 +419,7 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, segmentLength, wallThickness, - cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + cornerInnerRadiusFactor, haustrumInnerRadiusFactorSegmentList, radiusSegmentList, dRadiusSegmentList, tcWidthSegmentList, startPhase): """ @@ -440,7 +445,7 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, inter-haustral septa. Factor is multiplied by inner radius to get a radius of curvature at the corners. Only applicable for three tenia coli. Set to zero for two tenia coli. - :param haustrumInnerRadiusFactor: Factor is multiplied by inner + :param haustrumInnerRadiusFactorSegmentList: Factor is multiplied by inner radius to obtain radius of intersecting circles in the middle cross-section along a haustra segment. :param radiusSegmentList: List of inner radius defined from center of triangular @@ -551,6 +556,7 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, radius = radiusSegmentList[n2] sdRadius = dRadiusSegmentList[n2] tcWidth = tcWidthSegmentList[n2] + haustrumInnerRadiusFactor = haustrumInnerRadiusFactorSegmentList[n2] # Create segment of inner radius # Calculate x and d1 at the start, mid, and end faces From 137b752339f7c9aaec4a69f30eaac8a4114a8bd1 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 17 Jun 2020 17:18:31 +1200 Subject: [PATCH 05/12] Update unit test for small intestine --- tests/test_smallintestine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_smallintestine.py b/tests/test_smallintestine.py index 556a088f..ea1fff99 100644 --- a/tests/test_smallintestine.py +++ b/tests/test_smallintestine.py @@ -69,7 +69,7 @@ def test_smallintestine1(self): del tmpRegion annotationGroups = MeshType_3d_smallintestine1.generateBaseMesh(region, options) - self.assertEqual(0, len(annotationGroups)) + self.assertEqual(3, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) From 0176cd14b10088896869ffe3ca3fda2a71a953f4 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 2 Jul 2020 10:12:27 +1200 Subject: [PATCH 06/12] Correct how apex is handled in tubemesh --- src/scaffoldmaker/utils/tubemesh.py | 50 ++++++----------------------- 1 file changed, 9 insertions(+), 41 deletions(-) diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 372d2839..9185b9cb 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -135,9 +135,6 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, translateMatrix = [sx[nAlongSegment][j] - centroidRot[j] for j in range(3)] - if closedProximalEnd and nAlongSegment == 0: - translateMatrixForClosedEnd = copy.deepcopy(translateMatrix) - for n1 in range(elementsCountAround): x = xElementAlongSegment[n1] d1 = d1ElementAlongSegment[n1] @@ -171,10 +168,6 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, else: rotFrame2 = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] - # Store how second face is rotated to align sd2 to manipulate closed end later - if closedProximalEnd and nAlongSegment == 1: - rotFrame2ForClosedEnd = copy.deepcopy(rotFrame2) - xRot2 = [rotFrame2[j][0]*xRot1[0] + rotFrame2[j][1]*xRot1[1] + rotFrame2[j][2]*xRot1[2] for j in range(3)] d1Rot2 = [rotFrame2[j][0]*d1Rot1[0] + rotFrame2[j][1]*d1Rot1[1] + rotFrame2[j][2]*d1Rot1[2] for j in range(3)] d2Rot2 = [rotFrame2[j][0]*d2Rot1[0] + rotFrame2[j][1]*d2Rot1[1] + rotFrame2[j][2]*d2Rot1[2] for j in range(3)] @@ -184,34 +177,9 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, d1WarpedList.append(d1Rot2) d2WarpedList.append(d2Rot2) - # Manipulate closed end to align with how second face was rotated - if closedProximalEnd: - closedNodesNew = [] - closedD1New = [] - closedD2New = [] - for n1 in range(elementsCountAround): - closedNodesRot1 = [xWarpedList[n1][c] - translateMatrixForClosedEnd[c] for c in range(3)] - closedNodesRot2 = [rotFrame2ForClosedEnd[j][0] * closedNodesRot1[0] + - rotFrame2ForClosedEnd[j][1] * closedNodesRot1[1] + - rotFrame2ForClosedEnd[j][2] * closedNodesRot1[2] for j in range(3)] - closedD1Rot1 = d1WarpedList[n1] - closedD1Rot2 = [ rotFrame2ForClosedEnd[j][0] * closedD1Rot1[0] + - rotFrame2ForClosedEnd[j][1] * closedD1Rot1[1] + - rotFrame2ForClosedEnd[j][2] * closedD1Rot1[2] for j in range(3)] - closedD2Rot1 = d2WarpedList[n1] - closedD2Rot2 = [ rotFrame2ForClosedEnd[j][0] * closedD2Rot1[0] + - rotFrame2ForClosedEnd[j][1] * closedD2Rot1[1] + - rotFrame2ForClosedEnd[j][2] * closedD2Rot1[2] for j in range(3)] - closedNodesNew.append(closedNodesRot2) - closedD1New.append(closedD1Rot2) - closedD2New.append(closedD2Rot2) - xWarpedListNew = closedNodesNew + xWarpedList[elementsCountAround:] - d1WarpedListNew = closedD1New + d1WarpedList[elementsCountAround:] - d2WarpedListNew = closedD2New + d2WarpedList[elementsCountAround:] - else: - xWarpedListNew = xWarpedList - d1WarpedListNew = d1WarpedList - d2WarpedListNew = d2WarpedList + # xWarpedListNew = xWarpedList + # d1WarpedListNew = d1WarpedList + # d2WarpedListNew = d2WarpedList # Scale d2 with curvature of central path d2WarpedListScaled = [] @@ -221,7 +189,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, n = nAlongSegment * elementsCountAround + n1 # Calculate norm sd1Normalised = vector.normalise(sd1[nAlongSegment]) - v = [xWarpedListNew[n][c] - sx[nAlongSegment][c] for c in range(3)] + v = [xWarpedList[n][c] - sx[nAlongSegment][c] for c in range(3)] dp = vector.dotproduct(v, sd1Normalised) dpScaled = [dp * c for c in sd1Normalised] vProjected = [v[c] - dpScaled[c] for c in range(3)] @@ -245,7 +213,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, vProjectedNormlised, 0.0)) # Scale factor = 1.0 - curvature * innerRadiusAlong[nAlongSegment] - d2 = [factor * c for c in d2WarpedListNew[n]] + d2 = [factor * c for c in d2WarpedList[n]] d2WarpedListScaled.append(d2) # Smooth d2 for segment @@ -255,7 +223,7 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, nd2 = [] for n2 in range(elementsCountAlongSegment + 1): n = n2*elementsCountAround + n1 - nx.append(xWarpedListNew[n]) + nx.append(xWarpedList[n]) nd2.append(d2WarpedListScaled[n]) smoothd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixStartDerivative = True, fixEndDerivative = True) smoothd2Raw.append(smoothd2) @@ -266,12 +234,12 @@ def warpSegmentPoints(xList, d1List, d2List, segmentAxis, d2WarpedListFinal.append(smoothd2Raw[n1][n2]) # Calculate unit d3 - for n in range(len(xWarpedListNew)): - d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1WarpedListNew[n]), + for n in range(len(xWarpedList)): + d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1WarpedList[n]), vector.normalise(d2WarpedListFinal[n]))) d3WarpedUnitList.append(d3Unit) - return xWarpedListNew, d1WarpedListNew, d2WarpedListFinal, d3WarpedUnitList + return xWarpedList, d1WarpedList, d2WarpedListFinal, d3WarpedUnitList def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, wallThicknessList, elementsCountAround, From 94e1b6b2b7a00887f9288e205db53fdd1c0743ab Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 2 Jul 2020 10:36:27 +1200 Subject: [PATCH 07/12] Remove repetition for calculating d3 after warping --- src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 0d9f49d9..7a57a8af 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -401,18 +401,11 @@ def generateBaseMesh(cls, region, options): sd2ProjectedListRef, elementsCountAround, elementsCountAlong, zRefList, innerRadiusAlongCecum, closedProximalEnd=True) - # Calculate unit d3 - d3UnitWarpedList = [] - for n in range(len(xWarpedList)): - d3Unit = vector.normalise( - vector.crossproduct3(vector.normalise(d1WarpedList[n]), vector.normalise(d2WarpedList[n]))) - d3UnitWarpedList.append(d3Unit) - # Create coordinates and derivatives wallThicknessList = [wallThickness] * (elementsCountAlong + 1) xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, - d2WarpedList, d3UnitWarpedList, wallThicknessList, + d2WarpedList, d3WarpedUnitList, wallThicknessList, elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) # Deal with multiple nodes at end point for closed proximal end From 32bbbadbe109bbc60e0e8b851f1536a4c22a84d2 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 2 Jul 2020 11:51:34 +1200 Subject: [PATCH 08/12] Move cecum to align with origin before plane projection --- src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py | 9 +++++++++ src/scaffoldmaker/utils/tubemesh.py | 1 - 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 7a57a8af..8638f773 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -390,6 +390,15 @@ def generateBaseMesh(cls, region, options): elementsCountAroundHalfHaustrum, elementsCountAroundTC, elementsCountAround, elementsCountAlong, tcCount) + # Ensure cecum starts at z = 0.0 + minZ = xToWarp[0][2] + for n in range(len(xToWarp)): + if xToWarp[n][2] < minZ: + minZ = xToWarp[n][2] + + for n in range(len(xToWarp)): + xToWarp[n][2] = xToWarp[n][2] - minZ + # Project reference point for warping onto central path sxRefList, sd1RefList, sd2ProjectedListRef, zRefList = \ tubemesh.getPlaneProjectionOnCentralPath(xToWarp, elementsCountAround, elementsCountAlong, diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 9185b9cb..5d7342de 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -3,7 +3,6 @@ using a segment profile. ''' from __future__ import division -import copy import math from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldTextureCoordinates from opencmiss.zinc.element import Element From 93c3d58a04c6cc9b5c0c26f4f3f2756594ddbcf0 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 2 Jul 2020 12:17:35 +1200 Subject: [PATCH 09/12] Align with origin using min z-coordinate of the first node in each element group along instead --- src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 8638f773..eca5feb9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -392,9 +392,10 @@ def generateBaseMesh(cls, region, options): # Ensure cecum starts at z = 0.0 minZ = xToWarp[0][2] - for n in range(len(xToWarp)): - if xToWarp[n][2] < minZ: - minZ = xToWarp[n][2] + for n2 in range(elementsCountAlong + 1): + zFirstNodeAlong = xToWarp[n2 * elementsCountAround][2] + if zFirstNodeAlong < minZ: + minZ = zFirstNodeAlong for n in range(len(xToWarp)): xToWarp[n][2] = xToWarp[n][2] - minZ From 4d01a2ef91bcaa6bfb3cd90885c769ac3094a781 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Jul 2020 12:06:47 +1200 Subject: [PATCH 10/12] Add tenia coli for cecum --- .../meshtypes/meshtype_3d_cecum1.py | 86 ++- .../meshtypes/meshtype_3d_colon1.py | 9 +- .../meshtypes/meshtype_3d_colonsegment1.py | 686 ++++++++++++------ .../utils/eftfactory_bicubichermitelinear.py | 163 +++++ .../utils/eftfactory_tricubichermite.py | 219 +++++- 5 files changed, 919 insertions(+), 244 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index eca5feb9..24c632a1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -7,7 +7,8 @@ import copy import math from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion -from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import ColonSegmentTubeMeshInnerPoints, getFullProfileFromHalfHaustrum +from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import ColonSegmentTubeMeshInnerPoints, \ + getFullProfileFromHalfHaustrum, getTeniaColi, createNodesAndElementsTeniaColi from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1, generateOstiumMesh from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -320,7 +321,8 @@ def generateBaseMesh(cls, region, options): cecumLength += arcLength # Sample central path - smoothd1 = interp.smoothCubicHermiteDerivativesLine(cx, cd1, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + smoothd1 = interp.smoothCubicHermiteDerivativesLine(cx, cd1, + magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) sxCecum, sd1Cecum, se, sxi, ssf = interp.sampleCubicHermiteCurves(cx, smoothd1, elementsCountAlong) sd2Cecum, sd12Cecum = interp.interpolateSampleCubicHermite(cd2, cd12, se, sxi, ssf) @@ -332,6 +334,8 @@ def generateBaseMesh(cls, region, options): dInnerRadiusAlongCecum = [] tcWidthAlongCecum = [] + closedProximalEnd = True + for n2 in range(elementsCountAlongSegment*segmentCount + 1): xi = 1/(elementsCountAlongSegment*segmentCount) * n2 @@ -360,7 +364,7 @@ def generateBaseMesh(cls, region, options): segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongCecum, innerRadiusAlongCecum, dInnerRadiusAlongCecum, tcWidthAlongCecum, startPhase) - annotationArrayAlong = [''] * elementsCountAlongSegment + annotationArrayAlong = [''] * elementsCountAlongSegment * segmentCount for nSegment in range(segmentCount): # Make regular segments @@ -409,7 +413,7 @@ def generateBaseMesh(cls, region, options): xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = \ tubemesh.warpSegmentPoints(xToWarp, d1ToWarp, d2ToWarp, segmentAxis, sxRefList, sd1RefList, sd2ProjectedListRef, elementsCountAround, elementsCountAlong, - zRefList, innerRadiusAlongCecum, closedProximalEnd=True) + zRefList, innerRadiusAlongCecum, closedProximalEnd) # Create coordinates and derivatives wallThicknessList = [wallThickness] * (elementsCountAlong + 1) @@ -421,11 +425,13 @@ def generateBaseMesh(cls, region, options): # Deal with multiple nodes at end point for closed proximal end xApexInner = xList[0] # arclength between apex point and corresponding point on next face - mag = interp.getCubicHermiteArcLength(xList[0], d2List[0], xList[elementsCountAround*2], d2List[elementsCountAround*2]) + mag = interp.getCubicHermiteArcLength(xList[0], d2List[0], xList[elementsCountAround*2], + d2List[elementsCountAround*2]) d2ApexInner = vector.setMagnitude(sd2Cecum[0], mag) d1ApexInner = vector.crossproduct3(sd1Cecum[0], d2ApexInner) d1ApexInner = vector.setMagnitude(d1ApexInner, mag) - d3ApexUnit = vector.normalise(vector.crossproduct3(vector.normalise(d1ApexInner), vector.normalise(d2ApexInner))) + d3ApexUnit = vector.normalise(vector.crossproduct3(vector.normalise(d1ApexInner), + vector.normalise(d2ApexInner))) d3ApexInner = [d3ApexUnit[c] * wallThickness/elementsCountThroughWall for c in range(3)] xCecum = [] @@ -449,11 +455,26 @@ def generateBaseMesh(cls, region, options): xTexture = d1Texture = d2Texture = d3Texture = [] # Create nodes and elements - nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( - region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, - elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=True) + if tcThickness > 0: + tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() + xCecum, d1Cecum, d2Cecum, d3Cecum, annotationGroups, annotationArrayAround = getTeniaColi( + region, xCecum, d1Cecum, d2Cecum, d3Cecum, curvatureList, tcCount, elementsCountAroundTC, + elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, + tubeTCWidthList, tcThickness, sxRefList, annotationGroups, annotationArrayAround, + closedProximalEnd) + + nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + + else: + nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( + region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, + annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) # Add ostium on track surface between two tenia on the last segment elementsAroundTrackSurface = elementsCountAroundHaustrum @@ -464,14 +485,17 @@ def generateBaseMesh(cls, region, options): startIdxElementsAround = int((elementsCountAroundHaustrum + elementsCountAroundTC)*sectorIdx + elementsCountAroundTC*0.5) baseNodesIdx = (elementsCountThroughWall + 1) + \ - + elementsCountAround * (elementsCountThroughWall + 1) * (elementsCountAlongSegment * (segmentCount - 1) - 1) \ - + elementsCountAround + + (elementsCountAround * (elementsCountThroughWall + 1) + + ((elementsCountAroundTC - 1)*tcCount if tcThickness > 0.0 else 0)) * \ + (elementsCountAlongSegment * (segmentCount - 1) - 1) + elementsCountAround xTrackSurface = [] d1TrackSurface = [] d2TrackSurface = [] for n2 in range(elementsCountAlongSegment + 1): for n1 in range(elementsCountAroundHaustrum + 1): - idx = baseNodesIdx + elementsCountAround * (elementsCountThroughWall + 1) * n2 + \ + idx = baseNodesIdx + \ + (elementsCountAround * (elementsCountThroughWall + 1) + + ((elementsCountAroundTC - 1)*tcCount if tcThickness > 0.0 else 0)) * n2 + \ startIdxElementsAround + n1 xTrackSurface.append(xCecum[idx]) d1TrackSurface.append(d1Cecum[idx]) @@ -530,20 +554,26 @@ def generateBaseMesh(cls, region, options): ei2Bottom >= 0 and ei2Top < elementsAlongTrackSurface), \ 'cecum1.py: Insufficient elements around ostium on tracksurface to make annulus mesh.' - nodeStart = int(baseNodesIdx + elementsCountAround * (elementsCountThroughWall + 1) * ei2Bottom + ei1Centre + + nodeStart = int(baseNodesIdx + ei2Bottom * (elementsCountAround * (elementsCountThroughWall + 1) + \ + ((elementsCountAroundTC - 1)*tcCount if tcThickness > 0.0 else 0)) + ei1Centre + \ sectorIdx*(elementsCountAroundHaustrum + elementsCountAroundTC) + elementsCountAroundTC*0.5) - \ elementsCountAround + 1 # only for 1 layer through wall # Store elements and nodes to be deleted later from tracked surface deleteElementsCountAcross = ei1Right - ei1Left + 1 deleteElementsCountAlong = ei2Top - ei2Bottom + 1 - deleteElementIdxStart = int(elementsCountAround * (elementsCountAlong - elementsCountAlongSegment + ei2Bottom) - *elementsCountThroughWall + elementsCountAroundTC*0.5 + - (elementsCountAroundTC + elementsCountAroundHaustrum)*sectorIdx + ei1Left + 1) + deleteElementIdxStart = int((elementsCountAround * elementsCountThroughWall + + (elementsCountAroundTC * tcCount if tcThickness > 0.0 else 0.0)) * + (elementsCountAlong - elementsCountAlongSegment + ei2Bottom) + + elementsCountAroundTC * 0.5 + + (elementsCountAroundTC + elementsCountAroundHaustrum) * sectorIdx + ei1Left + 1) + deleteElementIdentifier = [] for n2 in range(deleteElementsCountAlong): for n1 in range(deleteElementsCountAcross): - elementIdx = deleteElementIdxStart + n1 + elementsCountAround*n2 + elementIdx = deleteElementIdxStart + n1 + n2 * (elementsCountAround + + (int(elementsCountAroundTC * tcCount) + if tcThickness > 0.0 else 0)) deleteElementIdentifier.append(elementIdx) deleteNodeIdxStart = nodeStart - int(deleteElementsCountAcross * 0.5) @@ -551,8 +581,9 @@ def generateBaseMesh(cls, region, options): for n2 in range(deleteElementsCountAlong + 1): for n3 in range(elementsCountThroughWall + 1): for n1 in range(deleteElementsCountAcross + 1): - nodeIdx = deleteNodeIdxStart + n1 + elementsCountAround*n3 + \ - elementsCountAround*(elementsCountThroughWall+1)*n2 + nodeIdx = deleteNodeIdxStart + n1 + elementsCountAround * n3 + \ + n2 * (elementsCountAround * (elementsCountThroughWall + 1) + + ((elementsCountAroundTC - 1) * tcCount if tcThickness > 0.0 else 0)) deleteNodeIdentifier.append(nodeIdx) innerEndPoints_Id = [] @@ -562,7 +593,7 @@ def generateBaseMesh(cls, region, options): innerEndPoints_Id.append(idx) endProportions.append([(ei1Centre - n1)/elementsAroundTrackSurface, ei2Bottom/elementsAlongTrackSurface]) for n2 in range(ei2Top - ei2Bottom + 1): - idx = idx + 2*elementsCountAround + idx = idx + 2 * elementsCountAround + ((elementsCountAroundTC - 1) * tcCount if tcThickness > 0.0 else 0) innerEndPoints_Id.append(idx) endProportions.append([ei1Left/elementsAroundTrackSurface, (ei2Bottom + n2 + 1)/elementsAlongTrackSurface]) for n1 in range(1, ei1Right - ei1Left + 2): @@ -570,7 +601,7 @@ def generateBaseMesh(cls, region, options): innerEndPoints_Id.append(idx) endProportions.append([(ei1Left + n1) / elementsAroundTrackSurface, (ei2Top+1) / elementsAlongTrackSurface]) for n2 in range(ei2Top - ei2Bottom + 1): - idx = idx - 2*elementsCountAround + idx = idx - 2 * elementsCountAround - ((elementsCountAroundTC - 1) * tcCount if tcThickness > 0.0 else 0) innerEndPoints_Id.append(idx) endProportions.append( [(ei1Right+1) / elementsAroundTrackSurface, (ei2Top - n2) / elementsAlongTrackSurface]) @@ -697,12 +728,13 @@ def getApexSegmentForCecum(xInner, d1Inner, d2Inner, elementsCountAroundHalfHaus # Compile nodes and d2 for sampling xFirstSegment += xInner[elementsCountAround * int(elementsCountAlongSegment * 0.5):] # second half of first regular segment d1FirstDirectionVector = vector.normalise(d1Inner[elementsCountAround]) # Store direction vector of first d1 intra-haustral for later - d2Vector = xInner[elementsCountAround*int(elementsCountAlongSegment*0.5):elementsCountAround*(int(elementsCountAlongSegment*0.5)+1)] # half face of segment - apex + d2Vector = xInner[elementsCountAround * int(elementsCountAlongSegment * 0.5): + elementsCountAround * (int(elementsCountAlongSegment * 0.5) + 1)] # half face of segment - apex d2FirstSegment = [] for c in range(elementsCountAround): d2 = [d2Vector[c][0], d2Vector[c][1], 0.0 ] # project onto x-y plane to get d2 pointing vertically d2FirstSegment.append(d2) - d2FirstSegment += d2Inner[elementsCountAround*int(elementsCountAlongSegment*0.5):] + d2FirstSegment += d2Inner[elementsCountAround * int(elementsCountAlongSegment*0.5):] # Sample along first segment xFirstSegmentSampledRaw = [] @@ -714,8 +746,8 @@ def getApexSegmentForCecum(xInner, d1Inner, d2Inner, elementsCountAroundHalfHaus for n1 in range(elementsCountAround): xForSamplingAlong = [] d2ForSamplingAlong = [] - for n2 in range(1 + elementsCountAlongSegment - int(elementsCountAlongSegment*0.5) + 1): - idx = elementsCountAround*n2 + n1 + for n2 in range(1 + elementsCountAlongSegment - int(elementsCountAlongSegment * 0.5) + 1): + idx = elementsCountAround * n2 + n1 xForSamplingAlong.append(xFirstSegment[idx]) d2ForSamplingAlong.append(d2FirstSegment[idx]) xResampled, d1Resampled, se, sxi, _ = interp.sampleCubicHermiteCurves(xForSamplingAlong, d2ForSamplingAlong, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 3fc8ef6f..f1f72762 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -516,25 +516,28 @@ def generateBaseMesh(cls, region, options): relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() + closedProximalEnd = False + if tcThickness > 0: tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi( region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroups, annotationArrayAround) + tubeTCWidthList, tcThickness, sxRefExtrudeList, annotationGroups, annotationArrayAround, + closedProximalEnd) # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( xiList, relaxedLengthList, length, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, - elementsCountThroughWall, transitElementList) + 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, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives) + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) else: # Create flat and texture coordinates diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 3ebb272d..11865541 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -288,25 +288,27 @@ def generateBaseMesh(cls, region, options): relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() + closedProximalEnd = False + if tcThickness > 0: tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi( region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, - tubeTCWidthList, tcThickness, sxRefList, annotationGroups, annotationArrayAround) + tubeTCWidthList, tcThickness, sxRefList, annotationGroups, annotationArrayAround, closedProximalEnd) # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( xiList, relaxedLengthList, segmentLength, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, - elementsCountThroughWall, transitElementList) + 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, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives) + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) else: # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( @@ -318,7 +320,7 @@ def generateBaseMesh(cls, region, options): region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) + firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) return annotationGroups @@ -378,8 +380,9 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment): tcWidthSegmentList = self._tcWidthAlongElementList[nSegment*self._elementsCountAlongSegment: (nSegment + 1)*self._elementsCountAlongSegment + 1] - haustrumInnerRadiusFactorSegmentList = self._haustrumInnerRadiusFactorAlongElementList[nSegment * self._elementsCountAlongSegment: - (nSegment + 1) * self._elementsCountAlongSegment + 1] + haustrumInnerRadiusFactorSegmentList = self._haustrumInnerRadiusFactorAlongElementList[ + nSegment * self._elementsCountAlongSegment: + (nSegment + 1) * self._elementsCountAlongSegment + 1] xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, \ contractedWallThicknessSegment, segmentAxis, annotationGroups, annotationArray\ @@ -1207,7 +1210,7 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tubeTCWidthList, tcThickness, sx, - annotationGroups, annotationArrayAround): + annotationGroups, annotationArrayAround, closedProximalEnd): """ Create equally spaced points for tenia coli over the outer surface of the colon. Points are sampled from a cubic @@ -1230,12 +1233,13 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, :param tubeTCWidthList: List of tenia coli width along tube length. :param tcThickness: Thickness of tenia coli at its thickest part. :param sx: Coordinates of central path. - :param annotationGroups, annotationArrayAround: annotations for tube + :param annotationGroups, annotationArrayAround: annotations for tube. + :param closedProximalEnd: True when proximal end of tube is closed. :return: coordinates, derivatives, annotationGroups for colon with tenia coli. """ - elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum) * tcCount TCEdgeFactor = 1.5 xTCArranged = [] d1TCArranged = [] @@ -1243,22 +1247,28 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, d3TCArranged = [] # Create points for tenia coli - for n2 in range(elementsCountAlong + 1): + for n2 in range(1 if closedProximalEnd else 0, elementsCountAlong + 1): xTCRaw = [] d1TCRaw = [] d2TCRaw = [] d3TCRaw = [] tcWidth = tubeTCWidthList[n2] for N in range(tcCount): - idxTCMid = n2*elementsCountAround*(elementsCountThroughWall + 1) + \ - elementsCountAround*elementsCountThroughWall + \ - N*(elementsCountAroundTC + elementsCountAroundHaustrum) + idxTCMid = elementsCountThroughWall + 1 + (n2 - 1) * elementsCountAround * (elementsCountThroughWall + 1) \ + + elementsCountAround * elementsCountThroughWall + \ + N * (elementsCountAroundTC + elementsCountAroundHaustrum) if closedProximalEnd \ + else n2 * elementsCountAround * (elementsCountThroughWall + 1) + \ + elementsCountAround * elementsCountThroughWall + \ + N * (elementsCountAroundTC + elementsCountAroundHaustrum) unitNorm = vector.normalise(d3List[idxTCMid]) xMid = [xList[idxTCMid][i] + unitNorm[i]*tcThickness for i in range(3)] d1Mid = d1List[idxTCMid] TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) if N > 0 else \ idxTCMid + tcCount*(elementsCountAroundTC + elementsCountAroundHaustrum) - int(elementsCountAroundTC*0.5) TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) + if closedProximalEnd: + tcWidth = vector.magnitude([xList[TCStartIdx][c] - xList[TCEndIdx][c] for c in range(3)]) + v1 = xList[TCStartIdx] v2 = xMid d1MidScaled = [c*tcWidth*TCEdgeFactor for c in vector.normalise(d1Mid)] @@ -1301,14 +1311,14 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, x, d1, d2, d3 = combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTCArranged, d1TCArranged, d2TCArranged, d3TCArranged, (elementsCountAroundTC - 1)*tcCount, elementsCountAround, - elementsCountAlong, elementsCountThroughWall) + elementsCountAlong, elementsCountThroughWall, closedProximalEnd) # Update annotation groups - if tcCount == 2: + if tcCount == 2 or closedProximalEnd: tcGroup = AnnotationGroup(region, get_colon_term("taenia coli")) annotationGroups += [tcGroup] annotationArrayAround += (['taenia coli']*elementsCountAroundTC*tcCount) - if tcCount == 3: + elif tcCount == 3: tlGroup = AnnotationGroup(region, get_colon_term("taenia libera")) tmGroup = AnnotationGroup(region, get_colon_term("taenia mesocolica")) toGroup = AnnotationGroup(region, get_colon_term("taenia omentalis")) @@ -1324,7 +1334,7 @@ def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, d3TC, nodesCountAroundTC, elementsCountAround, elementsCountAlong, - elementsCountThroughWall): + elementsCountThroughWall, closedProximalEnd): """ Arranges coordinates and derivatives around inner surface to outer surface, followed by tenia coli points before extending @@ -1335,6 +1345,7 @@ def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, :param elementsCountAround: Number of elements around colon. :param elementsCountAlong: Number of elements along colon. :param elementsCountThroughWall: Number of elements through wall. + :param closedProximalEnd: True when proximal end of tube is closed. : return: reordered coordinates and derivatives """ x = [] @@ -1345,30 +1356,42 @@ def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, # Add tenia coli points to coordinates list for n2 in range(elementsCountAlong + 1): for n3 in range(elementsCountThroughWall + 1): - for n1 in range(elementsCountAround): - # Append colon wall coordinates from inside to outside wall - n = n2*elementsCountAround*(elementsCountThroughWall + 1) + n3*elementsCountAround + n1 - x.append(xList[n]) - d1.append(d1List[n]) - d2.append(d2List[n]) + if closedProximalEnd and n2 == 0: + x.append(xList[n3]) + d1.append(d1List[n3]) + d2.append(d2List[n3]) if d3List: - d3.append(d3List[n]) + d3.append(d3List[n3]) + else: + for n1 in range(elementsCountAround): + # Append colon wall coordinates from inside to outside wall + n = (elementsCountThroughWall + 1) + \ + (n2 - 1) * elementsCountAround * (elementsCountThroughWall + 1) + \ + n3 * elementsCountAround + n1 if closedProximalEnd \ + else n2 * elementsCountAround * (elementsCountThroughWall + 1) + n3 * elementsCountAround + n1 + + x.append(xList[n]) + d1.append(d1List[n]) + d2.append(d2List[n]) + if d3List: + d3.append(d3List[n]) # Append tenia coli coordinates - for nTC in range(nodesCountAroundTC): - nTCCount = n2*nodesCountAroundTC + nTC - x.append(xTC[nTCCount]) - d1.append(d1TC[nTCCount]) - d2.append(d2TC[nTCCount]) - if d3TC: - d3.append(d3TC[nTCCount]) + if not (closedProximalEnd and n2 == 0): + for nTC in range(nodesCountAroundTC): + nTCCount = (n2 - 1 if closedProximalEnd else n2) * nodesCountAroundTC + nTC + x.append(xTC[nTCCount]) + d1.append(d1TC[nTCCount]) + d2.append(d2TC[nTCCount]) + if d3TC: + d3.append(d3TC[nTCCount]) return x, d1, d2, d3 def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, totalLengthAlong, wallThickness, tcCount, tcThickness, elementsCountAroundTC, elementsCountAroundHaustrum, - elementsCountAlong, elementsCountThroughWall, transitElementList): + elementsCountAlong, elementsCountThroughWall, transitElementList, closedProximalEnd): """ Calculates flat coordinates and texture coordinates for a colon scaffold with tenia coli when it is opened @@ -1387,6 +1410,7 @@ def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, :param elementsCountThroughWall: Number of elements through wall. :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. """ @@ -1501,12 +1525,12 @@ def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, xFlat, d1Flat, d2Flat, _ = combineTeniaColiWithColon(xFlatColon, d1FlatColon, d2FlatColon, [], xFlatListTC, d1FlatListTC, d2FlatListTC, [], (elementsCountAroundTC - 1)*tcCount + 1, - elementsCountAround + 1, elementsCountAlong, elementsCountThroughWall) + elementsCountAround + 1, elementsCountAlong, elementsCountThroughWall, closedProximalEnd) xTexture, d1Texture, d2Texture, _ = combineTeniaColiWithColon(xTextureColon, d1TextureColon, d2TextureColon, [], xTextureListTC, d1TextureListTC, d2TextureListTC, [], (elementsCountAroundTC - 1)*tcCount + 1, elementsCountAround + 1, elementsCountAlong, - elementsCountThroughWall) + elementsCountThroughWall, closedProximalEnd) return xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture @@ -1518,10 +1542,11 @@ def createNodesAndElementsTeniaColi(region, elementsCountAlong, elementsCountThroughWall, tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives): + useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ Create nodes and elements for the coordinates, flat coordinates, - and texture coordinates field. + and texture coordinates field. Note that flat and texture 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. @@ -1537,8 +1562,9 @@ def createNodesAndElementsTeniaColi(region, :param annotationArrayAlong: stores annotation names of elements along. :param firstNodeIdentifier, firstElementIdentifier: first node and element identifier to use. - :param useCubicHermiteThroughWall: use linear when false - :param useCrossDerivatives: use cross derivatives when true + :param useCubicHermiteThroughWall: use linear when false. + :param useCrossDerivatives: use cross derivatives when true. + :param closedProximalEnd: True when proximal end of tube is closed. :return nodeIdentifier, elementIdentifier, annotationGroups """ @@ -1591,24 +1617,6 @@ def createNodesAndElementsTeniaColi(region, eft2 = eftfactory.createEftWedgeXi1Zero() elementtemplate2.defineField(coordinates, -1, eft2) - # Create flat coordinates field - flatCoordinates = findOrCreateFieldCoordinates(fm, name="flat coordinates") - flatNodetemplate1 = nodes.createNodetemplate() - flatNodetemplate1.defineField(flatCoordinates) - flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) - flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - - flatNodetemplate2 = nodes.createNodetemplate() - flatNodetemplate2.defineField(flatCoordinates) - flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) - flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) - flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) - if useCrossDerivatives: - flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) - bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) eftTexture3 = bicubichermitelinear.createEftBasic() eftTexture4 = bicubichermitelinear.createEftOpenTube() @@ -1616,63 +1624,83 @@ def createNodesAndElementsTeniaColi(region, eftTexture6 = bicubichermitelinear.createEftWedgeXi1Zero() eftTexture7 = bicubichermitelinear.createEftWedgeXi1ZeroOpenTube() - flatElementtemplate1 = mesh.createElementtemplate() - flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate1.defineField(flatCoordinates, -1, eftTexture3) - - flatElementtemplate2 = mesh.createElementtemplate() - flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture4) - - flatElementtemplate3 = mesh.createElementtemplate() - flatElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate3.defineField(flatCoordinates, -1, eftTexture5) - - flatElementtemplate4 = mesh.createElementtemplate() - flatElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate4.defineField(flatCoordinates, -1, eftTexture6) - - flatElementtemplate5 = mesh.createElementtemplate() - flatElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate5.defineField(flatCoordinates, -1, eftTexture7) - - # 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) + if xFlat: + # Create flat coordinates field + flatCoordinates = findOrCreateFieldCoordinates(fm, name="flat coordinates") + flatNodetemplate1 = nodes.createNodetemplate() + flatNodetemplate1.defineField(flatCoordinates) + flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) + flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + flatNodetemplate1.setValueNumberOfVersions(flatCoordinates, -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) + flatNodetemplate2 = nodes.createNodetemplate() + flatNodetemplate2.defineField(flatCoordinates) + flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) + flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) + flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) + if useCrossDerivatives: + flatNodetemplate2.setValueNumberOfVersions(flatCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) + + flatElementtemplate1 = mesh.createElementtemplate() + flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + flatElementtemplate1.defineField(flatCoordinates, -1, eftTexture3) + + flatElementtemplate2 = mesh.createElementtemplate() + flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture4) + + flatElementtemplate3 = mesh.createElementtemplate() + flatElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) + flatElementtemplate3.defineField(flatCoordinates, -1, eftTexture5) + + flatElementtemplate4 = mesh.createElementtemplate() + flatElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) + flatElementtemplate4.defineField(flatCoordinates, -1, eftTexture6) + + 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) - textureElementtemplate1 = mesh.createElementtemplate() - textureElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - textureElementtemplate1.defineField(textureCoordinates, -1, eftTexture3) + 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) + 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) + 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) + 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) + textureElementtemplate5 = mesh.createElementtemplate() + textureElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate5.defineField(textureCoordinates, -1, eftTexture7) # create nodes for coordinates field for n in range(len(x)): @@ -1691,86 +1719,277 @@ def createNodesAndElementsTeniaColi(region, nodeIdentifier = nodeIdentifier + 1 # Create nodes for flat coordinates field - nodeIdentifier = firstNodeIdentifier - for n2 in range(elementsCountAlong + 1): - for n3 in range(elementsCountThroughWall + 1): - for n1 in range(elementsCountAround): - i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 + \ - n2*((elementsCountAroundTC - 1)*tcCount + 1) + if xFlat and xTexture: + nodeIdentifier = firstNodeIdentifier + for n2 in range(elementsCountAlong + 1): + for n3 in range(elementsCountThroughWall + 1): + for n1 in range(elementsCountAround): + i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + \ + (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 + for nTC in range((elementsCountAroundTC - 1)*tcCount): + j = i + 2 + nTC node = nodes.findNodeByIdentifier(nodeIdentifier) - node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) - node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) + node.merge(flatNodetemplate2 if nTC == 0 else flatNodetemplate1) + node.merge(textureNodetemplate2 if nTC == 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]) + 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 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 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 flat coordinates nodes for tenia coli - for nTC in range((elementsCountAroundTC - 1)*tcCount): - 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 - now = elementsCountAround*(elementsCountThroughWall+1) - tcOffset = (elementsCountAroundTC-1)*tcCount - for e2 in range(elementsCountAlong): - tcOffset1 = e2*(elementsCountAroundTC-1)*tcCount + if closedProximalEnd: + elementtemplate3 = mesh.createElementtemplate() + elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) + radiansPerElementAround = math.pi * 2.0 / elementsCountAround + + elementtemplate4 = mesh.createElementtemplate() + elementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + elementtemplate5 = mesh.createElementtemplate() + elementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + elementtemplate6 = mesh.createElementtemplate() + elementtemplate6.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + # Create apex for e3 in range(elementsCountThroughWall): for e1 in range(elementsCountAround): - bni11 = e2*now + e3*elementsCountAround + e1 + 1 + tcOffset1 - bni12 = e2*now + e3*elementsCountAround + (e1 + 1) % elementsCountAround + 1 + tcOffset1 - bni21 = e2*now + (e3 + 1)*elementsCountAround + e1 + 1 + tcOffset1 - bni22 = e2*now + (e3 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1 + tcOffset1 - nodeIdentifiers = [ bni11, bni12, bni11 + now + tcOffset, bni12 + now + tcOffset, - bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset] + va = e1 + vb = (e1 + 1) % elementsCountAround + eft3 = eftfactory.createEftShellPoleBottom(va * 100, vb * 100) + elementtemplate3.defineField(coordinates, -1, eft3) + element = mesh.createElement(elementIdentifier, elementtemplate3) + bni1 = e3 + 1 + bni2 = elementsCountThroughWall + 1 + elementsCountAround * e3 + e1 + 1 + bni3 = elementsCountThroughWall + 1 + elementsCountAround * e3 + (e1 + 1) % elementsCountAround + 1 + nodeIdentifiers = [bni1, bni2, bni3, bni1 + 1, bni2 + elementsCountAround, bni3 + elementsCountAround] + element.setNodesByIdentifier(eft3, nodeIdentifiers) + # set general linear map coefficients + radiansAround = e1 * radiansPerElementAround + radiansAroundNext = ((e1 + 1) % elementsCountAround) * radiansPerElementAround + scalefactors = [ + -1.0, + math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, + math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround, + math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, + math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround + ] + result = element.setScaleFactors(eft3, scalefactors) + elementIdentifier = elementIdentifier + 1 + for annotationGroup in annotationGroups: + if annotationArrayAround[e1] == annotationGroup._name or \ + annotationArrayAlong[0] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + # Create tenia coli + for eTC in range(int(elementsCountAroundTC * 0.5)): + va = eTC + vb = eTC + 1 + # set general linear map coefficients + radiansAround = va * radiansPerElementAround + radiansAroundNext = vb * radiansPerElementAround + scalefactors = [-1.0, + math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, + math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround] + + bni21 = elementsCountThroughWall + 1 + bni22 = bni21 + elementsCountAround * elementsCountThroughWall + 1 + eTC + bni23 = bni22 + 1 + bni31 = bni22 + elementsCountAround + + tetrahedralElement = (eTC == int(elementsCountAroundTC * 0.5) - 1) + if tetrahedralElement: + eft4 = eftfactory.createEftTetrahedronXi1One(va * 10000, vb * 10000) + elementtemplate4.defineField(coordinates, -1, eft4) + nodeIdentifiers = [bni21, bni22, bni23, bni31] + else: + eft6 = eftfactory.createEftPyramidBottomSimple(va * 10000, vb * 10000) + elementtemplate6.defineField(coordinates, -1, eft6) + nodeIdentifiers = [bni21, bni22, bni23, bni31, bni31 + 1] + + element = mesh.createElement(elementIdentifier, elementtemplate4 if tetrahedralElement else elementtemplate6) + element.setNodesByIdentifier(eft4 if tetrahedralElement else eft6, nodeIdentifiers) + element.setScaleFactors(eft4 if tetrahedralElement else eft6, scalefactors) + elementIdentifier = elementIdentifier + 1 + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + eTC] == annotationGroup._name or \ + annotationArrayAlong[0] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + for N in range(tcCount - 1): + for eTC in range(elementsCountAroundTC): + va = int(elementsCountAroundTC * 0.5) + elementsCountAroundHaustrum + \ + N * (elementsCountAroundHaustrum + elementsCountAroundTC) + eTC + vb = va + 1 + # set general linear map coefficients + radiansAround = va * radiansPerElementAround + radiansAroundNext = vb * radiansPerElementAround + scalefactors = [ + -1.0, + math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, + math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround] + + bni21 = elementsCountThroughWall + 1 + bni22 = bni21 + elementsCountAround * elementsCountThroughWall + int(elementsCountAroundTC * 0.5) + \ + elementsCountAroundHaustrum + N*(elementsCountAroundTC + elementsCountAroundHaustrum) + eTC + 1 + bni23 = bni22 + 1 + bni31 = bni21 + elementsCountAround*(elementsCountThroughWall + 1) + int(elementsCountAroundTC*0.5) + \ + N*(elementsCountAroundTC - 1) + eTC + 1 + + if eTC == 0: + eft5 = eftfactory.createEftTetrahedronXi1Zero(va * 10000, vb * 10000) + elementtemplate5.defineField(coordinates, -1, eft5) + nodeIdentifiers = [bni21, bni22, bni23, bni31] + element = mesh.createElement(elementIdentifier, elementtemplate5) + element.setNodesByIdentifier(eft5, nodeIdentifiers) + element.setScaleFactors(eft5, scalefactors) + + elif eTC > 0 and eTC < elementsCountAroundTC - 1: + eft6 = eftfactory.createEftPyramidBottomSimple(va * 10000, vb * 10000) + elementtemplate6.defineField(coordinates, -1, eft6) + nodeIdentifiers = [bni21, bni22, bni23, bni31-1, bni31] + element = mesh.createElement(elementIdentifier, elementtemplate6) + element.setNodesByIdentifier(eft6, nodeIdentifiers) + element.setScaleFactors(eft6, scalefactors) + + else: + eft4 = eftfactory.createEftTetrahedronXi1One(va * 10000, vb * 10000) + elementtemplate4.defineField(coordinates, -1, eft4) + nodeIdentifiers = [bni21, bni22, bni23, bni31 - 1] + element = mesh.createElement(elementIdentifier, elementtemplate4) + element.setNodesByIdentifier(eft4, nodeIdentifiers) + element.setScaleFactors(eft4, scalefactors) + + elementIdentifier = elementIdentifier + 1 + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + int(elementsCountAroundTC*0.5) + + N*elementsCountAroundTC + eTC] == annotationGroup._name or \ + annotationArrayAlong[0] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + for eTC in range(int(elementsCountAroundTC * 0.5)): + va = int(elementsCountAround - elementsCountAroundTC * 0.5) + eTC + vb = (va + 1) % elementsCountAround + # set general linear map coefficients + radiansAround = va * radiansPerElementAround + radiansAroundNext = vb * radiansPerElementAround + scalefactors = [ -1.0, + math.sin(radiansAround), math.cos(radiansAround), radiansPerElementAround, + math.sin(radiansAroundNext), math.cos(radiansAroundNext), radiansPerElementAround] + + bni21 = elementsCountThroughWall + 1 + bni22 = bni21 + elementsCountAround * (elementsCountThroughWall + 1) - \ + int(elementsCountAroundTC*0.5) + eTC + 1 + if elementsCountAroundTC > 2: + bni23 = bni21 + elementsCountAround * elementsCountThroughWall + 1 \ + if eTC == int(elementsCountAroundTC * 0.5 - 1) \ + else bni22 + 1 + bni31 = bni21 + elementsCountAround * (elementsCountThroughWall + 1) + \ + int(tcCount - 1) * (elementsCountAroundTC - 1) + \ + int(elementsCountAroundTC * 0.5) + eTC + (0 if eTC > 0 else 1) + else: + bni23 = bni21 + elementsCountAround * elementsCountThroughWall + eTC + 1 + bni31 = bni21 + elementsCountAround * (elementsCountThroughWall + 1) + 1 + bni32 = bni21 + elementsCountAround * (elementsCountThroughWall + 1) + 1 \ + if eTC == int(elementsCountAroundTC * 0.5 - 1) else bni31 + 1 + + tetrahedralElement = (eTC == 0) + if tetrahedralElement: + eft5 = eftfactory.createEftTetrahedronXi1Zero(va * 10000, vb * 10000) + elementtemplate5.defineField(coordinates, -1, eft5) + nodeIdentifiers = [bni21, bni22, bni23, bni31] + else: + eft6 = eftfactory.createEftPyramidBottomSimple(va * 10000, vb * 10000) + elementtemplate6.defineField(coordinates, -1, eft6) + nodeIdentifiers = [bni21, bni22, bni23, bni31, bni32] + + element = mesh.createElement(elementIdentifier, elementtemplate5 if tetrahedralElement else elementtemplate6) + element.setNodesByIdentifier(eft5 if tetrahedralElement else eft6, nodeIdentifiers) + element.setScaleFactors(eft5 if tetrahedralElement else eft6, scalefactors) + elementIdentifier = elementIdentifier + 1 + for annotationGroup in annotationGroups: + if annotationArrayAround[elementsCountAround + eTC] == annotationGroup._name or \ + annotationArrayAlong[0] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + # create regular elements + now = elementsCountAround * (elementsCountThroughWall + 1) + tcOffset = (elementsCountAroundTC - 1) * tcCount + for e2 in range(1 if closedProximalEnd else 0, elementsCountAlong): + tcOffset1 = (e2 - 1 if closedProximalEnd else e2) * (elementsCountAroundTC - 1) * tcCount + for e3 in range(elementsCountThroughWall): + for e1 in range(elementsCountAround): + if closedProximalEnd: + bni11 = (e2 - 1) * now + e3 * elementsCountAround + e1 + 1 + (elementsCountThroughWall + 1) + \ + tcOffset1 + bni12 = (e2 - 1) * now + e3 * elementsCountAround + (e1 + 1) % elementsCountAround + 1 + \ + (elementsCountThroughWall + 1) + tcOffset1 + bni21 = (e2 - 1) * now + (e3 + 1) * elementsCountAround + e1 + 1 + \ + (elementsCountThroughWall + 1) + tcOffset1 + bni22 = (e2 - 1) * now + (e3 + 1) * elementsCountAround + (e1 + 1) % elementsCountAround + 1 + \ + (elementsCountThroughWall + 1) + tcOffset1 + else: + bni11 = e2 * now + e3 * elementsCountAround + e1 + 1 + tcOffset1 + bni12 = e2 * now + e3 * elementsCountAround + (e1 + 1) % elementsCountAround + 1 + tcOffset1 + bni21 = e2 * now + (e3 + 1) * elementsCountAround + e1 + 1 + tcOffset1 + bni22 = e2 * now + (e3 + 1) * elementsCountAround + (e1 + 1) % elementsCountAround + 1 + tcOffset1 + nodeIdentifiers = [bni11, bni12, bni11 + now + tcOffset, bni12 + now + tcOffset, + bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset] onOpening = e1 > elementsCountAround - 2 element = mesh.createElement(elementIdentifier, elementtemplate) element.setNodesByIdentifier(eft, nodeIdentifiers) - element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - element.merge(textureElementtemplate2 if onOpening else textureElementtemplate1) - element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture3, nodeIdentifiers) + if xFlat and xTexture: + element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) + element.merge(textureElementtemplate2 if onOpening else textureElementtemplate1) + element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture3, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 if annotationGroups: for annotationGroup in annotationGroups: @@ -1781,11 +2000,20 @@ def createNodesAndElementsTeniaColi(region, # Add elements for tenia coli for eTC in range(int(elementsCountAroundTC*0.5)): - bni21 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 1 + tcOffset1 - bni22 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 2 + tcOffset1 - bni31 = (e2+1)*now + eTC + 1 + tcOffset1 - bni32 = (e2+1)*now + eTC + 2 + tcOffset1 - if eTC < int(elementsCountAroundTC*0.5) - 1: + if closedProximalEnd: + bni21 = elementsCountThroughWall + 1 + (e2 - 1) * now + elementsCountThroughWall * elementsCountAround + \ + eTC + 1 + tcOffset1 + bni22 = elementsCountThroughWall + 1 + (e2 - 1) * now + elementsCountThroughWall * elementsCountAround + \ + eTC + 2 + tcOffset1 + bni31 = elementsCountThroughWall + 1 + e2 * now + eTC + 1 + tcOffset1 + bni32 = elementsCountThroughWall + 1 + e2 * now + eTC + 2 + tcOffset1 + else: + bni21 = e2 * now + elementsCountThroughWall * elementsCountAround + eTC + 1 + tcOffset1 + bni22 = e2 * now + elementsCountThroughWall * elementsCountAround + eTC + 2 + tcOffset1 + bni31 = (e2 + 1) * now + eTC + 1 + tcOffset1 + bni32 = (e2 + 1) * now + eTC + 2 + tcOffset1 + + if eTC < int(elementsCountAroundTC * 0.5) - 1: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] else: @@ -1793,9 +2021,10 @@ 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) - 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) + if xFlat and xTexture: + 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) elementIdentifier = elementIdentifier + 1 for annotationGroup in annotationGroups: if annotationArrayAround[elementsCountAround + eTC] == annotationGroup._name or \ @@ -1805,38 +2034,57 @@ def createNodesAndElementsTeniaColi(region, for N in range(tcCount - 1): for eTC in range(elementsCountAroundTC): - bni21 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 1 + tcOffset1 + \ - int(elementsCountAroundTC*0.5) + (N+1)*elementsCountAroundHaustrum + N*elementsCountAroundTC - bni22 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 2 + tcOffset1 + \ - int(elementsCountAroundTC*0.5) + (N+1)*elementsCountAroundHaustrum + N*elementsCountAroundTC - bni31 = (e2+1)*now + eTC + 1 + tcOffset1 + int(elementsCountAroundTC*0.5) - 1 + \ - N*(elementsCountAroundTC-1) - bni32 = (e2+1)*now + eTC + 2 + tcOffset1 + int(elementsCountAroundTC*0.5) - 1 + \ - N*(elementsCountAroundTC-1) + if closedProximalEnd: + bni21 = elementsCountThroughWall + 1 + ( + e2 - 1) * now + elementsCountThroughWall * elementsCountAround \ + + eTC + 1 + tcOffset1 + int(elementsCountAroundTC * 0.5) + \ + (N + 1) * elementsCountAroundHaustrum + N * elementsCountAroundTC + bni22 = elementsCountThroughWall + 1 + ( + e2 - 1) * now + elementsCountThroughWall * elementsCountAround \ + + eTC + 2 + tcOffset1 + int(elementsCountAroundTC * 0.5) + \ + (N + 1) * elementsCountAroundHaustrum + N * elementsCountAroundTC + bni31 = elementsCountThroughWall + 1 + e2 * now + eTC + 1 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) - 1 + N * (elementsCountAroundTC - 1) + bni32 = elementsCountThroughWall + 1 + e2 * now + eTC + 2 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) - 1 + N * (elementsCountAroundTC - 1) + else: + bni21 = e2 * now + elementsCountThroughWall * elementsCountAround + eTC + 1 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) + ( + N + 1) * elementsCountAroundHaustrum + N * elementsCountAroundTC + bni22 = e2 * now + elementsCountThroughWall * elementsCountAround + eTC + 2 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) + ( + N + 1) * elementsCountAroundHaustrum + N * elementsCountAroundTC + bni31 = (e2 + 1) * now + eTC + 1 + tcOffset1 + int(elementsCountAroundTC * 0.5) - 1 + \ + N * (elementsCountAroundTC - 1) + bni32 = (e2 + 1) * now + eTC + 2 + tcOffset1 + int(elementsCountAroundTC * 0.5) - 1 + \ + N * (elementsCountAroundTC - 1) if eTC == 0: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni32, bni32 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate2) element.setNodesByIdentifier(eft2, nodeIdentifiers) - element.merge(flatElementtemplate4) - element.merge(textureElementtemplate4) - element.setNodesByIdentifier(eftTexture6, nodeIdentifiers) + if xFlat and xTexture: + element.merge(flatElementtemplate4) + element.merge(textureElementtemplate4) + element.setNodesByIdentifier(eftTexture6, 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) - element.merge(flatElementtemplate1) - element.merge(textureElementtemplate1) - element.setNodesByIdentifier(eftTexture3, nodeIdentifiers) + if xFlat and xTexture: + element.merge(flatElementtemplate1) + element.merge(textureElementtemplate1) + element.setNodesByIdentifier(eftTexture3, nodeIdentifiers) else: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni31 + now + tcOffset] element = mesh.createElement(elementIdentifier, elementtemplate1) element.setNodesByIdentifier(eft1, nodeIdentifiers) - element.merge(flatElementtemplate3) - element.merge(textureElementtemplate3) - element.setNodesByIdentifier(eftTexture5, nodeIdentifiers) + if xFlat and xTexture: + element.merge(flatElementtemplate3) + element.merge(textureElementtemplate3) + element.setNodesByIdentifier(eftTexture5, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 for annotationGroup in annotationGroups: if annotationArrayAround[elementsCountAround + int(elementsCountAroundTC*0.5) + @@ -1846,13 +2094,26 @@ def createNodesAndElementsTeniaColi(region, meshGroup.addElement(element) for eTC in range(int(elementsCountAroundTC*0.5)): - bni21 = e2*now + (elementsCountThroughWall)*elementsCountAround + eTC + 1 + tcOffset1 + \ - int(elementsCountAroundTC*0.5) + tcCount*elementsCountAroundHaustrum + \ - (tcCount - 1)*elementsCountAroundTC - bni22 = e2*now + (elementsCountThroughWall)*elementsCountAround + 1 + tcOffset1 if eTC == int(elementsCountAroundTC*0.5 - 1) else bni21 + 1 - bni31 = (e2+1)*now + eTC + 1 + tcOffset1 + int(elementsCountAroundTC*0.5) - 1 + \ - (tcCount-1)*(elementsCountAroundTC-1) - bni32 = (e2+1)*now + 1 + tcOffset1 if eTC == int(elementsCountAroundTC*0.5 - 1) else bni31 + 1 + if closedProximalEnd: + bni21 = elementsCountThroughWall + 1 + (e2 - 1) * now + elementsCountThroughWall * elementsCountAround + \ + eTC + 1 + tcOffset1 + int(elementsCountAroundTC * 0.5) + tcCount * elementsCountAroundHaustrum + \ + (tcCount - 1) * elementsCountAroundTC + bni22 = elementsCountThroughWall + 1 + (e2 - 1) * now + elementsCountThroughWall * elementsCountAround + \ + 1 + tcOffset1 if eTC == int(elementsCountAroundTC * 0.5 - 1) else bni21 + 1 + bni31 = elementsCountThroughWall + 1 + e2 * now + eTC + 1 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) - 1 + (tcCount - 1) * (elementsCountAroundTC - 1) + bni32 = elementsCountThroughWall + 1 + e2 * now + 1 + tcOffset1 if eTC == int( + elementsCountAroundTC * 0.5 - 1) \ + else bni31 + 1 + else: + bni21 = e2 * now + (elementsCountThroughWall) * elementsCountAround + eTC + 1 + tcOffset1 + \ + int(elementsCountAroundTC * 0.5) + tcCount * elementsCountAroundHaustrum + \ + (tcCount - 1) * elementsCountAroundTC + bni22 = e2 * now + (elementsCountThroughWall) * elementsCountAround + 1 + tcOffset1 if eTC == int( + elementsCountAroundTC * 0.5 - 1) else bni21 + 1 + bni31 = (e2 + 1) * now + eTC + 1 + tcOffset1 + int(elementsCountAroundTC * 0.5) - 1 + \ + (tcCount - 1) * (elementsCountAroundTC - 1) + bni32 = (e2 + 1) * now + 1 + tcOffset1 if eTC == int(elementsCountAroundTC * 0.5 - 1) else bni31 + 1 if eTC > 0: nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] @@ -1862,14 +2123,15 @@ 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 eTC > 0: - element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - element.merge(textureElementtemplate2 if onOpening else textureElementtemplate1) - element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture3, nodeIdentifiers) - else: - element.merge(flatElementtemplate5 if onOpening else flatElementtemplate4) - element.merge(textureElementtemplate5 if onOpening else textureElementtemplate4) - element.setNodesByIdentifier(eftTexture7 if onOpening else eftTexture6, nodeIdentifiers) + if xFlat and xTexture: + 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) + else: + element.merge(flatElementtemplate5 if onOpening else flatElementtemplate4) + element.merge(textureElementtemplate5 if onOpening else textureElementtemplate4) + element.setNodesByIdentifier(eftTexture7 if onOpening else eftTexture6, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 for annotationGroup in annotationGroups: if annotationArrayAround[elementsCountAround + int(elementsCountAroundTC*(tcCount - 0.5)) + eTC] \ diff --git a/src/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py b/src/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py index a9307803..ed7c5146 100644 --- a/src/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py +++ b/src/scaffoldmaker/utils/eftfactory_bicubichermitelinear.py @@ -268,3 +268,166 @@ def createEftWedgeXi1ZeroOpenTube(self): remapEftLocalNodes(eft, 6, ln_map) assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1ZeroOpenTube: Failed to validate eft' return eft + + def createEftTetrahedronXi1One(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a bicubic hermite linear element field for a solid tetrahedron for the apex of cecum, + with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 1 and xi2 = 1. + Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full bicubic hermite linear + eft = self._mesh.createElementfieldtemplate(self._basis) + + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 4 + 4, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ] ) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), + (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])]) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), + (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + + ln_map = [ 1, 1, 2, 3, 1, 1, 4, 3] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftTetrahedronXi1One: Failed to validate eft' + return eft + + def createEftTetrahedronXi1Zero(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a bicubic hermite linear element field for a solid tetrahedron for the apex of cecum, + with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 0, xi2 = 1. + Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full bicubic hermite linear + eft = self._mesh.createElementfieldtemplate(self._basis) + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 4 + 4, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), + (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])]) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), + (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + + ln_map = [ 1, 1, 2, 3, 1, 1, 2, 4] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftTetrahedronXi1Zero: Failed to validate eft' + return eft + + def createEftPyramidBottomSimple(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a bicubic hermite linear element field for a solid pyramid for elements within + a tenia coli joining to the cecal apex, with xi1 and xi3 collapsed on xi2 = 0. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of the + radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Need to create a new template for each + sector around axis giving common nodeScaleFactorOffset values on common faces. + Suggestion is to start at 0 and add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full bicubic hermite linear + eft = self._mesh.createElementfieldtemplate(self._basis) + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 4 + 4, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3]) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])]) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + + ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ] + remapEftLocalNodes(eft, 5, ln_map) + + assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftPyramidBottomSimple: Failed to validate eft' + return eft diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 0934420a..3f141ccd 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -355,7 +355,7 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs ln_map = [ 1, 1, 2, 2, 1, 1, 3, 4 ] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronTop: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' return eft def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): @@ -437,6 +437,157 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronTop: Failed to validate eft' return eft + def createEftTetrahedronXi1One(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for a solid tetrahedron for the apex of cecum, + with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 1 and xi2 = 1. + Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 8 + 4, 0) + eft.setFunctionNumberOfTerms(n * 8 + 6, 0) + eft.setFunctionNumberOfTerms(n * 8 + 7, 0) + eft.setFunctionNumberOfTerms(n * 8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ] ) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), + (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), + (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2), [1]]) + + # remap parameters on xi1 = 1, xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [4, 8], Node.VALUE_LABEL_D_DS3, []) + remapEftNodeValueLabel(eft, [4, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [4, 8], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [4, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 1, 2, 3, 1, 1, 4, 3] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronXi1One: Failed to validate eft' + return eft + + def createEftTetrahedronXi1Zero(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for a solid tetrahedron for the apex of cecum, + with xi1 and xi3 collapsed on xi2 = 0, and xi3 collapsed on xi1 = 0, xi2 = 1. + Each collapsed node on xi2 = 0 has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 8 + 4, 0) + eft.setFunctionNumberOfTerms(n * 8 + 6, 0) + eft.setFunctionNumberOfTerms(n * 8 + 7, 0) + eft.setFunctionNumberOfTerms(n * 8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), + (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3])]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), + (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ln], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2), [1]]) + + # remap parameters on xi1 = 0, xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS3, []) + remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 1, 2, 3, 1, 1, 2, 4] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronXi1Zero: Failed to validate eft' + return eft + def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' @@ -499,7 +650,7 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ] remapEftLocalNodes(eft, 5, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidBottom: Failed to validate eft' return eft def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): @@ -566,6 +717,70 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' return eft + def createEftPyramidBottomSimple(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for a solid pyramid for elements within + a tenia coli joining to the cecal apex, with xi1 and xi3 collapsed on xi2 = 0. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of the + radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Need to create a new template for each + sector around axis giving common nodeScaleFactorOffset values on common faces. + Suggestion is to start at 0 and add 10000 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n * 8 + 4, 0) + eft.setFunctionNumberOfTerms(n * 8 + 6, 0) + eft.setFunctionNumberOfTerms(n * 8 + 7, 0) + eft.setFunctionNumberOfTerms(n * 8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3]) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) + + for layer in range(2): + soAround = 1 + ln = layer * 4 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 1]), (Node.VALUE_LABEL_D_DS2, [soAround + 2])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer * 4 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 4]), (Node.VALUE_LABEL_D_DS2, [soAround + 5])]) + # 2 terms for cross derivative 1 2 to correct circular apex: cos(theta).phi, -sin(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, + [(Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6]), + (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6])]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ] + remapEftLocalNodes(eft, 5, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidBottomSimple: Failed to validate eft' + return eft + def createEftWedgeXi1One(self): ''' Create a basic tricubic hermite element template for elements From 811e9cc3495a3290885d363389d5269629d93495 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Jul 2020 12:14:41 +1200 Subject: [PATCH 11/12] Update unittest for cecum with addition of tenia coli --- tests/test_cecum.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/test_cecum.py b/tests/test_cecum.py index 9aa92374..3083a00a 100644 --- a/tests/test_cecum.py +++ b/tests/test_cecum.py @@ -48,25 +48,25 @@ def test_cecum1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_cecum1.generateBaseMesh(region, options) - self.assertEqual(0, len(annotationGroups)) + self.assertEqual(1, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(1252, mesh3d.getSize()) + self.assertEqual(1492, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(5017, mesh2d.getSize()) + self.assertEqual(5617, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(6287, mesh1d.getSize()) + self.assertEqual(6767, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(2522, nodes.getSize()) + self.assertEqual(2642, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) - assertAlmostEqualList(self, minimums, [-49.01602978165163, -46.896253729833276, -2.0], 1.0E-6) + assertAlmostEqualList(self, minimums, [-49.01602978165163, -46.896253729833276, -2.343705425292791], 1.0E-6) assertAlmostEqualList(self, maximums, [42.18397974244375, 54.906723342953256, 180.0], 1.0E-6) with ChangeManager(fieldmodule): @@ -79,10 +79,10 @@ def test_cecum1(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 65800.2656669658, delta=1.0E-6) + self.assertAlmostEqual(surfaceArea, 65958.89379225131, delta=1.0E-6) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 127049.69543528013, delta=1.0E-6) + self.assertAlmostEqual(volume, 127903.31774270078, delta=1.0E-6) if __name__ == "__main__": unittest.main() From f0502851fab8c81567d65f447f1010d664f1022e Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Jul 2020 16:35:11 +1200 Subject: [PATCH 12/12] Allow for annotation of elements through wall --- src/scaffoldmaker/annotation/colon_terms.py | 4 ++ .../meshtypes/meshtype_3d_cecum1.py | 36 ++++++++---- .../meshtypes/meshtype_3d_colon1.py | 12 ++-- .../meshtypes/meshtype_3d_colonsegment1.py | 58 ++++++++++++------- .../meshtypes/meshtype_3d_smallintestine1.py | 6 +- src/scaffoldmaker/utils/tubemesh.py | 13 ++++- tests/test_cecum.py | 2 +- 7 files changed, 89 insertions(+), 42 deletions(-) diff --git a/src/scaffoldmaker/annotation/colon_terms.py b/src/scaffoldmaker/annotation/colon_terms.py index 9f172a35..8250e7de 100644 --- a/src/scaffoldmaker/annotation/colon_terms.py +++ b/src/scaffoldmaker/annotation/colon_terms.py @@ -6,12 +6,16 @@ colon_terms = [ ( "ascending colon", "UBERON:0001156", "FMA:14545", "ILX:0734393"), ( "descending colon", "UBERON:0001158", "FMA:14547", "ILX:0724444"), + ( "caecum", "UBERON:0001153", "FMA:14541", "ILX:0732270"), ( "colon", "UBERON:0001155", "FMA:14543", "ILX:0736005"), + ( "colonic mucosa", "UBERON:0000317", "FMA:14984", "ILX:0731046"), ( "distal colon", "UBERON:0008971", "ILX:0727523"), ( "mesenteric zone", None), ( "non-mesenteric zone", None), ( "proximal colon", "UBERON:0008972", "ILX:0733240"), + ( "serosa of colon", "UBERON:0003335", "FMA:14990", "ILX:0736932"), ( "spiral colon", "UBERON:0010239", "ILX:0735018"), + ( "submucosa of colon", "UBERON:0003331", "FMA:14985", "ILX:0728041"), ( "taenia coli", "UBERON:0012419", "FMA:15041", "ILX:0731555"), ( "taenia libera", "ILX:0739285"), ( "taenia mesocolica", "ILX:0739284"), diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 24c632a1..bf7dddb7 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -6,6 +6,8 @@ import copy import math +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +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 ColonSegmentTubeMeshInnerPoints, \ getFullProfileFromHalfHaustrum, getTeniaColi, createNodesAndElementsTeniaColi @@ -364,11 +366,16 @@ def generateBaseMesh(cls, region, options): segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongCecum, innerRadiusAlongCecum, dInnerRadiusAlongCecum, tcWidthAlongCecum, startPhase) - annotationArrayAlong = [''] * elementsCountAlongSegment * segmentCount + # Create annotation + cecumGroup = AnnotationGroup(region, get_colon_term("caecum")) + annotationGroups = [cecumGroup] + annotationArrayAlong = (['caecum'] * elementsCountAlong) + + annotationArrayThroughWall = [''] * elementsCountThroughWall for nSegment in range(segmentCount): # Make regular segments - xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround \ + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroupsAround, annotationArrayAround \ = colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) # Replace first half of first segment with apex and sample along apex and second half of segment @@ -416,6 +423,7 @@ def generateBaseMesh(cls, region, options): zRefList, innerRadiusAlongCecum, closedProximalEnd) # Create coordinates and derivatives + annotationGroups += annotationGroupsAround wallThicknessList = [wallThickness] * (elementsCountAlong + 1) xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, @@ -466,15 +474,17 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) else: nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xCecum, d1Cecum, d2Cecum, d3Cecum, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) # Add ostium on track surface between two tenia on the last segment elementsAroundTrackSurface = elementsCountAroundHaustrum @@ -665,14 +675,16 @@ def generateBaseMesh(cls, region, options): ostiumSettings['Number of elements around ostium'] = len(innerEndPoints_Id) - nextNodeIdentifier, nextElementIdentifier, (o1_x, o1_d1, o1_d2, o1_d3, o1_NodeId, o1_Positions) = \ - generateOstiumMesh(region, ostiumSettings, trackSurfaceOstium, centrePosition, axis1, - nextNodeIdentifier, nextElementIdentifier) - fm = region.getFieldmodule() mesh = fm.findMeshByDimension(3) nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + cecumMeshGroup = cecumGroup.getMeshGroup(mesh) + + nextNodeIdentifier, nextElementIdentifier, (o1_x, o1_d1, o1_d2, o1_d3, o1_NodeId, o1_Positions) = \ + generateOstiumMesh(region, ostiumSettings, trackSurfaceOstium, centrePosition, axis1, + nextNodeIdentifier, nextElementIdentifier, ostiumMeshGroups= [cecumMeshGroup] ) + startProportions = [] for n in range(len(innerEndPoints_Id)): startProportions.append(trackSurfaceOstium.getProportion(o1_Positions[n])) @@ -681,8 +693,8 @@ def generateBaseMesh(cls, region, options): nodes, mesh, nextNodeIdentifier, nextElementIdentifier, o1_x, o1_d1, o1_d2, None, o1_NodeId, None, endPoints_x, endPoints_d1, endPoints_d2, None, endPoints_Id, endDerivativesMap, - elementsCountRadial = 2, tracksurface = trackSurfaceOstium, startProportions = startProportions, - endProportions = endProportions) + elementsCountRadial = 2, meshGroups = [cecumMeshGroup], tracksurface = trackSurfaceOstium, + startProportions = startProportions, endProportions = endProportions) # Delete elements under annulus mesh deleteElementsAndNodesUnderAnnulusMesh(fm, nodes, mesh, deleteElementIdentifier, deleteNodeIdentifier) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index f1f72762..d42d1771 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -465,6 +465,8 @@ def generateBaseMesh(cls, region, options): ['transverse colon'] * elementsAlongInTransverse + ['descending colon'] * elementsAlongInDistal) + annotationArrayThroughWall = ([''] * elementsCountThroughWall) + xExtrude = [] d1Extrude = [] d2Extrude = [] @@ -536,8 +538,9 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) else: # Create flat and texture coordinates @@ -549,8 +552,9 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 11865541..e66f579a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -260,14 +260,26 @@ def generateBaseMesh(cls, region, options): segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactorAlongSegment, radiusAlongSegment, dRadiusAlongSegment, tcWidthAlongSegment, startPhase) + # Create annotation + annotationGroups = [] annotationArrayAlong = [''] * elementsCountAlongSegment + # 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 = [serosaGroup, submucosaGroup, mucosaGroup] + # annotationGroups += annotationGroupsThroughWall + # annotationArrayThroughWall = ['colonic mucosa', 'submucosa of colon', 'serosa of colon'] + annotationArrayThroughWall = [''] * elementsCountThroughWall # Create inner points nSegment = 0 + closedProximalEnd = False - xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArrayAround = \ + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroupsAround, annotationArrayAround = \ colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) + annotationGroups += annotationGroupsAround + # Project reference point for warping onto central path sxRefList, sd1RefList, sd2ProjectedListRef, zRefList = \ tubemesh.getPlaneProjectionOnCentralPath(xInner, elementsCountAround, elementsCountAlongSegment, @@ -277,7 +289,7 @@ def generateBaseMesh(cls, region, options): xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = tubemesh.warpSegmentPoints( xInner, d1Inner, d2Inner, segmentAxis, sxRefList, sd1RefList, sd2ProjectedListRef, elementsCountAround, elementsCountAlongSegment, zRefList, radiusAlongSegment, - closedProximalEnd=False) + closedProximalEnd) contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() @@ -288,8 +300,6 @@ def generateBaseMesh(cls, region, options): relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() - closedProximalEnd = False - if tcThickness > 0: tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() xList, d1List, d2List, d3List, annotationGroups, annotationArrayAround = getTeniaColi( @@ -307,8 +317,9 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = createNodesAndElementsTeniaColi( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, - tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + tcCount, annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) else: # Create flat and texture coordinates xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( @@ -319,8 +330,9 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, - firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd) + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd) return annotationGroups @@ -384,9 +396,9 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment): nSegment * self._elementsCountAlongSegment: (nSegment + 1) * self._elementsCountAlongSegment + 1] - xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, \ - contractedWallThicknessSegment, segmentAxis, annotationGroups, annotationArray\ - = getColonSegmentInnerPoints(self._region, + xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, contractedWallThicknessSegment, \ + segmentAxis, annotationGroupsAround, annotationArrayAround = \ + getColonSegmentInnerPoints(self._region, self._elementsCountAroundTC, self._elementsCountAroundHaustrum, self._elementsCountAlongSegment, self._tcCount, self._segmentLengthEndDerivativeFactor, self._segmentLengthMidDerivativeFactor, self._segmentLength, self._wallThickness, @@ -407,7 +419,8 @@ def getColonSegmentTubeMeshInnerPoints(self, nSegment): contractedWallThickness = contractedWallThicknessSegment[startIdx:self._elementsCountAlongSegment + 1] self._contractedWallThicknessList += contractedWallThickness - return xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray + return xInner, d1Inner, d2Inner, transitElementList, segmentAxis, \ + annotationGroupsAround, annotationArrayAround def getTubeTCWidthList(self): return self._tubeTCWidthList @@ -471,7 +484,7 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, along colon segment. Assume incompressiblity and a shortened length around will result in a thicker wall and vice-versa. :return segmentAxis: Axis of segment. - :return annotationGroups, annotationArray: annotationArray stores annotation + :return annotationGroupsAround, annotationArrayAround: annotationArray stores annotation names of elements around. """ @@ -545,8 +558,8 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, # Create annotation groups for mouse colon mzGroup = AnnotationGroup(region, get_colon_term("mesenteric zone")) nonmzGroup = AnnotationGroup(region, get_colon_term("non-mesenteric zone")) - annotationGroups = [mzGroup, nonmzGroup] - annotationArray = (['mesenteric zone']*int(elementsCountAroundTC*0.5) + + annotationGroupsAround = [mzGroup, nonmzGroup] + annotationArrayAround = (['mesenteric zone']*int(elementsCountAroundTC*0.5) + ['non-mesenteric zone']*elementsCountAroundHaustrum + ['mesenteric zone']*int(elementsCountAroundTC*0.5)) @@ -746,11 +759,11 @@ def getColonSegmentInnerPoints(region, elementsCountAroundTC, d2Final = d2Final + d2AlongList # Create annotation groups - annotationGroups = [] - annotationArray = ['']*(elementsCountAround) + annotationGroupsAround = [] + annotationArrayAround = ['']*(elementsCountAround) return xFinal, d1Final, d2Final, transitElementList, xiList, relaxedLengthList, contractedWallThicknessList, \ - segmentAxis, annotationGroups, annotationArray + segmentAxis, annotationGroupsAround, annotationArrayAround def createHalfSetInterHaustralSegment(elementsCountAroundTC, elementsCountAroundHaustrum, tcCount, tcWidth, radius, cornerInnerRadiusFactor, sampleElementOut): @@ -1540,7 +1553,7 @@ def createNodesAndElementsTeniaColi(region, xTexture, d1Texture, d2Texture, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, tcCount, - annotationGroups, annotationArrayAround, annotationArrayAlong, + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ @@ -1560,6 +1573,7 @@ def createNodesAndElementsTeniaColi(region, :param annotationGroups: stores information about annotation groups. :param annotationArrayAround: stores annotation names of elements around. :param annotationArrayAlong: stores annotation names of elements along. + :param annotationArrayThroughWall: stores annotation names of elements through wall. :param firstNodeIdentifier, firstElementIdentifier: first node and element identifier to use. :param useCubicHermiteThroughWall: use linear when false. @@ -1822,7 +1836,8 @@ def createNodesAndElementsTeniaColi(region, elementIdentifier = elementIdentifier + 1 for annotationGroup in annotationGroups: if annotationArrayAround[e1] == annotationGroup._name or \ - annotationArrayAlong[0] == annotationGroup._name: + annotationArrayAlong[0] == annotationGroup._name or \ + annotationArrayThroughWall[e3] == annotationGroup._name: meshGroup = annotationGroup.getMeshGroup(mesh) meshGroup.addElement(element) @@ -1994,7 +2009,8 @@ def createNodesAndElementsTeniaColi(region, if annotationGroups: for annotationGroup in annotationGroups: if annotationArrayAround[e1] == annotationGroup._name or \ - annotationArrayAlong[e2] == annotationGroup._name: + annotationArrayAlong[e2] == annotationGroup._name or \ + annotationArrayThroughWall[e3] == annotationGroup._name: meshGroup = annotationGroup.getMeshGroup(mesh) meshGroup.addElement(element) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 1a969158..b0b72817 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -272,6 +272,7 @@ def generateBaseMesh(cls, region, options): ['jejunum'] * elementsAlongJejunum + ['ileum'] * elementsAlongIleum) annotationArrayAround = [''] * (elementsCountAround) + annotationArrayThroughWall = [''] * elementsCountThroughWall xExtrude = [] d1Extrude = [] @@ -348,8 +349,9 @@ def generateBaseMesh(cls, region, options): nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, firstNodeIdentifier, firstElementIdentifier, - useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd=False) + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, + firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, + closedProximalEnd=False) return annotationGroups diff --git a/src/scaffoldmaker/utils/tubemesh.py b/src/scaffoldmaker/utils/tubemesh.py index 5d7342de..23d14393 100644 --- a/src/scaffoldmaker/utils/tubemesh.py +++ b/src/scaffoldmaker/utils/tubemesh.py @@ -438,7 +438,7 @@ def createNodesAndElements(region, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, elementsCountAround, elementsCountAlong, elementsCountThroughWall, - annotationGroups, annotationArrayAround, annotationArrayAlong, + annotationGroups, annotationArrayAround, annotationArrayAlong, annotationArrayThroughWall, firstNodeIdentifier, firstElementIdentifier, useCubicHermiteThroughWall, useCrossDerivatives, closedProximalEnd): """ @@ -455,6 +455,7 @@ def createNodesAndElements(region, :param annotationGroups: stores information about annotation groups. :param annotationArrayAround: stores annotation names of elements around. :param annotationArrayAlong: stores annotation names of elements along. + :param annotationArrayThroughWall: stores annotation names of elements through wall. :param firstNodeIdentifier, firstElementIdentifier: first node and element identifier to use. :param useCubicHermiteThroughWall: use linear when false @@ -639,6 +640,13 @@ def createNodesAndElements(region, ] result = element.setScaleFactors(eft1, scalefactors) elementIdentifier = elementIdentifier + 1 + if annotationGroups: + for annotationGroup in annotationGroups: + if annotationArrayAround[e1] == annotationGroup._name or \ + annotationArrayAlong[0] == annotationGroup._name or\ + annotationArrayThroughWall[e3] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) # Create regular elements now = elementsCountAround * (elementsCountThroughWall + 1) @@ -669,7 +677,8 @@ def createNodesAndElements(region, if annotationGroups: for annotationGroup in annotationGroups: if annotationArrayAround[e1] == annotationGroup._name or \ - annotationArrayAlong[e2] == annotationGroup._name: + annotationArrayAlong[e2] == annotationGroup._name or\ + annotationArrayThroughWall[e3] == annotationGroup._name: meshGroup = annotationGroup.getMeshGroup(mesh) meshGroup.addElement(element) diff --git a/tests/test_cecum.py b/tests/test_cecum.py index 3083a00a..afa48e1c 100644 --- a/tests/test_cecum.py +++ b/tests/test_cecum.py @@ -48,7 +48,7 @@ def test_cecum1(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_cecum1.generateBaseMesh(region, options) - self.assertEqual(1, len(annotationGroups)) + self.assertEqual(2, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())