diff --git a/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 13e32185..083d0f67 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -6,13 +6,13 @@ import copy from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion -from scaffoldmaker.meshtypes.meshtype_3d_colonsegmentsimplemesentery1 import MeshType_3d_colonsegmentsimplemesentery1, getColonSegmentInnerPointsNoTeniaColi -from scaffoldmaker.meshtypes.meshtype_3d_colonsegmentteniacoli1 import MeshType_3d_colonsegmentteniacoli1, getColonSegmentInnerPointsTeniaColi, getTeniaColi +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 import vector from scaffoldmaker.utils import zinc_utils from opencmiss.zinc.node import Node @@ -34,15 +34,15 @@ class MeshType_3d_colon1(Scaffold_base): }, 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ 0.0, 0.0, 0.0 ], [ -1.6, 5.8, 0.0 ], [ -2.4, -0.6, -1.2 ], [ -1.4, -0.1, -1.2 ] ], - [ [ -1.6, 6.1, 0.0 ], [ -0.8, 6.2, 0.0 ], [ -2.2, -0.4, -0.8 ], [ -0.4, 1.9, 2.2 ] ], - [ [ -0.3, 12.3, 0.0 ], [ 4.8, 1.2, 0.0 ], [ -1.0, 2.0, 0.8 ], [ -0.6, 0.0, 5.1 ] ], - [ [ 3.0, 12.2, 0.0 ], [ 4.0, -0.3, 0.0 ], [ -0.5, 0.4, 2.9 ], [ 0.0, 0.1, 2.4 ] ], - [ [ 7.3, 12.3, 0.0 ], [ 4.6, 1.1, 0.0 ], [ -0.2, 1.0, 2.2 ], [ 0.5, 2.5, -2.0 ] ], - [ [ 12.2, 12.9, 0.0 ], [ 5.1, -3.6, 0.0 ], [ 1.0, 1.7, 0.6 ], [ 0.1, -0.6, -3.5 ] ], - [ [ 13.9, 6.4, 0.0 ], [ 0.2, -4.0, 0.0 ], [ 2.0, 0.0, -2.0 ], [ 1.5, -0.1, -1.0 ] ], - [ [ 12.9, -0.3, 0.0 ], [ -4.1, -3.0, 0.0 ], [ 0.6, -0.9, -1.4 ], [ 0.8, -1.1, -1.3 ] ], - [ [ 7.7, 0.3, 0.0 ], [ -3.0, 1.7, 0.0 ], [ 0.1, -1.1, -1.8 ], [ 0.4, -1.2, -1.2 ] ] ] ) + [ [ 0.0, 0.0, 0.0 ], [ -50.7, 178.2, 0.0 ], [ -24.0, -6.0, -12.0 ], [ -14.0, -1.0, -12.0 ] ], + [ [ -47.4, 188.6, 0.0 ], [ -19.3, 177.1, 0.0 ], [ -22.0, -4.0, -8.0 ], [ -4.0, 19.0, 22.0 ] ], + [ [ -4.4, 396.5, 0.0 ], [ 206.0, 40.1, 0.0 ], [ -10.0, 20.0, 8.0 ], [ -6.0, 0.0, 51.0 ] ], + [ [ 130.0, 384.1, 0.0 ], [ 130.8, -40.5, 0.0 ], [ -5.0, 4.0, 29.0 ], [ 0.0, 1.0, 24.0 ] ], + [ [ 279.4, 383.0, 0.0 ], [ 118.0, 48.7, 0.0 ], [ -2.0, 10.0, 22.0 ], [ 5.0, 25.0, -20.0 ] ], + [ [ 443.9, 390.8, 0.0 ], [ 111.3, -97.0, 0.0 ], [ 10.0, 17.0, 6.0 ], [ 1.0, -6.0, -35.0 ] ], + [ [ 475.2, 168.0, 0.0 ], [ -0.8, -112.4, 0.0 ], [ 20.0, 0.0, -20.0 ], [ 15.0, -1.0, -10.0 ] ], + [ [ 432.6, -32.3, 0.0 ], [ -90.5, -59.0, 0.0 ], [ 6.0, -9.0, -14.0 ], [ 8.0, -11.0, -13.0 ] ], + [ [ 272.4, 7.5, 0.0 ], [ -79.0, 47.4, 0.0 ], [ 1.0, -11.0, -18.0 ], [ 4.0, -12.0, -12.0 ] ] ] ) } ), 'Human 2' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { @@ -52,86 +52,92 @@ class MeshType_3d_colon1(Scaffold_base): }, 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ 0.0, 0.0, 0.0 ], [ -1.6, 5.8, 0.0 ], [ -2.4, -0.6, -1.2 ], [ -1.4, -0.1, -1.2 ] ], - [ [ -1.6, 6.1, 0.0 ], [ -0.8, 6.2, 0.0 ], [ -2.2, -0.4, -0.8 ], [ -0.4, 1.9, 2.2 ] ], - [ [ -0.3, 12.3, 0.0 ], [ 3.8, -0.4, 2.6 ], [ -1.0, 2.0, 0.8 ], [ -0.6, 0.0, 5.1 ] ], - [ [ 3.6, 10.4, 2.3 ], [ 4.4, 3.9, 1.4 ], [ -0.5, 0.4, 2.9 ], [ 0.0, 0.1, 2.4 ] ], - [ [ 9.5, 11.7, 2.4 ], [ 5.9, -0.6, 2.0 ], [ -0.2, 1.0, 2.2 ], [ 0.5, 2.5, -2.0 ] ], - [ [ 12.9, 17.3, -1.1 ], [ 1.8, -1.5, -5.5 ], [ 1.0, 1.7, 0.6 ], [ 0.1, -0.6, -3.5 ] ], - [ [ 14.6, 7.8, -1.5 ], [ -0.2, -3.8, -0.1 ], [ 2.0, 0.0, -2.0 ], [ 1.5, -0.1, -1.0 ] ], - [ [ 12.9, -0.3, 0.0 ], [ -4.1, -3.0, 0.0 ], [ 0.6, -0.9, -1.4 ], [ 0.8, -1.1, -1.3 ] ], - [ [ 7.7, 0.3, 0.0 ], [ -3.0, 1.7, 0.0 ], [ 0.1, -1.1, -1.8 ], [ 0.4, -1.2, -1.2 ] ] ] ) + [ [ 0.0, 0.0, 0.0 ], [ -34.7, 104.1, -18.1 ], [ -24.0, -6.0, -12.0 ], [ -14.0, -1.0, -12.0 ] ], + [ [ -34.5, 114.0, -18.1 ], [ 1.2, 86.6, -3.4 ], [ -22.0, -4.0, -8.0 ], [ -4.0, 19.0, 22.0 ] ], + [ [ -19.1, 218.5, 5.5 ], [ 78.7, -7.1, 94.5 ], [ -10.0, 20.0, 8.0 ], [ -6.0, 0.0, 51.0 ] ], + [ [ 82.5, 189.1, 94.2 ], [ 84.5, 7.1, 71.6 ], [ -5.0, 4.0, 29.0 ], [ 0.0, 1.0, 24.0 ] ], + [ [ 226.6, 218.7, 85.7 ], [ 95.0, 91.3, -58.5 ], [ -2.0, 10.0, 22.0 ], [ 5.0, 25.0, -20.0 ] ], + [ [ 325.5, 381.7, -57.9 ], [ 229.2, -66.7, -20.4 ], [ 10.0, 17.0, 6.0 ], [ 1.0, -6.0, -35.0 ] ], + [ [ 354.0, 105.3, -24.4 ], [ -6.3, -143.7, 20.3 ], [ 20.0, 0.0, -20.0 ], [ 15.0, -1.0, -10.0 ] ], + [ [ 296.5, -121.2, -0.6 ], [ -90.5, -59.0, 0.0 ], [ 6.0, -9.0, -14.0 ], [ 8.0, -11.0, -13.0 ] ], + [ [ 169.8, -73.4, -33.5 ], [ -72.2, 43.4, -27.4 ], [ 1.0, -11.0, -18.0 ], [ 4.0, -12.0, -12.0 ] ] ] ) } ), 'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, 'Length' : 1.0, - 'Number of elements' : 8 + 'Number of elements' : 7 + }, + 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( + [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ + [ [ 0.0, 0.0, 0.0 ], [ 6.0, 12.0, -2.0 ], [ 2.0, 1.0, 2.0 ], [ 6.0, 0.0, 3.0 ] ], + [ [ -2.0, 11.0, -3.0 ], [ -8.0, 4.0, 9.0 ], [ 2.0, 2.0, 1.0 ], [ 0.0, 1.0, 2.0 ] ], + [ [ -3.0, 2.0, 3.0 ], [ -4.0, -8.0, 0.0 ], [ 2.0, -1.0, 2.0 ], [ 1.0, 0.0, 2.0 ] ], + [ [ -11.0, -3.0, -4.0 ], [ -8.0, -3.0, -7.0 ], [ 1.0, -2.0, 1.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -16.0, -4.0, 0.0 ], [ 4.0, -3.0, 14.0 ], [ 1.0, -3.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -7.0, -8.0, 0.0 ], [ 5.0, -1.0, -14.0 ], [ 0.0, -3.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -1.0, -6.0, -1.0 ], [ 2.0, -2.0, 9.0 ], [ 1.0, -3.0, -1.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -2.0, -14.0, 5.0 ], [ -2.0, -4.0, 2.0 ], [ 1.0, -2.0, -2.0 ], [ 0.0, 0.0, 0.5 ] ] ] ) + } ), + 'Mouse 2' : ScaffoldPackage(MeshType_1d_path1, { + 'scaffoldSettings' : { + 'Coordinate dimensions' : 3, + 'Length' : 1.0, + 'Number of elements' : 4 }, 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ 0.3, -0.2, -0.6 ], [ 0.4, 1.3, -0.1 ], [ -0.1, 0.1, 0.4 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 0.2, 1.4, -0.7 ], [ -1.9, -0.3, -0.5 ], [ -0.1, 0.0, 0.4 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -0.6, -0.3, -0.9 ], [ -2.8, 0.0, -1.4 ], [ -0.3, 0.0, 0.2 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.0, 0.8, -2.2 ], [ -1.8, 0.1, 0.0 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.6, 0.0, -1.6 ], [ 0.7, -1.2, 1.0 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -0.9, -0.9, -1.6 ], [ 1.3, 0.5, -1.6 ], [ 0.0, 0.3, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -0.3, 0.6, -1.7 ], [ 3.1, 0.3, 0.5 ], [ -0.2, 0.0, 0.2 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 0.1, -1.0, -0.7 ], [ -0.3, -0.7, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -0.4, -2.3, -0.1 ], [ -0.4, -0.6, 0.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ] ] ) + [ [ 0.0, 0.0, 0.0 ], [ 0.0, 0.0, 13.0 ], [ 0.0, -10.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 0.0, 0.0, 13.0 ], [ 0.0, 2.0, 28.0 ], [ 0.0, -10.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -14.0, -2.0, 13.0 ], [ 0.0, -3.0, -19.0 ], [ 0.0, -10.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -14.0, -1.0, -10.0 ], [ 1.0, 1.0, -17.0 ], [ 0.0, -10.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -14.0, 0.0, -28.0 ], [ 0.0, 0.0, -11.0 ], [ 0.0, -10.0, 0.0 ], [ 0.0, 0.0, 0.5 ] ] ] ) } ), 'Pig 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, 'Length' : 1.0, - 'Number of elements' : 43 + 'Number of elements' : 36 }, 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ 8.1, 2.4, 5.1 ], [ -1.5, 0.6, -3.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 6.2, 3.0, 1.8 ], [ -3.0, 0.6, -2.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.2, 3.6, 0.5 ], [ -5.5, 1.0, -0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.0, 3.5, 0.8 ], [ -3.6, -2.1, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -4.0, 0.0, 1.3 ], [ 0.0, -4.2, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.0, -3.5, 1.7 ], [ 3.6, -2.1, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.0, -3.5, 2.1 ], [ 3.6, 2.1, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 4.0, 0.0, 2.5 ], [ 0.0, 4.2, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.3, 3.9, 3.1 ], [ -4.1, 2.4, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.2, 3.9, 3.5 ], [ -4.1, -2.4, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -4.5, 0.0, 3.9 ], [ 0.0, -4.7, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.2, -3.9, 4.4 ], [ 4.1, -2.4, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.3, -3.9, 4.8 ], [ 4.1, 2.4, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 4.5, 0.0, 5.1 ], [ 0.3, 3.6, 0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.8, 4.8, 5.5 ], [ -5.1, 2.8, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.7, 4.8, 6.5 ], [ -5.0, -2.9, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -5.5, 0.0, 6.3 ], [ 0.0, -5.8, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.8, -4.8, 6.7 ], [ 5.0, -2.9, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.7, -4.8, 7.1 ], [ 5.0, 2.9, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 5.5, 0.0, 7.5 ], [ 0.0, 5.8, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 3.0, 4.2, 8.2 ], [ -2.5, 1.4, 0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.3, 4.9, 8.8 ], [ -5.0, -1.3, 0.0 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -4.7, 1.5, 8.8 ], [ -1.4, -5.1, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -3.9, -2.5, 9.5 ], [ 4.4, -0.8, 0.8 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.8, 0.3, 10.0 ], [ 0.0, 3.6, 0.9 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -0.4, 3.4, 10.5 ], [ 4.1, 2.7, 0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 3.1, 3.0, 10.3 ], [ 2.2, -1.4, -0.2 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 4.4, -0.1, 9.6 ], [ -0.6, -4.6, -1.0 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.0, -3.1, 8.5 ], [ -2.7, -1.4, -0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.6, -3.6, 7.6 ], [ -2.8, 1.5, -0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -3.9, -0.5, 7.0 ], [ -0.2, 3.5, -0.2 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -2.5, 2.6, 7.0 ], [ 1.8, 1.7, -0.2 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 1.4, 3.8, 6.8 ], [ 5.3, -1.9, -0.8 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 3.6, 1.1, 6.5 ], [ 0.2, -2.4, -0.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 1.9, -2.8, 6.2 ], [ -4.1, -3.7, -0.8 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.7, -3.5, 5.7 ], [ -2.0, 1.0, -0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -3.7, -0.2, 5.1 ], [ 0.2, 4.4, -0.8 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.9, 3.2, 4.6 ], [ 2.3, 1.1, -0.6 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 1.3, 3.6, 4.4 ], [ 3.7, -1.5, -0.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 3.3, 0.9, 4.0 ], [ -0.4, -2.4, -0.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 2.0, -1.9, 3.6 ], [ -2.7, -1.5, -0.2 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.2, -2.7, 3.1 ], [ -3.0, 0.8, -0.4 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ -1.9, 1.5, 2.2 ], [ 5.5, 4.5, -1.9 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ], - [ [ 14.3, -5.9, -3.0 ], [ 5.2, -4.1, -0.3 ], [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.05 ] ] ] ) + [ [ 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 ] ] ] ) } ), 'Test Line' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { @@ -141,8 +147,8 @@ class MeshType_3d_colon1(Scaffold_base): }, 'meshEdits' : zinc_utils.exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ -4.0, 1.0, 3.0 ], [ 5.0, 1.0, -3.0 ], [ 0.8, -4.0, 0.0 ], [0.0, 0.0, 0.0 ] ], - [ [ 1.0, 2.0, 0.0 ], [ 5.0, 1.0, -3.0 ], [ 0.8, -4.0, 0.0 ], [0.0, 0.0, 0.0 ] ] ]) + [ [ -40.0, 10.0, 30.0 ], [ 50.0, 10.0, -30.0 ], [ 8.0, -40.0, 0.0 ], [0.0, 0.0, 0.0 ] ], + [ [ 10.0, 20.0, 0.0 ], [ 50.0, 10.0, -30.0 ], [ 8.0, -40.0, 0.0 ], [0.0, 0.0, 0.0 ] ] ]) } ) } @@ -157,35 +163,77 @@ def getParameterSetNames(): 'Human 1', 'Human 2', 'Mouse 1', + 'Mouse 2', 'Pig 1'] @classmethod def getDefaultOptions(cls, parameterSetName='Default'): if 'Human 2' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 2'] - elif 'Mouse' in parameterSetName: + elif 'Mouse 1' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 1'] + elif 'Mouse 2' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 2'] elif 'Pig' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 1'] else: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 1'] if 'Mouse' in parameterSetName: - segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegmentsimplemesentery1, defaultParameterSetName = 'Mouse 1') + segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Mouse 1') elif 'Pig' in parameterSetName: - segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegmentteniacoli1, defaultParameterSetName = 'Pig 1') + segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Pig 1') else: - segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegmentteniacoli1, defaultParameterSetName = 'Human 1') + segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Human 1') options = { 'Central path' : copy.deepcopy(centralPathOption), 'Segment profile' : segmentProfileOption, 'Number of segments': 30, + 'Proximal length': 420.0, + 'Transverse length': 460.0, + 'Distal length': 620.0, + 'Proximal inner radius': 43.5, + 'Proximal tenia coli width': 10.0, + 'Proximal-transverse inner radius': 33.0, + 'Proximal-transverse tenia coli width': 10.0, + 'Transverse-distal inner radius': 29.0, + 'Transverse-distal tenia coli width': 10.0, + 'Distal inner radius': 31.5, + 'Distal tenia coli width': 10.0, 'Refine' : False, 'Refine number of elements around' : 1, 'Refine number of elements along' : 1, 'Refine number of elements through wall' : 1 } - if 'Pig' in parameterSetName: + if 'Human 2' in parameterSetName: + options['Proximal length'] = 180.0 + options['Transverse length'] = 620.0 + options['Distal length'] = 700.0 + elif 'Mouse' in parameterSetName: + options['Number of segments'] = 10 + options['Proximal length'] = 30.0 + options['Transverse length'] = 20.0 + options['Distal length'] = 25.0 + options['Proximal inner radius'] = 1.0 + options['Proximal tenia coli width'] = 0.5 + options['Proximal-transverse inner radius'] = 0.9 + options['Proximal-transverse tenia coli width'] = 0.7 + options['Transverse-distal inner radius'] = 0.7 + options['Transverse-distal tenia coli width'] = 1.0 + options['Distal inner radius'] = 0.7 + options['Distal tenia coli width'] = 1.0 + elif 'Pig' in parameterSetName: options['Number of segments'] = 120 + options['Proximal length'] = 2610.0 + options['Transverse length'] = 200.0 + options['Distal length'] = 200.0 + options['Proximal inner radius'] = 20.0 + options['Proximal tenia coli width'] = 5.0 + options['Proximal-transverse inner radius'] = 10.5 + options['Proximal-transverse tenia coli width'] = 5.0 + options['Transverse-distal inner radius'] = 8.0 + options['Transverse-distal tenia coli width'] = 5.0 + options['Distal inner radius'] = 8.0 + options['Distal tenia coli width'] = 5.0 return options @staticmethod @@ -194,6 +242,17 @@ def getOrderedOptionNames(): 'Central path', 'Segment profile', 'Number of segments', + 'Proximal length', + 'Transverse length', + 'Distal length', + 'Proximal inner radius', + 'Proximal tenia coli width', + 'Proximal-transverse inner radius', + 'Proximal-transverse tenia coli width', + 'Transverse-distal inner radius', + 'Transverse-distal tenia coli width', + 'Distal inner radius', + 'Distal tenia coli width', 'Refine', 'Refine number of elements around', 'Refine number of elements along', @@ -204,8 +263,7 @@ def getOptionValidScaffoldTypes(cls, optionName): if optionName == 'Central path': return [ MeshType_1d_path1 ] if optionName == 'Segment profile': - return [ MeshType_3d_colonsegmentsimplemesentery1, - MeshType_3d_colonsegmentteniacoli1 ] + return [ MeshType_3d_colonsegment1 ] return [] @classmethod @@ -248,6 +306,20 @@ def checkOptions(cls, options): 'Refine number of elements through wall']: if options[key] < 1: options[key] = 1 + for key in [ + 'Proximal length', + 'Transverse length', + 'Distal length', + 'Proximal inner radius', + 'Proximal tenia coli width', + 'Proximal-transverse inner radius', + 'Proximal-transverse tenia coli width', + 'Transverse-distal inner radius', + 'Transverse-distal tenia coli width', + 'Distal inner radius', + 'Distal tenia coli width']: + if options[key] < 0.0: + options[key] = 0.0 @staticmethod def generateBaseMesh(region, options): @@ -260,37 +332,53 @@ def generateBaseMesh(region, options): centralPath = options['Central path'] segmentProfile = options['Segment profile'] segmentCount = options['Number of segments'] + proximalLength = options['Proximal length'] + transverseLength = options['Transverse length'] + distalLength = options['Distal length'] + proximalInnerRadius = options['Proximal inner radius'] + proximalTCWidth = options['Proximal tenia coli width'] + proximalTransverseInnerRadius = options['Proximal-transverse inner radius'] + proximalTransverseTCWidth = options['Proximal-transverse tenia coli width'] + transverseDistalInnerRadius = options['Transverse-distal inner radius'] + transverseDistalTCWidth = options['Transverse-distal tenia coli width'] + distalInnerRadius = options['Distal inner radius'] + distalTCWidth = options['Distal tenia coli width'] segmentScaffoldType = segmentProfile.getScaffoldType() segmentSettings = segmentProfile.getScaffoldSettings() - if segmentScaffoldType == MeshType_3d_colonsegmentsimplemesentery1: - elementsCountAroundMZ = segmentSettings['Number of elements around mesenteric zone'] - elementsCountAroundNonMZ = segmentSettings['Number of elements around non-mesenteric zone'] - elementsCountAround = elementsCountAroundMZ + elementsCountAroundNonMZ - mzWidth = segmentSettings['Mesenteric zone width'] - else: # segmentScaffoldType == MeshType_3d_colonsegmentteniacoli1: - elementsCountAroundTC = segmentSettings['Number of elements around tenia coli'] - elementsCountAroundHaustrum = segmentSettings['Number of elements around haustrum'] - cornerInnerRadiusFactor = segmentSettings['Corner inner radius factor'] - haustrumInnerRadiusFactor = segmentSettings['Haustrum inner radius factor'] - segmentLengthEndDerivativeFactor = segmentSettings['Segment length end derivative factor'] - segmentLengthMidDerivativeFactor = segmentSettings['Segment length mid derivative factor'] - tcCount = segmentSettings['Number of tenia coli'] - tcWidth = segmentSettings['Tenia coli width'] - tcThickness = segmentSettings['Tenia coli thickness'] - elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount + elementsCountAroundTC = segmentSettings['Number of elements around tenia coli'] + elementsCountAroundHaustrum = segmentSettings['Number of elements around haustrum'] + cornerInnerRadiusFactor = segmentSettings['Corner inner radius factor'] + haustrumInnerRadiusFactor = segmentSettings['Haustrum inner radius factor'] + segmentLengthEndDerivativeFactor = segmentSettings['Segment length end derivative factor'] + segmentLengthMidDerivativeFactor = segmentSettings['Segment length mid derivative factor'] + tcCount = segmentSettings['Number of tenia coli'] + tcThickness = segmentSettings['Tenia coli thickness'] + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount elementsCountAlongSegment = segmentSettings['Number of elements along segment'] elementsCountThroughWall = segmentSettings['Number of elements through wall'] - radius = segmentSettings['Inner radius'] + startRadius = segmentSettings['Start inner radius'] + startRadiusDerivative = segmentSettings['Start inner radius derivative'] + endRadius = segmentSettings['End inner radius'] + endRadiusDerivative = segmentSettings['End inner radius derivative'] wallThickness = segmentSettings['Wall thickness'] useCrossDerivatives = segmentSettings['Use cross derivatives'] useCubicHermiteThroughWall = not(segmentSettings['Use linear through wall']) elementsCountAlong = int(elementsCountAlongSegment*segmentCount) + firstNodeIdentifier = 1 + firstElementIdentifier = 1 + + # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion) + # for i in range(len(cx)): + # print('cx = ', i+1, cx[i]) + # print('cd1 = ', i+1, cd1[i]) + # print('cd2 = ', i+1, cd2[i]) + # print('cd12 = ', i+1, cd12[i]) del tmpRegion # find arclength of colon @@ -300,32 +388,111 @@ def generateBaseMesh(region, options): magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) for e in range(elementsCountIn): arcLength = interp.getCubicHermiteArcLength(cx[e], sd1[e], cx[e + 1], sd1[e + 1]) + # print(e+1, arcLength) length += arcLength segmentLength = length / segmentCount + # print('Length = ', length) + + # Sample central path + sx, sd1, se, sxi, ssf = interp.sampleCubicHermiteCurves(cx, cd1, elementsCountAlongSegment*segmentCount) + sd2 = interp.interpolateSampleCubicHermite(cd2, cd12, se, sxi, ssf)[0] + + # Generate variation of radius & tc width along length + lengthList = [0.0, proximalLength, proximalLength + transverseLength, length] + innerRadiusList = [proximalInnerRadius, proximalTransverseInnerRadius, transverseDistalInnerRadius, distalInnerRadius] + innerRadiusSegmentList, dInnerRadiusSegmentList = interp.sampleParameterAlongLine(lengthList, innerRadiusList, segmentCount) + + tcWidthList = [proximalTCWidth, proximalTransverseTCWidth, transverseDistalTCWidth, distalTCWidth] + tcWidthSegmentList, dTCWidthSegmentList = interp.sampleParameterAlongLine(lengthList, tcWidthList, segmentCount) + + xExtrude = [] + d1Extrude = [] + d2Extrude = [] + d3UnitExtrude = [] + + # Create object + colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( + region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, + tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, + segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + innerRadiusSegmentList, dInnerRadiusSegmentList, tcWidthSegmentList, dTCWidthSegmentList) + + for nSegment in range(segmentCount): + # Create inner points + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray = \ + colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) + + # Warp segment points + xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = tubemesh.warpSegmentPoints( + xInner, d1Inner, d2Inner, segmentAxis, segmentLength, sx, sd1, sd2, + elementsCountAround, elementsCountAlongSegment, nSegment) + + # Store points along length + xExtrude = xExtrude + (xWarpedList if nSegment == 0 else xWarpedList[elementsCountAround:]) + d1Extrude = d1Extrude + (d1WarpedList if nSegment == 0 else d1WarpedList[elementsCountAround:]) + + # Smooth d2 for nodes between segments and recalculate d3 + if nSegment == 0: + d2Extrude = d2Extrude + (d2WarpedList[:-elementsCountAround]) + d3UnitExtrude = d3UnitExtrude + (d3WarpedUnitList[:-elementsCountAround]) + else: + xSecondFace = xWarpedList[elementsCountAround:elementsCountAround*2] + d1SecondFace = d1WarpedList[elementsCountAround:elementsCountAround*2] + d2SecondFace = d2WarpedList[elementsCountAround:elementsCountAround*2] + for n1 in range(elementsCountAround): + nx = [xLastTwoFaces[n1], xLastTwoFaces[n1 + elementsCountAround], xSecondFace[n1]] + nd2 = [d2LastTwoFaces[n1], d2LastTwoFaces[n1 + elementsCountAround], d2SecondFace[n1]] + d2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixStartDerivative = True, fixEndDerivative = True)[1] + d2Extrude.append(d2) + d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1LastTwoFaces[n1 + elementsCountAround]), vector.normalise(d2))) + d3UnitExtrude.append(d3Unit) + d2Extrude = d2Extrude + (d2WarpedList[elementsCountAround:-elementsCountAround] if nSegment < segmentCount - 1 else d2WarpedList[elementsCountAround:]) + d3UnitExtrude = d3UnitExtrude + (d3WarpedUnitList[elementsCountAround:-elementsCountAround] if nSegment < segmentCount - 1 else d3WarpedUnitList[elementsCountAround:]) + xLastTwoFaces = xWarpedList[-elementsCountAround*2:] + d1LastTwoFaces = d1WarpedList[-elementsCountAround*2:] + d2LastTwoFaces = d2WarpedList[-elementsCountAround*2:] + + contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() + + # Create coordinates and derivatives + xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xExtrude, d1Extrude, + d2Extrude, d3UnitExtrude, sx, contractedWallThicknessList, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList) + + relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() + + if tcThickness > 0: + tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() + xList, d1List, d2List, d3List, annotationGroups, annotationArray = getTeniaColi( + region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, + elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, + tubeTCWidthList, tcThickness, sx, annotationGroups, annotationArray) + + # Create flat and texture coordinates + xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( + xiList, relaxedLengthList, length, wallThickness, tcCount, tcThickness, + elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, + elementsCountThroughWall, transitElementList) + + # 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, annotationArray, firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives) + + else: + # Create flat and texture coordinates + xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( + xiList, relaxedLengthList, length, wallThickness, elementsCountAround, + elementsCountAlong, elementsCountThroughWall, transitElementList) - # Generate inner surface of a colon segment - if segmentScaffoldType == MeshType_3d_colonsegmentsimplemesentery1: - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength, xInner, d1Inner, d2Inner, segmentAxis = getColonSegmentInnerPointsNoTeniaColi(region, elementsCountAroundMZ, elementsCountAroundNonMZ, - elementsCountAlongSegment, mzWidth, radius, segmentLength, wallThickness) - else: # segmentScaffoldType == MeshType_3d_colonsegmentteniacoli1: - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength, xInner, d1Inner, d2Inner, segmentAxis = getColonSegmentInnerPointsTeniaColi(region, elementsCountAroundTC, elementsCountAroundHaustrum, - elementsCountAlongSegment, tcCount, tcWidth, radius, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, - segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, segmentLength, wallThickness) - - # Generate tube mesh - annotationGroups, nextNodeIdentifier, nextElementIdentifier, xList, d1List, d2List, d3List, sx, curvatureAlong, factorList = tubemesh.generatetubemesh(region, - elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, segmentCount, cx, cd1, cd2, cd12, - xInner, d1Inner, d2Inner, wallThickness, segmentAxis, segmentLength, useCrossDerivatives, useCubicHermiteThroughWall, - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength) - - # Generate tenia coli - if segmentScaffoldType == MeshType_3d_colonsegmentteniacoli1: - annotationGroupsTC, nextNodeIdentifier, nextElementIdentifier = getTeniaColi(region, nextNodeIdentifier, nextElementIdentifier, - useCrossDerivatives, useCubicHermiteThroughWall, xList, d1List, d2List, d3List, elementsCountAroundTC, elementsCountAroundHaustrum, - elementsCountAlong, elementsCountThroughWall, wallThickness, tcWidth, tcThickness, sx, curvatureAlong, factorList, uList, - arcLengthOuterMidLength, length, tcCount) - - annotationGroups += annotationGroupsTC + # 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, + useCubicHermiteThroughWall, useCrossDerivatives) return annotationGroups diff --git a/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py new file mode 100644 index 00000000..04a1ccfb --- /dev/null +++ b/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -0,0 +1,1779 @@ +""" +Generates a single 3-D colon segment mesh along a central +line, with variable numbers of elements around, along and +through wall and number of tenia coli, with variable radius +and thickness along. +""" + +import math +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear +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 +from scaffoldmaker.utils import zinc_utils +from opencmiss.zinc.element import Element, Elementbasis +from opencmiss.zinc.field import Field +from opencmiss.zinc.node import Node + +class MeshType_3d_colonsegment1(Scaffold_base): + ''' + Generates a single 3-D colon segment mesh with variable + numbers of tenia coli, elements around, along the central + line, and through wall. The cross-section profile of the colon + segment varies with species and is dependent on the number + of tenia coli. + Mouse: One "tenia coli", circular profile + Pig: Two tenia coli, bow tie profile + Human (Default): Three tenia coli, triangular profile with + rounded corners at the inter-haustral septa, and a clover + profile in the intra-haustral region. + ''' + @staticmethod + def getName(): + return '3D Colon Segment 1' + + @staticmethod + def getParameterSetNames(): + return [ + 'Default', + 'Human 1', + 'Mouse 1', + 'Pig 1'] + + @staticmethod + def getDefaultOptions(parameterSetName='Default'): + options = { + 'Number of elements around tenia coli' : 2, + 'Number of elements around haustrum' : 8, + 'Number of elements along segment' : 4, + 'Number of elements through wall' : 1, + 'Start inner radius': 43.5, + 'Start inner radius derivative': 0.0, + 'End inner radius': 33.0, + 'End inner radius derivative': 0.0, + 'Corner inner radius factor': 0.5, + 'Haustrum inner radius factor': 0.5, + 'Segment length end derivative factor': 0.5, + 'Segment length mid derivative factor': 3.0, + 'Segment length': 50.0, + 'Number of tenia coli': 3, + 'Start tenia coli width': 10.0, + 'Start tenia coli width derivative': 0.0, + 'End tenia coli width': 10.0, + 'End tenia coli width derivative': 0.0, + 'Tenia coli thickness': 1.6, + 'Wall thickness': 1.6, + 'Use cross derivatives' : False, + 'Use linear through wall' : True, + 'Refine' : False, + 'Refine number of elements around' : 1, + 'Refine number of elements along segment' : 1, + 'Refine number of elements through wall' : 1 + } + if 'Mouse' in parameterSetName: + options['Start inner radius'] = 0.94 + options['End inner radius'] = 0.94 + options['Corner inner radius factor'] = 0.0 + options['Haustrum inner radius factor'] = 0.0 + options['Segment length end derivative factor'] = 0.0 + options['Segment length mid derivative factor'] = 0.0 + options['Number of tenia coli'] = 1 + options['Start tenia coli width'] = 0.8 + options['End tenia coli width'] = 0.8 + options['Tenia coli thickness'] = 0.0 + options['Wall thickness'] = 0.55 + elif 'Pig' in parameterSetName: + 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['Segment length end derivative factor'] = 0.8 + options['Segment length mid derivative factor'] = 2.0 + options['Segment length'] = 25.0 + options['Number of tenia coli'] = 2 + options['Start tenia coli width'] = 5.0 + options['End tenia coli width'] = 5.0 + options['Tenia coli thickness'] = 0.5 + options['Wall thickness'] = 2.0 + return options + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Number of elements around tenia coli', + 'Number of elements around haustrum', + 'Number of elements along segment', + 'Number of elements through wall', + 'Start inner radius', + 'Start inner radius derivative', + 'End inner radius', + 'End inner radius derivative', + 'Corner inner radius factor', + 'Haustrum inner radius factor', + 'Segment length end derivative factor', + 'Segment length mid derivative factor', + 'Segment length', + 'Number of tenia coli', + 'Start tenia coli width', + 'Start tenia coli width derivative', + 'End tenia coli width', + 'End tenia coli width derivative', + 'Tenia coli thickness', + 'Wall thickness', + 'Use cross derivatives', + 'Use linear through wall', + 'Refine', + 'Refine number of elements around', + 'Refine number of elements along segment', + 'Refine number of elements through wall' + ] + + @staticmethod + def checkOptions(options): + for key in [ + 'Number of elements through wall', + 'Refine number of elements around', + 'Refine number of elements along segment', + 'Refine number of elements through wall']: + if options[key] < 1: + options[key] = 1 + for key in [ + 'Number of elements around tenia coli', + 'Number of elements along segment']: + if options[key] < 2: + options[key] = 2 + if options['Number of elements around haustrum'] < 4: + options['Number of elements around haustrum'] = 4 + for key in [ + 'Number of elements around tenia coli', + 'Number of elements around haustrum']: + if options[key] % 2 > 0: + options[key] = options[key] + 1 + for key in [ + 'Start inner radius', + 'End inner radius', + 'Haustrum inner radius factor', + 'Segment length end derivative factor', + 'Segment length mid derivative factor', + 'Segment length', + 'Tenia coli thickness', + 'Wall thickness']: + if options[key] < 0.0: + options[key] = 0.0 + if options['Corner inner radius factor'] < 0.1: + options['Corner inner radius factor'] = 0.1 + for key in [ + 'Corner inner radius factor', + 'Segment length end derivative factor']: + if options[key] > 1.0: + options[key] = 1.0 + if options['Number of tenia coli'] < 1: + options['Number of tenia coli'] = 1 + elif options['Number of tenia coli'] > 3: + options['Number of tenia coli'] = 3 + for key in [ + 'Start tenia coli width', + 'End tenia coli width']: + if options[key] < 0.2*min(options['Start inner radius'], options['End inner radius']): + options[key] = round(0.2*min(options['Start inner radius'], options['End inner radius']), 2) + if options['Start tenia coli width'] > round(math.sqrt(3)*0.5*options['Start inner radius'], 2): + options['Start tenia coli width'] = round(math.sqrt(3)*0.5*options['Start inner radius'], 2) + if options['End tenia coli width'] > round(math.sqrt(3)*0.5*options['End inner radius'], 2): + options['End tenia coli width'] = round(math.sqrt(3)*0.5*options['End inner radius'], 2) + + @staticmethod + def generateBaseMesh(region, options): + """ + Generate the base tricubic Hermite mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: None + """ + elementsCountAroundTC = options['Number of elements around tenia coli'] + elementsCountAroundHaustrum = options['Number of elements around haustrum'] + elementsCountAlongSegment = options['Number of elements along segment'] + elementsCountThroughWall = options['Number of elements through wall'] + startRadius = options['Start inner radius'] + startRadiusDerivative = options['Start inner radius derivative'] + endRadius = options['End inner radius'] + endRadiusDerivative = options['End inner radius derivative'] + cornerInnerRadiusFactor = options['Corner inner radius factor'] + haustrumInnerRadiusFactor = options['Haustrum inner radius factor'] + segmentLengthEndDerivativeFactor = options['Segment length end derivative factor'] + segmentLengthMidDerivativeFactor = options['Segment length mid derivative factor'] + segmentLength = options['Segment length'] + tcCount = options['Number of tenia coli'] + startTCWidth = options['Start tenia coli width'] + startTCWidthDerivative = options['Start tenia coli width derivative'] + endTCWidth = options['End tenia coli width'] + endTCWidthDerivative = options['End tenia coli width derivative'] + tcThickness = options['Tenia coli thickness'] + wallThickness = options['Wall thickness'] + useCrossDerivatives = options['Use cross derivatives'] + useCubicHermiteThroughWall = not(options['Use linear through wall']) + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount + segmentCount = 1 + firstNodeIdentifier = 1 + firstElementIdentifier = 1 + + # Central path + cx = [ [ 0.0, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] + cd1 = [ [ segmentLength, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] + cd2 = [ [ 0.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ] ] + cd12 = [ [0.0, 0.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ] + + # Sample central path + sx, sd1, se, sxi, ssf = interp.sampleCubicHermiteCurves(cx, cd1, elementsCountAlongSegment*segmentCount) + sd2 = interp.interpolateSampleCubicHermite(cd2, cd12, se, sxi, ssf)[0] + + # Radius and tenia coli width along segment length + radiusList = [startRadius, endRadius] + dRadiusList = [startRadiusDerivative, endRadiusDerivative] + tcWidthList = [startTCWidth, endTCWidth] + dTCWidthList = [startTCWidthDerivative, endTCWidthDerivative] + + colonSegmentTubeMeshInnerPoints = ColonSegmentTubeMeshInnerPoints( + region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, + tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, + segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + radiusList, dRadiusList, tcWidthList, dTCWidthList) + + # Create inner points + nSegment = 0 + xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray = \ + colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment) + + # Warp segment points + xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = tubemesh.warpSegmentPoints( + xInner, d1Inner, d2Inner, segmentAxis, segmentLength, sx, sd1, sd2, + elementsCountAround, elementsCountAlongSegment, nSegment) + + contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList() + + # Create coordinates and derivatives + xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xWarpedList, d1WarpedList, + d2WarpedList, d3WarpedUnitList, sx, contractedWallThicknessList, + elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, transitElementList) + + relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList() + + if tcThickness > 0: + tubeTCWidthList = colonSegmentTubeMeshInnerPoints.getTubeTCWidthList() + xList, d1List, d2List, d3List, annotationGroups, annotationArray = getTeniaColi( + region, xList, d1List, d2List, d3List, curvatureList, tcCount, elementsCountAroundTC, + elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, + tubeTCWidthList, tcThickness, sx, annotationGroups, annotationArray) + + # Create flat and texture coordinates + xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = createFlatAndTextureCoordinatesTeniaColi( + xiList, relaxedLengthList, segmentLength, wallThickness, tcCount, tcThickness, + elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, + elementsCountThroughWall, transitElementList) + + # 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, annotationArray, firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives) + + else: + # Create flat and texture coordinates + xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture = tubemesh.createFlatAndTextureCoordinates( + xiList, relaxedLengthList, segmentLength, wallThickness, elementsCountAround, + elementsCountAlongSegment, elementsCountThroughWall, transitElementList) + + # Create nodes and elements + nextNodeIdentifier, nextElementIdentifier, annotationGroups = tubemesh.createNodesAndElements( + region, xList, d1List, d2List, d3List, xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture, + elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, + annotationGroups, annotationArray, firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives) + + return annotationGroups + + @classmethod + def generateMesh(cls, region, options): + """ + Generate base or refined mesh. + :param region: Zinc region to create mesh in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + """ + if not options['Refine']: + return cls.generateBaseMesh(region, options) + + refineElementsCountAround = options['Refine number of elements around'] + refineElementsCountAlong = options['Refine number of elements along segment'] + refineElementsCountThroughWall = options['Refine number of elements through wall'] + + baseRegion = region.createRegion() + baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) + + meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) + return meshrefinement.getAnnotationGroups() + +class ColonSegmentTubeMeshInnerPoints: + """ + Generates inner profile of a colon segment for use by tubemesh. + """ + + def __init__(self, region, elementsCountAroundTC, elementsCountAroundHaustrum, + elementsCountAlongSegment, tcCount, segmentLengthEndDerivativeFactor, + segmentLengthMidDerivativeFactor, segmentLength, wallThickness, + cornerInnerRadiusFactor, haustrumInnerRadiusFactor, innerRadiusSegmentList, + dInnerRadiusSegmentList, tcWidthSegmentList, dTCWidthSegmentList): + + self._region = region + self._elementsCountAroundTC = elementsCountAroundTC + self._elementsCountAroundHaustrum = elementsCountAroundHaustrum + self._elementsCountAlongSegment = elementsCountAlongSegment + self._tcCount = tcCount + self._segmentLengthEndDerivativeFactor = segmentLengthEndDerivativeFactor + self._segmentLengthMidDerivativeFactor = segmentLengthMidDerivativeFactor + self._segmentLength = segmentLength + self._wallThickness = wallThickness + self._cornerInnerRadiusFactor = cornerInnerRadiusFactor + self._haustrumInnerRadiusFactor = haustrumInnerRadiusFactor + self._innerRadiusSegmentList = innerRadiusSegmentList + self._dInnerRadiusSegmentList = dInnerRadiusSegmentList + self._tcWidthSegmentList = tcWidthSegmentList + self._dTCWidthSegmentList = dTCWidthSegmentList + self._tubeTCWidthList = [] + self._xiList = [] + self._relaxedLengthList = [] + self._contractedWallThicknessList = [] + + def getColonSegmentTubeMeshInnerPoints(self, nSegment): + + # Unpack radius and rate of change of inner radius + startRadius = self._innerRadiusSegmentList[nSegment] + startRadiusDerivative = self._dInnerRadiusSegmentList[nSegment] + endRadius = self._innerRadiusSegmentList[nSegment+1] + endRadiusDerivative = self._dInnerRadiusSegmentList[nSegment+1] + + # Unpack tcWidth and rate of change of tcWidth + startTCWidth = self._tcWidthSegmentList[nSegment] + startTCWidthDerivative = self._dTCWidthSegmentList[nSegment] + endTCWidth = self._tcWidthSegmentList[nSegment+1] + endTCWidthDerivative = self._dTCWidthSegmentList[nSegment+1] + + xInner, d1Inner, d2Inner, transitElementList, xiSegment, relaxedLengthSegment, \ + contractedWallThicknessSegment, sTCWidth, 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, + startRadius, startRadiusDerivative, endRadius, endRadiusDerivative, + startTCWidth, startTCWidthDerivative, endTCWidth, endTCWidthDerivative) + + startIdx = 0 if nSegment == 0 else 1 + for i in range(startIdx, self._elementsCountAlongSegment + 1): + self._tubeTCWidthList.append(sTCWidth[i]) + + xi = xiSegment[startIdx:self._elementsCountAlongSegment + 1] + self._xiList += xi + + relaxedLength = relaxedLengthSegment[startIdx:self._elementsCountAlongSegment + 1] + self._relaxedLengthList += relaxedLength + + contractedWallThickness = contractedWallThicknessSegment[startIdx:self._elementsCountAlongSegment + 1] + self._contractedWallThicknessList += contractedWallThickness + + return xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray + + def getTubeTCWidthList(self): + return self._tubeTCWidthList + + def getRelaxedLengthAndXiList(self): + return self._relaxedLengthList, self._xiList + + def getContractedWallThicknessList(self): + return self._contractedWallThicknessList + +def getColonSegmentInnerPoints(region, elementsCountAroundTC, + elementsCountAroundHaustrum, elementsCountAlongSegment, + tcCount, segmentLengthEndDerivativeFactor, + segmentLengthMidDerivativeFactor, segmentLength, wallThickness, + cornerInnerRadiusFactor, haustrumInnerRadiusFactor, + startRadius, startRadiusDerivative, endRadius, endRadiusDerivative, + startTCWidth, startTCWidthDerivative, endTCWidth, endTCWidthDerivative): + """ + Generates a 3-D colon segment mesh with variable numbers of tenia coli, + numbers of elements around, along the central path, and through wall. + Colon segment with one "tenia coli" (mouse) has a circular profile along + its length. Colon segment with two tenia coli (pig) has a circular profile + at the inter-haustral septa, and a bowtie profile in the intra-haustral + region. Colon segment with three tenia coli (human) has a triangular profile + with rounded corners at the inter-haustral septa, and a clover profile + in the intra-haustral region. + :param elementsCountAroundTC: Number of elements around each tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param elementsCountAlongSegment: Number of elements along colon segment. + :param tcCount: Number of tenia coli. + :param segmentLengthEndDerivativeFactor: Factor is multiplied by segment + length to scale derivative along the end of a segment length. + :param segmentLengthMidDerivativeFactor: Factor is multiplied by segment + length to scale derivative along the mid length of the segment. + :param segmentLength: Length of a colon segment. + :param wallThickness: Thickness of wall. + :param cornerInnerRadiusFactor: Roundness of triangular corners of + 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 + radius to obtain radius of intersecting circles in the middle cross-section + along a haustra segment. + :param startRadius: Inner radius defined from center of triangular + profile to vertex of the triangle at proximal end of the colon segment. + :param startRadiusDerivative: Rate of change of inner radius at proximal end. + :param endRadius: Inner radius defined from center of triangular + profile to vertex of the triangle at distal end of the colon segment. + :param endRadiusDerivative: Rate of change of inner radius at distal end. + :param startTCWidth: Width of tenia coli at proximal end of the colon segment. + :param startTCWidthDerivative: Rate of change of tenia coli width at proximal end. + :param endTCWidth: Width of tenia coli at distal end of the colon segment. + :param endTCWidthDerivative: Rate of change of tenia coli width at distal end. + :return coordinates, derivatives on inner surface of a colon segment. + :return transitElementList: stores true if element around is an element that + transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. + :return xiList: List of xi for each node around. xi refers to node position + along the relaxed length when colon segment is opened into a flat preparation, + nominally in [0.0, 1.0]. + :return relaxedLengthList: List of relaxed length around elements for each element + along colon segment when the segment is opened into a flat preparation. + :return contractedWallThicknessList: List of wall thickness for each element + along colon segment. Assume incompressiblity and a shortened length around will + result in a thicker wall and vice-versa. + :return sTCWidthAlongSegment: width of tenia coli along segment. + :return segmentAxis: Axis of segment. + :return annotationGroups, annotationArray: annotationArray stores annotation + names of elements around. + """ + + transitElementListHaustrum = ([0]*int(elementsCountAroundTC*0.5) + [1] + + [0]*int(elementsCountAroundHaustrum - 2) + [1] + + [0]*int(elementsCountAroundTC*0.5)) + transitElementList = transitElementListHaustrum * tcCount + + # Find parameter variation along elementsCountAlongSegment + sRadiusAlongSegment = [] + sTCWidthAlongSegment = [] + for n2 in range(elementsCountAlongSegment + 1): + xi = 1/elementsCountAlongSegment * n2 + radius = interp.interpolateCubicHermite([startRadius], [startRadiusDerivative], [endRadius], [endRadiusDerivative], xi)[0] + sRadiusAlongSegment.append(radius) + tcWidth = interp.interpolateCubicHermite([startTCWidth], [startTCWidthDerivative], [endTCWidth], [endTCWidthDerivative], xi)[0] + sTCWidthAlongSegment.append(tcWidth) + + # create nodes + x = [ 0.0, 0.0, 0.0 ] + dx_ds1 = [ 0.0, 0.0, 0.0 ] + dx_ds2 = [ 0.0, 0.0, 0.0 ] + sampleElementOut = 20 + segmentAxis = [0.0, 0.0, 1.0] + + d2HalfSet = [] + d2Raw = [] + xInnerRaw = [] + dx_ds2InnerRaw = [] + xFinal = [] + d1Final = [] + d2Final = [] + xiList = [] + relaxedLengthList = [] + contractedWallThicknessList = [] + + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount + if tcCount == 1: + for n2 in range(elementsCountAlongSegment + 1): + xHalfSet, d1HalfSet = createHalfSetInterHaustralSegment(elementsCountAroundTC, + elementsCountAroundHaustrum, tcCount, sTCWidthAlongSegment[n2], + sRadiusAlongSegment[n2], cornerInnerRadiusFactor, sampleElementOut) + for i in range(len(xHalfSet)): + d2HalfSet.append([0.0, 0.0, 0.0]) + x, d1, _ = getFullProfileFromHalfHaustrum(xHalfSet, d1HalfSet, d2HalfSet, tcCount) + z = segmentLength / elementsCountAlongSegment * n2 + d1Final = d1Final + d1 + xFace = [] + for n1 in range(elementsCountAroundTC + elementsCountAroundHaustrum): + xFinal.append([x[n1][0], x[n1][1], z]) + xFace.append([x[n1][0], x[n1][1], z]) + xiFace, lengthAroundFace = getXiListFromOuterLengthProfile(xFace, d1, segmentAxis, + wallThickness, transitElementList) + xiList.append(xiFace) + relaxedLengthList.append(lengthAroundFace) + contractedWallThicknessList.append(wallThickness) + + for n1 in range(elementsCountAround): + xUp = [] + d2Up = [] + for n2 in range(elementsCountAlongSegment + 1): + n = elementsCountAround * n2 + n1 + xUp.append(xFinal[n]) + d2 = [ xFinal[n + elementsCountAround][i] - xFinal[n][i] if n2 < elementsCountAlongSegment else xFinal[n][i] - xFinal[n - elementsCountAround][i] for i in range(3)] + d2Up.append(d2) + d2Smoothed = interp.smoothCubicHermiteDerivativesLine(xUp, d2Up) + d2Raw.append(d2Smoothed) + + # Re-arrange d2Raw + for n2 in range(elementsCountAlongSegment + 1): + for n1 in range(elementsCountAround): + d2Final.append(d2Raw[n1][n2]) + + # Create annotation groups for mouse colon + mzGroup = AnnotationGroup(region, 'mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + nonmzGroup = AnnotationGroup(region, 'non-mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + annotationGroups = [mzGroup, nonmzGroup] + annotationArray = (['mesenteric zone']*int(elementsCountAroundTC*0.5) + + ['non-mesenteric zone']*elementsCountAroundHaustrum + + ['mesenteric zone']*int(elementsCountAroundTC*0.5)) + + else: + # Find parameter variation at start, mid and end of segment + sRadius = [] + sdRadius = [] + sTCWidth = [] + for n2 in range(3): + xi = 1/2 * n2 + radius = interp.interpolateCubicHermite([startRadius], [startRadiusDerivative], [endRadius], [endRadiusDerivative], xi)[0] + sRadius.append(radius) + dRadius = interp.interpolateCubicHermiteDerivative([startRadius], [startRadiusDerivative], [endRadius], [endRadiusDerivative], xi)[0] + sdRadius.append(dRadius) + tcWidth = interp.interpolateCubicHermite([startTCWidth], [startTCWidthDerivative], [endTCWidth], [endTCWidthDerivative], xi)[0] + sTCWidth.append(tcWidth) + + # Calculate x and d1 at the start, mid and end faces + xHalfSetStart, d1HalfSetStart = createHalfSetInterHaustralSegment( + elementsCountAroundTC, elementsCountAroundHaustrum, tcCount, sTCWidth[0], sRadius[0], + cornerInnerRadiusFactor, sampleElementOut) + + xHalfSetMid, d1HalfSetMid = createHalfSetIntraHaustralSegment( + elementsCountAroundTC, elementsCountAroundHaustrum, tcCount, sTCWidth[1], sRadius[1], + cornerInnerRadiusFactor, sampleElementOut, haustrumInnerRadiusFactor) + + xHalfSetEnd, d1HalfSetEnd = createHalfSetInterHaustralSegment( + elementsCountAroundTC, elementsCountAroundHaustrum, tcCount, sTCWidth[2], sRadius[2], + cornerInnerRadiusFactor, sampleElementOut) + + # Sample arclength of haustra segment + elementsCountAroundHalfHaustrum = int((elementsCountAroundTC + elementsCountAroundHaustrum)*0.5) + + for n1 in range(elementsCountAroundHalfHaustrum + 1): + radiansAround = (math.pi*2/(tcCount*2)) / elementsCountAroundHalfHaustrum * n1 + d2Start = [ sdRadius[0]*0.5, 0.0, segmentLength*0.5 ] + d2Mid = [ sdRadius[1]*0.5, 0.0, segmentLength*0.5 ] + d2End = [ sdRadius[2]*0.5, 0.0, segmentLength*0.5 ] + + if n1 > elementsCountAroundTC*0.5: + # Non tenia coli + # Rotate about segment axis (z-axis) + d2StartRot = matrix.rotateAboutZAxis(d2Start, radiansAround) + d2MidRot = matrix.rotateAboutZAxis(d2Mid, radiansAround) + d2EndRot = matrix.rotateAboutZAxis(d2End, radiansAround) + + d1 = [ c*segmentLengthEndDerivativeFactor for c in d2StartRot ] + d2 = [ c*segmentLengthMidDerivativeFactor for c in d2MidRot ] + d3 = [ c*segmentLengthEndDerivativeFactor for c in d2EndRot ] + + else: + # Tenia coli do not have haustra so not subjected to derivative scaling along segment + d1 = d2Start + d2 = d2Mid + d3 = d2End + + v1 = [xHalfSetStart[n1][0], xHalfSetStart[n1][1], 0.0] + v2 = [xHalfSetMid[n1][0], xHalfSetMid[n1][1], segmentLength/2] + v3 = [xHalfSetEnd[n1][0], xHalfSetEnd[n1][1], segmentLength] + + nx = [v1, v2, v3] + nd1 = [d1, d2, d3] + + sx, sd1, se, sxi, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAlongSegment) + xInnerRaw.append(sx) + dx_ds2InnerRaw.append(sd1) + + # Re-arrange sample order & calculate dx_ds1 and dx_ds3 from dx_ds2 + for n2 in range(elementsCountAlongSegment + 1): + xAround = [] + unitdx_ds1Around = [] + d2Around = [] + + for n1 in range(elementsCountAroundHalfHaustrum+1): + x = xInnerRaw[n1][n2] + dx_ds2 = dx_ds2InnerRaw[n1][n2] + unitTangent = vector.normalise(dx_ds2) + + # Intra-Haustra segments + if n1 == 0: + unitdx_ds1 = vector.normalise(d1HalfSetStart[n1]) + else: # points on clover + if n2 <= int(elementsCountAlongSegment/2): # first half of segmentLength + axisRot = vector.crossproduct3(segmentAxis, unitTangent) + elif n2 > int(elementsCountAlongSegment/2): # second half of segmentLength + axisRot = vector.crossproduct3(unitTangent, segmentAxis) + rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, math.pi/2) + rotNormal = [rotFrame[j][0]*unitTangent[0] + rotFrame[j][1]*unitTangent[1] + rotFrame[j][2]*unitTangent[2] for j in range(3)] + unitdx_ds3 = vector.normalise(rotNormal) + unitdx_ds1 = vector.crossproduct3(unitTangent, unitdx_ds3) + + xAround.append(x) + d2Around.append(dx_ds2) + unitdx_ds1Around.append(unitdx_ds1) + + if n2 > 0 and n2 < elementsCountAlongSegment: + dx_ds1InnerAroundList = [] + if elementsCountAlongSegment%2 == 0 and n2 == int(elementsCountAlongSegment*0.5): + dx_ds1InnerAroundList = dx_ds1InnerAroundList + d1HalfSetMid + else: + for n1 in range(elementsCountAroundHalfHaustrum): + v1 = xAround[n1] + d1 = unitdx_ds1Around[n1] + v2 = xAround[n1+1] + d2 = unitdx_ds1Around[n1+1] + arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) + dx_ds1 = [c*arcLengthAround for c in d1] + dx_ds1InnerAroundList.append(dx_ds1) + # Account for d1 of node sitting on half haustrum + dx_ds1 = [c*arcLengthAround for c in unitdx_ds1Around[elementsCountAroundHalfHaustrum]] + dx_ds1InnerAroundList.append(dx_ds1) + d1Smoothed = interp.smoothCubicHermiteDerivativesLine(xAround, dx_ds1InnerAroundList, fixStartDerivative = True) + d1TCEdge = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC*0.5)], vector.magnitude(d1Smoothed[int(elementsCountAroundTC*0.5 - 1)])) + d1Transition = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC*0.5 + 1)], vector.magnitude(d1Smoothed[int(elementsCountAroundTC*0.5 + 2)])) + d1Corrected = [] + d1Corrected = d1Corrected + d1Smoothed[:int(elementsCountAroundTC*0.5)] + d1Corrected.append(d1TCEdge) + d1Corrected.append(d1Transition) + d1Corrected = d1Corrected + d1Smoothed[int(elementsCountAroundTC*0.5 + 2):] + elif n2 < 1: + d1Corrected = d1HalfSetStart + elif n2 > elementsCountAlongSegment - 1: + d1Corrected = d1HalfSetEnd + + xAlongList, d1AlongList, d2AlongList = getFullProfileFromHalfHaustrum(xAround, d1Corrected, d2Around, tcCount) + + # Calculate xiList for relaxed state + xHalfSetRelaxed, d1HalfSetRelaxed = createHalfSetIntraHaustralSegment(elementsCountAroundTC, + elementsCountAroundHaustrum, tcCount, sTCWidthAlongSegment[n2], sRadiusAlongSegment[n2], + cornerInnerRadiusFactor, sampleElementOut, haustrumInnerRadiusFactor) + xRelaxed, d1Relaxed, _ = getFullProfileFromHalfHaustrum(xHalfSetRelaxed, d1HalfSetRelaxed, d2Around, tcCount) + xiFace, relaxedLengthAroundFace = getXiListFromOuterLengthProfile(xRelaxed, d1Relaxed, segmentAxis, + wallThickness, transitElementList) + xiList.append(xiFace) + relaxedLengthList.append(relaxedLengthAroundFace) + + contractedLengthAroundFace = getXiListFromOuterLengthProfile(xAlongList, d1AlongList, segmentAxis, wallThickness, transitElementList)[1] + contractedWallThickness = relaxedLengthAroundFace * wallThickness / contractedLengthAroundFace + contractedWallThicknessList.append(contractedWallThickness) + + xFinal = xFinal + xAlongList + d1Final = d1Final + d1AlongList + d2Final = d2Final + d2AlongList + + # Create annotation groups + annotationGroups = [] + annotationArray = ['']*(elementsCountAround) + + return xFinal, d1Final, d2Final, transitElementList, xiList, relaxedLengthList, contractedWallThicknessList, sTCWidthAlongSegment, segmentAxis, annotationGroups, annotationArray + +def createHalfSetInterHaustralSegment(elementsCountAroundTC, elementsCountAroundHaustrum, + tcCount, tcWidth, radius, cornerInnerRadiusFactor, sampleElementOut): + """ + Find locations and derivative of nodes in half of an + inter-haustral segment. Circular profile for segment + with two or less tenia coli and triangular profile + with round corners for segment with three tenia coli. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param tcCount: Number of tenia coli. + :param tcWidth: Width of tenia coli. + :param radius: Inner radius of circular profile with two tenia coli, + radius of circle enclosing triangle for profile with three tenia coli. + :param cornerInnerRadiusFactor: Roundness of triangular corners of + 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 sampleElementOut: Number of sample points used to set up profile + :return: Node location and derivative on half of a haustrum segment. + """ + + xAround = [] + d1Around = [] + xHalfSetInterHaustra = [] + d1HalfSetInterHaustra = [] + + # Set up profile + if tcCount < 3: # Circular profile + xLoop, d1Loop = createCirclePoints([ 0.0, 0.0, 0.0 ], [ radius, 0.0, 0.0 ], [ 0.0, radius, 0.0 ], + sampleElementOut, startRadians = 0.0) + + else: # tcCount == 3, Triangular profile + cornerRC = cornerInnerRadiusFactor*radius + radiansRangeRC = [7*math.pi/4, 0.0, math.pi/4] + + for n1 in range(3): + radiansAround = n1*2.0*math.pi / 3.0 + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + xc = [(radius - cornerRC) * cosRadiansAround, (radius - cornerRC) * sinRadiansAround, 0.0] + + for n in range(3): + radiansRC = radiansAround + radiansRangeRC[n] + cosRadiansRC = math.cos(radiansRC) + sinRadiansRC = math.sin(radiansRC) + x = [xc[0] + cornerRC*cosRadiansRC, xc[1] + cornerRC*sinRadiansRC, 0.0] + xAround.append(x) + d1 = [ cornerRC*math.pi/4.0 * -sinRadiansRC, cornerRC*math.pi/4.0 * cosRadiansRC, 0.0] + d1Around.append(d1) + + xSample = xAround[1:9] +[xAround[0], xAround[1]] + d1Sample = d1Around[1:9] +[d1Around[0], d1Around[1]] + sx, sd1, _, _, _= interp.sampleCubicHermiteCurves(xSample, d1Sample, sampleElementOut) + xLoop = sx[:-1] + d1Loop = interp.smoothCubicHermiteDerivativesLoop(sx[:-1], sd1[:-1]) + + # Calculate arc length + arcLength = 0.0 + arcDistance = interp.getCubicHermiteArcLength(xLoop[0], d1Loop[0], xLoop[1], d1Loop[1]) + arcLength = arcDistance * sampleElementOut + arcStart = 0.0 + numberOfHalves = tcCount * 2 + arcEnd = arcLength/numberOfHalves + + # Find edge of TC + arcDistanceTCEdge = findEdgeOfTeniaColi(xLoop, d1Loop, tcWidth, arcStart, arcEnd) + + # Sample TC into equally sized elements + xTC, d1TC, _ = sampleTeniaColi(xLoop, d1Loop, arcDistanceTCEdge, elementsCountAroundTC) + + # Sample haustrum into equally sized elements + xHaustrum, d1Haustrum, _, _ = sampleHaustrum(xLoop, d1Loop, xTC[-1], d1TC[-1], + arcLength/numberOfHalves, arcDistanceTCEdge, elementsCountAroundHaustrum) + + xHalfSetInterHaustra = xHalfSetInterHaustra + xTC + xHaustrum[1:] + d1HalfSetInterHaustra = d1HalfSetInterHaustra + d1TC + d1Haustrum[1:] + + return xHalfSetInterHaustra, d1HalfSetInterHaustra + +def createHalfSetIntraHaustralSegment(elementsCountAroundTC, elementsCountAroundHaustrum, + tcCount, tcWidth, radius, cornerInnerRadiusFactor, sampleElementOut, haustrumInnerRadiusFactor): + """ + Find locations and derivative of nodes in half of an intra-haustral + segment. Bow-tie profile for segment with two tenia coli and + clover profile for segment with three tenia coli. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param tcCount: Number of tenia coli. + :param tcWidth: Width of tenia coli. + :param radius: Inner radius of circular inter-haustral profile with two + tenia coli, radius of circle enclosing triangle for inter-haustral profile + with three tenia coli. + :param cornerInnerRadiusFactor: Roundness of triangular corners of + 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 sampleElementOut: Number of sample points used to set up profile + :param haustrumInnerRadiusFactor: Factor is multiplied by inner + radius to obtain radius of intersecting circles in the middle cross-section + along a haustra segment. + :return: Node location and derivative on half of a haustrum segment. + """ + + nxHaustrum = [] + nd1Haustrum = [] + xHalfSetIntraHaustra = [] + d1HalfSetIntraHaustra = [] + + # Set up profile + cornerRC = cornerInnerRadiusFactor*radius + haustrumRadius = (haustrumInnerRadiusFactor + 1)*radius + if tcCount == 2: # Bow-tie profile + originRC = (radius*radius - haustrumRadius*haustrumRadius)/(-2.0*haustrumRadius) + RC = haustrumRadius - originRC + thetaStart = math.asin(-originRC/ RC) + thetaEnd = math.pi - thetaStart + rotOriginRC = [ 0.0, -originRC, 0.0 ] + + else: # tcCount == 3, Clover profile + xc = [(radius - cornerRC)* math.cos(0.0), (radius - cornerRC)*math.sin(0.0), 0.0] + pt1 = [xc[0] + cornerRC*math.cos(0.0), xc[1] + cornerRC*math.sin(0.0), 0.0] + xTC2 = radius* math.cos(2.0*math.pi/3.0) + yTC2 = radius* math.sin(2.0*math.pi/3.0) + originRC = (xTC2*xTC2 + yTC2*yTC2 - haustrumRadius*haustrumRadius) / (2*(-xTC2 - haustrumRadius)) + RC = haustrumRadius - originRC + + # Rotate to find originRC of 1st haustrum + yTC1 = pt1[1] + rotOriginRC = [ originRC*math.cos(-2.0/3.0*math.pi), originRC*math.sin(-2.0/3.0*math.pi), 0.0] + + thetaStart = math.asin((yTC1 + rotOriginRC[1]) / RC) + thetaEnd = math.pi - math.asin((yTC2 + rotOriginRC[1])/ RC) + + thetaHalfHaustrum = (thetaEnd + thetaStart)*0.5 + concaveFactor = 0.15 + thetaConcave = (thetaEnd - thetaStart)*concaveFactor + thetaStart + thetaCircularStart = (thetaEnd - thetaStart)*(concaveFactor + 0.5)*0.5 + thetaStart + + thetaSet = [thetaStart, thetaConcave] + xConcave, _ = getCircleXandD1FromRadians(thetaSet, RC, rotOriginRC) + d1 = [0.0, xConcave[1][1], 0.0] + nxHaustrum = nxHaustrum + xConcave + nd1Haustrum = nd1Haustrum + [d1, d1] + thetaCircularSet = [thetaCircularStart, thetaHalfHaustrum] + xCircular, d1Circular = getCircleXandD1FromRadians(thetaCircularSet, RC, rotOriginRC) + nxHaustrum = nxHaustrum + xCircular + nd1Haustrum = nd1Haustrum + d1Circular + smoothd1 = interp.smoothCubicHermiteDerivativesLine(nxHaustrum, nd1Haustrum, fixStartDirection = True, fixEndDirection = True) + + # Find max x along path + sxHaustrum, sd1Haustrum, _, _ , _ = interp.sampleCubicHermiteCurves(nxHaustrum, smoothd1, sampleElementOut) + xVal = [] + for i in range(len(sxHaustrum)): + xVal.append(sxHaustrum[i][0]) + maxValue = max(xVal) + maxIndex = xVal.index(maxValue) + yAtMaxValue = sxHaustrum[maxIndex][1] + + arcLength = 0.0 + for e in range(len(nxHaustrum)-1): + arcDistance = interp.getCubicHermiteArcLength(nxHaustrum[e], smoothd1[e], nxHaustrum[e+1], smoothd1[e+1]) + arcLength = arcLength + arcDistance + + arcDistanceAtMaxValue = arcLength / sampleElementOut * maxIndex + arcStart = 0.0 if yAtMaxValue > tcWidth*0.5 else arcDistanceAtMaxValue + arcEnd = arcDistanceAtMaxValue if yAtMaxValue > tcWidth*0.5 else arcLength + + # Find edge of TC + arcDistanceTCEdge = findEdgeOfTeniaColi(nxHaustrum, smoothd1, tcWidth, arcStart, arcEnd) + + # Sample TC into equally sized elements + xTC, d1TC, arcLengthPerTC = sampleTeniaColi(nxHaustrum, smoothd1, arcDistanceTCEdge, elementsCountAroundTC) + + # Sample haustrum into equally sized elements + xHaustrum, d1Haustrum, arcLengthPerHaustrum, arcLengthPerTransition = sampleHaustrum(nxHaustrum, smoothd1, xTC[-1], d1TC[-1], arcLength, arcDistanceTCEdge, elementsCountAroundHaustrum) + + xHalfSetIntraHaustra = xHalfSetIntraHaustra + xTC + xHaustrum[1:] + d1HalfSetIntraHaustra = d1HalfSetIntraHaustra + d1TC + d1Haustrum[1:] + + return xHalfSetIntraHaustra, d1HalfSetIntraHaustra + +def findEdgeOfTeniaColi(nx, nd1, tcWidth, arcStart, arcEnd): + """ + Locate edge of tenia coli on a cubic hermite interpolated + curve defined by nodes nx and derivatives nd1. + :param nx: Coordinates of nodes along curve. + :param nd1: Derivatives of nodes along curve. + :param tcWidth: Width of tenia coli. + :param arcStart: Lower limit of arc distance to search for + edge of tenia coli. + :param arcEnd: Upper limit of arc distance to search for + edge of tenia coli. + :return: arc distance covered by tenia coli. + """ + xTol = 1.0E-6 + for iter in range(100): + arcDistance = (arcStart + arcEnd)*0.5 + x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) + diff = x[1] - tcWidth*0.5 + if abs(diff) > xTol: + if diff < 0.0: + arcStart = arcDistance + else: + arcEnd = arcDistance + else: + arcDistanceTCEdge = arcDistance + break + if iter > 99: + print('Search for TC boundary - Max iters reached:',iter) + + return arcDistanceTCEdge + +def sampleTeniaColi(nx, nd1, arcDistanceTCEdge, elementsCountAroundTC): + """ + Get systematically spaced points and derivatives over the width of + tenia coli over a cubic hermite interpolated curves defined by nodes + nx and derivatives nd1 lying along half of a haustrum. + :param nx: Coordinates of nodes along curve. + :param nd1: Derivatives of nodes along curve. + :param arcDistanceTCEdge: Arc distance covered by tenia coli. + :param elementsCountAroundTC: Number of elements around tenia coli. + :return: coordinates, derivatives of tenia coli, and arclength of + each tenia coli element. + """ + xTC = [] + d1TC = [] + arcDistancePerElementTC = arcDistanceTCEdge / (elementsCountAroundTC*0.5) + for e in range(int(elementsCountAroundTC*0.5)+1): + arcDistance = arcDistancePerElementTC * e + x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) + d1Scaled = vector.setMagnitude(d1, arcDistancePerElementTC) + xTC.append(x) + d1TC.append(d1Scaled) + + return xTC, d1TC, arcDistancePerElementTC + +def sampleHaustrum(nx, nd1, xTCLast, d1TCLast, arcLength, arcDistanceTCEdge, + elementsCountAroundHaustrum): + """ + Get systematically spaced points and derivatives over the width of + haustrum over a cubic hermite interpolated curves defined by nodes + nx and derivatives nd1 lying along half of a haustrum. + :param nx: Coordinates of nodes along curve. + :param nd1: Derivatives of nodes along curve. + :param xTCLast: Coordinates of node lying on edge of tenia coli. + :param d1TCLast: Derivative of node lying on edge of tenia coli. + :param arcLength: Arc length of curve. + :param arcDistanceTCEdge: Arc distance covered by tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :return: coordinates, derivatives on haustrum. + :return: arclength of haustrum element and transition elements. + """ + xHaustrum = [] + d1Haustrum = [] + elementLengths = [] + length = arcLength - arcDistanceTCEdge + elementsCountOut = int(elementsCountAroundHaustrum * 0.5) + addLengthStart = 0.5 * vector.magnitude(d1TCLast) + lengthFractionStart = 0.5 + proportionStart = 1.0 + elementLengthMid = (length - addLengthStart) / (elementsCountOut - 1.0 + proportionStart*lengthFractionStart ) + elementLengthProportionStart = proportionStart*lengthFractionStart*elementLengthMid + elementLengthProportionEnd = elementLengthMid + + for e in range(elementsCountOut): + xi = e/(elementsCountOut - 1) + elementLengths.append(elementLengthMid) + elementLengths[0] = addLengthStart + elementLengthProportionStart + elementLengths[-1] = 0.0 + elementLengthProportionEnd + haustrumElementLength = elementLengthMid + transitionElementLength = elementLengths[0] + + arcDistance = arcDistanceTCEdge + xHaustrum.append(xTCLast) + d1Scaled = vector.setMagnitude(d1TCLast, elementLengths[0]) + d1Haustrum.append(d1Scaled) + + for e in range(elementsCountOut): + arcDistance = arcDistance + elementLengths[e] + x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) + d1Scaled = vector.setMagnitude(d1, elementLengths[e] if e > 0 else elementLengths[e+1]) + xHaustrum.append(x) + d1Haustrum.append(d1Scaled) + + return xHaustrum, d1Haustrum, haustrumElementLength, transitionElementLength + +def getCircleXandD1FromRadians(thetaSet, radius, origin): + """ + Gets the coordinates and derivatives along a circular path + based on an angular range. + :param thetaSet: Lower and upper limit of theta. + :param radius: Radius of circle. + :param origin: Origin of circle. + :return: coordinates, derivatives on lower and upper + limit of thetaSet. + """ + nx = [] + nd1 = [] + dTheta = thetaSet[1] - thetaSet[0] + for n in range(2): + theta = thetaSet[n] + x = [radius*math.cos(theta) - origin[0], + radius*math.sin(theta) - origin[1], + 0.0] + d1 = [-radius*math.sin(theta)*dTheta, + radius*math.cos(theta)*dTheta, + 0.0] + nx.append(x) + nd1.append(d1) + + return nx, nd1 + +def getFullProfileFromHalfHaustrum(xHaustrumHalfSet, d1HaustrumHalfSet, + d2HaustrumHalfSet, tcCount): + """ + Gets the coordinates and derivatives of the entire profile + using points from first half of the first sector. The first + sector starts from the x-axis. The first half set of points + are reflected across the x-axis followed by rotation to get + the points in the second half of the first sector. The full + set of points in the first sector are then rotated to obtain + points in the other two sectors. + :param xHaustrumHalfSet: Coordinates of points in first + half of the first sector. + :param d1HaustrumHalfSet: Derivatives of points around first + half of the first sector. + :param d2HaustrumHalfSet: Derivatives of points along first + half of the first sector. + :param tcCount: Number of tenia coli. + :return: coordinates, derivatives of points over entire profile. + """ + xHaustrumHalfSet2 = [] + d1HaustrumHalfSet2 = [] + d2HaustrumHalfSet2 = [] + xHaustrum = [] + d1Haustrum = [] + d2Haustrum = [] + xHaustra = [] + d1Haustra = [] + d2Haustra = [] + rotAng = 2*math.pi/tcCount + + for n in range(1,len(xHaustrumHalfSet)): + idx = -n + len(xHaustrumHalfSet) - 1 + x = xHaustrumHalfSet[idx] + d1 = d1HaustrumHalfSet[idx] + xReflect = [x[0], -x[1], x[2]] + d1Reflect = [d1[0], -d1[1], d1[2]] + xRot = [xReflect[0]*math.cos(rotAng) - xReflect[1]*math.sin(rotAng), + xReflect[0]*math.sin(rotAng) + xReflect[1]*math.cos(rotAng), + xReflect[2]] + d1Rot = [-(d1Reflect[0]*math.cos(rotAng) - d1Reflect[1]*math.sin(rotAng)), + -(d1Reflect[0]*math.sin(rotAng) + d1Reflect[1]*math.cos(rotAng)), + -d1Reflect[2]] + d2 = d2HaustrumHalfSet[idx] + d2Reflect = [d2[0], -d2[1], d2[2]] + d2Rot = [(d2Reflect[0]*math.cos(rotAng) - d2Reflect[1]*math.sin(rotAng)), + (d2Reflect[0]*math.sin(rotAng) + d2Reflect[1]*math.cos(rotAng)), + d2Reflect[2]] + xHaustrumHalfSet2.append(xRot) + d1HaustrumHalfSet2.append(d1Rot) + d2HaustrumHalfSet2.append(d2Rot) + + xHaustrum = xHaustrumHalfSet + xHaustrumHalfSet2 + d1Haustrum = d1HaustrumHalfSet + d1HaustrumHalfSet2 + d2Haustrum = d2HaustrumHalfSet + d2HaustrumHalfSet2 + + # Rotate to get all 3 sectors + xHaustra = xHaustra + xHaustrum[:-1] + d1Haustra = d1Haustra + d1Haustrum[:-1] + d2Haustra = d2Haustra + d2Haustrum[:-1] + + ang = [ 2/3*math.pi, -2/3*math.pi] if tcCount == 3 else [math.pi] + for i in range(tcCount - 1): + rotAng = ang[i] + cosRotAng = math.cos(rotAng) + sinRotAng = math.sin(rotAng) + for n in range(len(xHaustrum)- 1): + x = xHaustrum[n] + d1 = d1Haustrum[n] + d2 = d2Haustrum[n] + x = [ x[0]*cosRotAng - x[1]*sinRotAng, x[0]*sinRotAng + x[1]*cosRotAng, x[2]] + xHaustra.append(x) + dx_ds1 = [ d1[0]*cosRotAng - d1[1]*sinRotAng, d1[0]*sinRotAng + d1[1]*cosRotAng, d1[2]] + d1Haustra.append(dx_ds1) + dx_ds2 = [ d2[0]*cosRotAng - d2[1]*sinRotAng, d2[0]*sinRotAng + d2[1]*cosRotAng, d2[2]] + d2Haustra.append(dx_ds2) + + return xHaustra, d1Haustra, d2Haustra + +def getXiListFromOuterLengthProfile(xInner, d1Inner, segmentAxis, + wallThickness, transitElementList): + """ + Gets a list of xi for flat and texture coordinates calculated + from outer arclength of elements around a segment (most relaxed state). + :param xInner: Coordinates of points on inner surface around segment. + :param d1Inner: Derivatives of points on inner surface around segment. + :param segmentAxis: Axis of segment. + :param wallThickness: Thickness of wall. + :param transitElementList: stores true if element around is an element that + transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. + :return xiList: List containing xi for each point on outer surface extruded by + wall thickness from inner points. + :return totalArcLengthOuter: Total arclength around outer surface of elements. + """ + unitNormList = [] + xOuter = [] + curvatureInner = [] + d1Outer = [] + + for n in range(len(xInner)): + unitNormList.append(vector.normalise(vector.crossproduct3(d1Inner[n], segmentAxis))) + + for n in range(len(xInner)): + norm = unitNormList[n] + # Calculate outer coordinates + x = [xInner[n][i] + norm[i]*wallThickness for i in range(3)] + xOuter.append(x) + # Calculate curvature along elements around + prevIdx = n - 1 if (n != 0) else len(xInner) - 1 + nextIdx = n + 1 if (n < (len(xInner) - 1)) else 0 + kappam = interp.getCubicHermiteCurvatureSimple(xInner[prevIdx], d1Inner[prevIdx], xInner[n], d1Inner[n], 1.0) + kappap = interp.getCubicHermiteCurvatureSimple(xInner[n], d1Inner[n], xInner[nextIdx], d1Inner[nextIdx], 0.0) + if not transitElementList[n] and not transitElementList[(n-1)%(len(xInner))]: + curvatureAround = 0.5*(kappam + kappap) + elif transitElementList[n]: + curvatureAround = kappam + elif transitElementList[(n-1)%(len(xInner))]: + curvatureAround = kappap + curvatureInner.append(curvatureAround) + + for n in range(len(xOuter)): + factor = 1.0 + wallThickness * curvatureInner[n] + dx_ds1 = [ factor*c for c in d1Inner[n]] + d1Outer.append(dx_ds1) + + arcLengthList = [] + for n1 in range(len(xOuter)): + arcLengthPerElement = interp.getCubicHermiteArcLength(xOuter[n1], d1Outer[n1], + xOuter[(n1+1) % len(xOuter)], d1Outer[(n1+1)% len(xOuter)]) + arcLengthList.append(arcLengthPerElement) + + # Total arcLength + totalArcLengthOuter = 0.0 + for n1 in range(len(arcLengthList)): + totalArcLengthOuter += arcLengthList[n1] + + xiList = [0.0] + arcDistance = 0 + for n in range(len(arcLengthList)): + arcDistance = arcDistance + arcLengthList[n] + xi = arcDistance / totalArcLengthOuter + xiList.append(xi) + + return xiList, totalArcLengthOuter + +def getTeniaColi(region, xList, d1List, d2List, d3List, curvatureList, + tcCount, elementsCountAroundTC, elementsCountAroundHaustrum, + elementsCountAlong, elementsCountThroughWall, + tubeTCWidthList, tcThickness, sx, + annotationGroups, annotationArray): + """ + 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 + boundary on the outer surface of the haustra, a midpoint + lying at tenia coli thickness normal to the midpoint tenia coli + boundary on the outer surface, and the right edge of tenia coli + boundary. Function incorporates annotation groups for tenia coli, + and arranges coordinates, derivatives of tenia coli to follow + after points on outer surface of the haustra as the mesh extends + along its length. + :param xList, d1List, d2List, d3List: Coordinates and derivatives of nodes + on haustra. + :param curvatureList: Curvature of points along length. + :param tcCount: Number of tenia coli. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param elementsCountAlong: Number of elements along colon. + :param elementsCountThroughWall: Number of elements through wall. + :param 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 + :return: coordinates, derivatives, annotationGroups for colon with tenia + coli. + """ + + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount + TCEdgeFactor = 1.5 + xTCArranged = [] + d1TCArranged = [] + d2TCArranged = [] + d3TCArranged = [] + + # Create points for tenia coli + for n2 in range(elementsCountAlong + 1): + xTCRaw = [] + d1TCRaw = [] + d2TCRaw = [] + d3TCRaw = [] + tcWidth = tubeTCWidthList[n2] + for N in range(tcCount): + idxTCMid = 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) + v1 = xList[TCStartIdx] + v2 = xMid + d1MidScaled = [c*tcWidth*TCEdgeFactor for c in vector.normalise(d1Mid)] + v3 = xList[TCEndIdx] + nx = [v1, v2, v3] + nd1 = [d1List[TCStartIdx], d1MidScaled, d1List[TCEndIdx]] + sxTC, sd1TC, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) + xTCRaw = xTCRaw + sxTC[1:-1] + if elementsCountAroundTC == 2: + p = [v2[i] - v1[i] for i in range(3)] + A = vector.dotproduct(unitNorm, p) # A<0 if v2 is higher than v1 + d1 = [c*tcWidth*0.5 for c in vector.normalise(d1Mid)] if A < 0 else d1MidScaled + d1TCRaw = d1TCRaw + sd1TC[1:-1] if elementsCountAroundTC > 2 else d1TCRaw + [d1] + xTCInnerSet = list(range(TCStartIdx+1, TCEndIdx)) if N > 0 else list(range(TCStartIdx + 1, TCStartIdx + int(elementsCountAroundTC * 0.5))) + list(range(idxTCMid, idxTCMid + int(elementsCountAroundTC * 0.5))) + + for n in range(elementsCountAroundTC - 1): + xTCInner = xList[xTCInnerSet[n]] + xTCOuter = sxTC[n + 1] + d3 = [xTCOuter[i] - xTCInner[i] for i in range(3)] + d3TCRaw.append(d3) + d3List[xTCInnerSet[n]] = d3 + + innerIdx = xTCInnerSet[n] - elementsCountThroughWall*elementsCountAround + curvature = curvatureList[innerIdx] + distanceToInnerIdx = vector.magnitude([xTCOuter[i] - xList[innerIdx][i] for i in range(3)]) + factor = 1.0 - curvature*distanceToInnerIdx + d2 = [factor*c for c in d2List[innerIdx]] + d2TCRaw.append(d2) + + xTCArranged = xTCArranged + xTCRaw[int(elementsCountAroundTC*0.5 -1):] + xTCRaw[:int(elementsCountAroundTC*0.5 -1)] + d1TCArranged = d1TCArranged + d1TCRaw[int(elementsCountAroundTC*0.5 -1):] + d1TCRaw[:int(elementsCountAroundTC*0.5 -1)] + d2TCArranged = d2TCArranged + d2TCRaw[int(elementsCountAroundTC*0.5 -1):] + d2TCRaw[:int(elementsCountAroundTC*0.5 -1)] + d3TCArranged = d3TCArranged + d3TCRaw[int(elementsCountAroundTC*0.5 -1):] + d3TCRaw[:int(elementsCountAroundTC*0.5 -1)] + + x, d1, d2, d3 = combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTCArranged, d1TCArranged, + d2TCArranged, d3TCArranged, (elementsCountAroundTC - 1)*tcCount, elementsCountAround, + elementsCountAlong, elementsCountThroughWall) + + # Update annotation groups + if tcCount == 3: + tlGroup = AnnotationGroup(region, 'tenia libera', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + tmGroup = AnnotationGroup(region, 'tenia mesocolica', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + toGroup = AnnotationGroup(region, 'tenia omentalis', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') + annotationGroupsTC = [tlGroup, tmGroup, toGroup] + annotationArrayTC = (['tenia omentalis']*int(elementsCountAroundTC*0.5) + + ['tenia libera']*elementsCountAroundTC + + ['tenia mesocolica']*elementsCountAroundTC + + ['tenia omentalis']*int(elementsCountAroundTC*0.5)) + annotationGroups += annotationGroupsTC + annotationArray += annotationArrayTC + + return x, d1, d2, d3, annotationGroups, annotationArray + +def combineTeniaColiWithColon(xList, d1List, d2List, d3List, xTC, d1TC, d2TC, + d3TC, nodesCountAroundTC, elementsCountAround, elementsCountAlong, + elementsCountThroughWall): + """ + Arranges coordinates and derivatives around inner surface to + outer surface, followed by tenia coli points before extending + along length of colon. + :param xList, d1List, d2List, d3List: coordinates and derivatives of colon. + :param xTC, d1TC, d2TC, d3TC: coordinates and derivatives of tenia coli. + :param nodesCountAroundTC: Number of nodes around tenia coli. + :param elementsCountAround: Number of elements around colon. + :param elementsCountAlong: Number of elements along colon. + :param elementsCountThroughWall: Number of elements through wall. + : return: reordered coordinates and derivatives + """ + x = [] + d1 = [] + d2 = [] + d3 = [] + + # 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 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]) + + return x, d1, d2, d3 + +def createFlatAndTextureCoordinatesTeniaColi(xiList, relaxedLengthList, + totalLengthAlong, wallThickness, tcCount, tcThickness, + elementsCountAroundTC, elementsCountAroundHaustrum, + elementsCountAlong, elementsCountThroughWall, transitElementList): + """ + Calculates flat coordinates and texture coordinates + for a colon scaffold with tenia coli when it is opened + up into a flat preparation. + :param xiList: List containing xi for each point around the outer surface of + colon in its most relaxed state. + :param relaxedLengthList: List of total arclength around the outer surface in + its most relaxed state for each element along. + :param totalLengthAlong: Total length along colon. + :param wallThickness: Thickness of wall. + :param tcCount: Number of tenia coli. + :param tcThickness: Thickness of tenia coli at its thickest region. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param elementsCountAlong: Number of elements along colon. + :param elementsCountThroughWall: Number of elements through wall. + :param transitElementList: stores true if element around is an element that + transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. + :return: coordinates and derivatives of flat and texture coordinates fields. + """ + + # Calculate flat coordinates + factor = 3.0 if tcCount == 3 else 2.0 + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount + + # Find flat and texture coordinates for colon + xFlatColon, d1FlatColon, d2FlatColon, xTextureColon, d1TextureColon, d2TextureColon \ + = tubemesh.createFlatAndTextureCoordinates(xiList, relaxedLengthList, + totalLengthAlong, wallThickness, elementsCountAround, elementsCountAlong, + elementsCountThroughWall, transitElementList) + + # Find flat coordinates for tenia coli + xFlatListTC = [] + d1FlatListTC = [] + d2FlatListTC = [] + + for n2 in range(elementsCountAlong + 1): + xiFace = xiList[n2] + relaxedLength = relaxedLengthList[n2] + xPad = (relaxedLengthList[0] - relaxedLength)*0.5 + for N in range(tcCount + 1): + idxTCMid = N*(elementsCountAroundTC + elementsCountAroundHaustrum) + TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) + TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) + dTC = (xiFace[idxTCMid] - xiFace[TCStartIdx])*relaxedLength if N > 0 else (xiFace[TCEndIdx] - xiFace[idxTCMid])*relaxedLength + v1 = [xiFace[TCStartIdx]*relaxedLength, 0.0, wallThickness] if N > 0 else [-dTC, 0.0, wallThickness] + v2 = [ xiFace[idxTCMid]*relaxedLength, 0.0, wallThickness + tcThickness] + v3 = [ xiFace[TCEndIdx]*relaxedLength, 0.0, wallThickness] if N < tcCount else [ relaxedLength + dTC, 0.0, wallThickness] + d1 = d3 = [dTC, 0.0, 0.0] + d2 = [c*factor for c in [dTC, 0.0, 0.0]] + nx = [v1, v2, v3] + nd1 = [d1, d2, d3] + sx, sd1, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) + if N > 0 and N < tcCount: + w = sx[1:-1] + dw = sd1[1:-1] if elementsCountAroundTC > 2 else nd1[1:-1] + elif N == 0: + w = sx[int(elementsCountAroundTC*0.5):-1] + dw = sd1[int(elementsCountAroundTC*0.5):-1] if elementsCountAroundTC > 2 else nd1[int(elementsCountAroundTC*0.5):-1] + else: + w = sx[1:int(elementsCountAroundTC*0.5)+1] + dw = sd1[1:int(elementsCountAroundTC*0.5)+1] if elementsCountAroundTC > 2 else nd1[1:int(elementsCountAroundTC*0.5)+1] + for n in range(len(w)): + x = [ xPad + w[n][0], + totalLengthAlong / elementsCountAlong * n2, + w[n][2]] + xFlatListTC.append(x) + d1FlatListTC.append(dw[n]) + + for n2 in range(elementsCountAlong): + for n1 in range((elementsCountAroundTC - 1)*tcCount + 1): + nodeIdx = n2*((elementsCountAroundTC - 1)*tcCount + 1) + n1 + nodeNextElementAlong = nodeIdx + ((elementsCountAroundTC - 1)*tcCount + 1) + v1 = xFlatListTC[nodeNextElementAlong] + v2 = xFlatListTC[nodeIdx] + d1 = d2 = [v1[i] - v2[i] for i in range(3)] + arclength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) + d2Flat = vector.setMagnitude(d1, arclength) + d2FlatListTC.append(d2Flat) + d2FlatListTC = d2FlatListTC + d2FlatListTC[-((elementsCountAroundTC - 1)*tcCount + 1):] + + # Find texture coordinates for tenia coli + wList = [] + dwList = [] + xiTexture = xiList[0] + xTextureListTC = [] + d1TextureListTC = [] + d2TextureListTC = [] + + for N in range(tcCount + 1): + idxTCMid = N*(elementsCountAroundTC + elementsCountAroundHaustrum) + TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) + TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) + dTC = xiTexture[idxTCMid] - xiTexture[TCStartIdx] if N > 0 else xiTexture[TCEndIdx] - xiTexture[idxTCMid] + v1 = [xiTexture[TCStartIdx], 0.0, 1.0] if N > 0 else [-dTC, 0.0, 1.0] + v2 = [ xiTexture[idxTCMid], 0.0, 2.0] + v3 = [ xiTexture[TCEndIdx], 0.0, 1.0] if N < tcCount else [ 1 + dTC , 0.0, 1.0] + d1 = d2 = d3 = [dTC, 0.0, 0.0] + nx = [v1, v2, v3] + nd1 = [d1, d2, d3] + sx, sd1, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) + if N > 0 and N < tcCount: + w = sx[1:-1] + dw = sd1[1:-1] if elementsCountAroundTC > 2 else nd1[1:-1] + elif N == 0: + w = sx[int(elementsCountAroundTC*0.5):-1] + dw = sd1[int(elementsCountAroundTC*0.5):-1] if elementsCountAroundTC > 2 else nd1[int(elementsCountAroundTC*0.5):-1] + else: + w = sx[1:int(elementsCountAroundTC*0.5)+1] + dw = sd1[1:int(elementsCountAroundTC*0.5)+1] if elementsCountAroundTC > 2 else nd1[1:int(elementsCountAroundTC*0.5)+1] + wList = wList + w + dwList = dwList + dw + + d2 = [0.0, 1.0 / elementsCountAlong, 0.0] + for n2 in range(elementsCountAlong + 1): + for n1 in range((elementsCountAroundTC-1) * tcCount + 1): + u = [ wList[n1][0], + 1.0 / elementsCountAlong * n2, + wList[n1][2]] + xTextureListTC.append(u) + d1TextureListTC.append(dwList[n1]) + d2TextureListTC.append(d2) + + xFlat, d1Flat, d2Flat, _ = combineTeniaColiWithColon(xFlatColon, d1FlatColon, d2FlatColon, [], + xFlatListTC, d1FlatListTC, d2FlatListTC, [], (elementsCountAroundTC - 1)*tcCount + 1, + elementsCountAround + 1, elementsCountAlong, elementsCountThroughWall) + + xTexture, d1Texture, d2Texture, _ = combineTeniaColiWithColon(xTextureColon, d1TextureColon, + d2TextureColon, [], xTextureListTC, d1TextureListTC, d2TextureListTC, [], + (elementsCountAroundTC - 1)*tcCount + 1, elementsCountAround + 1, elementsCountAlong, + elementsCountThroughWall) + + return xFlat, d1Flat, d2Flat, xTexture, d1Texture, d2Texture + +def createNodesAndElementsTeniaColi(region, + x, d1, d2, d3, + xFlat, d1Flat, d2Flat, + xTexture, d1Texture, d2Texture, + elementsCountAroundTC, elementsCountAroundHaustrum, + elementsCountAlong, elementsCountThroughWall, tcCount, + annotationGroups, annotationArray, + firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives): + """ + Create nodes and elements for the coordinates, flat coordinates, + and texture coordinates field. + :param x, d1, d2, d3: coordinates and derivatives of coordinates field. + :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of + flat coordinates field. + :param xTexture, d1Texture, d2Texture, d3Texture: coordinates and derivatives + of texture coordinates field. + :param elementsCountAroundTC: Number of elements around tenia coli. + :param elementsCountAroundHaustrum: Number of elements around haustrum. + :param elementsCountAlong: Number of elements along colon. + :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 firstNodeIdentifier, firstElementIdentifier: first node and element + identifier to use. + :param useCubicHermiteThroughWall: use linear when false + :param useCrossDerivatives: use cross derivatives when true + :return nodeIdentifier, elementIdentifier, annotationGroups + """ + + nodeIdentifier = firstNodeIdentifier + elementIdentifier = firstNodeIdentifier + elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum )*tcCount + + # Create coordinates field + zero = [ 0.0, 0.0, 0.0 ] + fm = region.getFieldmodule() + fm.beginChange() + cache = fm.createFieldcache() + coordinates = zinc_utils.getOrCreateCoordinateField(fm) + + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + if useCubicHermiteThroughWall: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) + + mesh = fm.findMeshByDimension(3) + + if useCubicHermiteThroughWall: + eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives) + else: + eftfactory = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) + eft = eftfactory.createEftBasic() + + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + result = elementtemplate.defineField(coordinates, -1, eft) + + # Tenia coli edge elements + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + eft1 = eftfactory.createEftWedgeXi1One() + elementtemplate1.defineField(coordinates, -1, eft1) + + elementtemplate2 = mesh.createElementtemplate() + elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + eft2 = eftfactory.createEftWedgeXi1Zero() + elementtemplate2.defineField(coordinates, -1, eft2) + + # Create flat coordinates field + flatCoordinates = zinc_utils.getOrCreateFlatCoordinateField(fm) + 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() + eftTexture5 = bicubichermitelinear.createEftWedgeXi1One() + 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 = zinc_utils.getOrCreateTextureCoordinateField(fm) + textureNodetemplate1 = nodes.createNodetemplate() + textureNodetemplate1.defineField(textureCoordinates) + textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) + textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + + textureNodetemplate2 = nodes.createNodetemplate() + textureNodetemplate2.defineField(textureCoordinates) + textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) + textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) + textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) + if useCrossDerivatives: + textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) + + textureElementtemplate1 = mesh.createElementtemplate() + textureElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate1.defineField(textureCoordinates, -1, eftTexture3) + + textureElementtemplate2 = mesh.createElementtemplate() + textureElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate2.defineField(textureCoordinates, -1, eftTexture4) + + textureElementtemplate3 = mesh.createElementtemplate() + textureElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate3.defineField(textureCoordinates, -1, eftTexture5) + + textureElementtemplate4 = mesh.createElementtemplate() + textureElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate4.defineField(textureCoordinates, -1, eftTexture6) + + textureElementtemplate5 = mesh.createElementtemplate() + textureElementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) + textureElementtemplate5.defineField(textureCoordinates, -1, eftTexture7) + + # create nodes for coordinates field + for n in range(len(x)): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, d3[n]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + # print('NodeIdentifier = ', nodeIdentifier, xList[n]) + 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) + 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 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 + 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] + 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) + elementIdentifier = elementIdentifier + 1 + if annotationGroups: + for annotationGroup in annotationGroups: + if annotationArray[e1] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + # 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: + nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] + else: + nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, 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) + elementIdentifier = elementIdentifier + 1 + if tcCount == 3: + for annotationGroup in annotationGroups: + if annotationArray[elementsCountAround + eTC] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) + + 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 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) + 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) + 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) + 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 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 eTC > 0: + nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni31, bni32, bni31 + now + tcOffset, bni32 + now + tcOffset] + else: + nodeIdentifiers = [bni21, bni22, bni21 + now + tcOffset, bni22 + now + tcOffset, bni32, bni32 + now + tcOffset] + 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) + 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) + + fm.endChange() + + return nodeIdentifier, elementIdentifier, annotationGroups diff --git a/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentsimplemesentery1.py b/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentsimplemesentery1.py deleted file mode 100644 index 44de0f52..00000000 --- a/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentsimplemesentery1.py +++ /dev/null @@ -1,258 +0,0 @@ -""" -Generates a single 3-D colon segment with a simple -mesentery mesh along a central line, with variable -numbers of elements around, along and through wall, -with variable radius and thickness along. -""" - -import math -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup -from scaffoldmaker.meshtypes.meshtype_3d_colonsegmentteniacoli1 import sampleHaustrum, getuListFromOuterMidLengthProfile -from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils.geometry import createCirclePoints -from scaffoldmaker.utils.meshrefinement import MeshRefinement -from scaffoldmaker.utils import tubemesh - -class MeshType_3d_colonsegmentsimplemesentery1(Scaffold_base): - ''' - Generates a single 3-D colon segment with simple mesentery mesh - with variable numbers of elements around, along the central line, - and through wall. The cross-section profile of the colon - segment is largely cylindrical due to the lack of a tenia coli. - ''' - @staticmethod - def getName(): - return '3D Colon Segment Simple Mesentery 1' - - @staticmethod - def getParameterSetNames(): - return [ - 'Default', - 'Mouse 1'] - - @staticmethod - def getDefaultOptions(parameterSetName='Default'): - return { - 'Number of elements around mesenteric zone' : 2, - 'Number of elements around non-mesenteric zone' : 10, - 'Number of elements along segment' : 4, - 'Number of elements through wall' : 1, - 'Inner radius': 0.2, - 'Mesenteric zone width': 0.08, - 'Segment length': 1.5, - 'Wall thickness': 0.02, - 'Use cross derivatives' : False, - 'Use linear through wall' : True, - 'Refine' : False, - 'Refine number of elements around' : 1, - 'Refine number of elements along segment' : 1, - 'Refine number of elements through wall' : 1 - } - - @staticmethod - def getOrderedOptionNames(): - return [ - 'Number of elements around mesenteric zone', - 'Number of elements around non-mesenteric zone', - 'Number of elements along segment', - 'Number of elements through wall', - 'Inner radius', - 'Mesenteric zone width', - 'Segment length', - 'Wall thickness', - 'Use cross derivatives', - 'Use linear through wall', - 'Refine', - 'Refine number of elements around', - 'Refine number of elements along segment', - 'Refine number of elements through wall' - ] - - @staticmethod - def checkOptions(options): - for key in [ - 'Number of elements along segment', - 'Number of elements through wall', - 'Refine number of elements around', - 'Refine number of elements along segment', - 'Refine number of elements through wall']: - if options[key] < 1: - options[key] = 1 - if options['Number of elements around mesenteric zone'] < 2: - options['Number of elements around mesenteric zone'] = 2 - if options['Number of elements around non-mesenteric zone'] < 4: - options['Number of elements around non-mesenteric zone'] = 4 - for key in [ - 'Number of elements around mesenteric zone', - 'Number of elements around non-mesenteric zone']: - if options[key] % 2 > 0: - options[key] = options[key] + 1 - for key in [ - 'Inner radius', - 'Mesenteric zone width', - 'Segment length', - 'Wall thickness']: - if options[key] < 0.0: - options[key] = 0.0 - if options['Mesenteric zone width'] < 10.0*math.pi/180.0*options['Inner radius']: - options['Mesenteric zone width'] = round(10.0*math.pi/180.0*options['Inner radius'], 2) - if options['Mesenteric zone width'] > math.pi*0.5*options['Inner radius']: - options['Mesenteric zone width'] = round(math.pi*0.5*options['Inner radius'], 2) - - @staticmethod - def generateBaseMesh(region, options): - """ - Generate the base tricubic Hermite mesh. See also generateMesh(). - :param region: Zinc region to define model in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - :return: None - """ - elementsCountAroundMZ = options['Number of elements around mesenteric zone'] - elementsCountAroundNonMZ = options['Number of elements around non-mesenteric zone'] - elementsCountAround = elementsCountAroundMZ + elementsCountAroundNonMZ - elementsCountAlongSegment = options['Number of elements along segment'] - elementsCountThroughWall = options['Number of elements through wall'] - radius = options['Inner radius'] - segmentLength = options['Segment length'] - mzWidth = options['Mesenteric zone width'] - wallThickness = options['Wall thickness'] - useCrossDerivatives = options['Use cross derivatives'] - useCubicHermiteThroughWall = not(options['Use linear through wall']) - segmentCount = 1 - - cx = [ [ 0.0, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] - cd1 = [ [ segmentLength, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] - cd2 = [ [ 0.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ] ] - cd12 = [ [0.0, 0.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ] - - # Generate inner surface of a colon segment - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength, xInner, d1Inner, d2Inner, segmentAxis = getColonSegmentInnerPointsNoTeniaColi(region, elementsCountAroundMZ, - elementsCountAroundNonMZ, elementsCountAlongSegment, mzWidth, radius, segmentLength, wallThickness) - - # Generate tube mesh - annotationGroups, nextNodeIdentifier, nextElementIdentifier, xList, d1List, d2List, d3List, sx, curvatureAlong, factorList = tubemesh.generatetubemesh(region, - elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, segmentCount, cx, cd1, cd2, cd12, - xInner, d1Inner, d2Inner, wallThickness, segmentAxis, segmentLength, useCrossDerivatives, useCubicHermiteThroughWall, - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength) - - return annotationGroups - - @classmethod - def generateMesh(cls, region, options): - """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - - refineElementsCountAround = options['Refine number of elements around'] - refineElementsCountAlong = options['Refine number of elements along segment'] - refineElementsCountThroughWall = options['Refine number of elements through wall'] - - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) - return meshrefinement.getAnnotationGroups() - -def getColonSegmentInnerPointsNoTeniaColi(region, elementsCountAroundMZ, elementsCountAroundNonMZ, elementsCountAlongSegment, - mzWidth, radius, segmentLength, wallThickness): - """ - Generates a 3-D colon segment mesh with a simple mesentery - (no tenia coli) with variable numbers of elements around, - along the central line, and through wall. The colon segment - has a cylindrical profile. - :param elementsCountAroundMZ: Number of elements around mesenteric zone. - :param elementsCountAroundNonMZ: Number of elements around non-mesenteric zone. - :param elementsCountAlongSegment: Number of elements along colon segment. - :param mzWidth: Width of mesenteric zone in flat preparation. - :param radius: Inner radius of colon segment. - :param segmentLength: Length of a colon segment. - :param wallThickness: Thickness of wall. - :return annotationGroups, annotationArray: annotationArray stores annotation - names of elements around - :return transitElementList: stores true if element around is an element that - transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. - :return uList: List of xi for node around - : return totalArcLengthOuterMidLength: total arclength of elements on outer - surface along mid-length of segment. - :return coordinates, derivatives on inner surface of a colon segment. - """ - mzGroup = AnnotationGroup(region, 'mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - nonmzGroup = AnnotationGroup(region, 'non-mesenteric zone', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - annotationGroups = [mzGroup, nonmzGroup] - annotationArray = ['mesenteric zone']*int(elementsCountAroundMZ*0.5) + ['non-mesenteric zone']*elementsCountAroundNonMZ + ['mesenteric zone']*int(elementsCountAroundMZ*0.5) - transitElementList = [0]*int(elementsCountAroundMZ*0.5) + [1] + [0]*int(elementsCountAroundNonMZ - 2) + [1] + [0]*int(elementsCountAroundMZ*0.5) - - # create nodes - x = [ 0.0, 0.0, 0.0 ] - d1 = [ 0.0, 0.0, 0.0 ] - segmentAxis = [0.0, 0.0, 1.0] - - xMZ = [] - d1MZ = [] - xListBaseLayer = [] - d1ListBaseLayer = [] - xList = [] - d1List = [] - d2List = [] - uList = [] - - #Set up profile - sampleElementOut = 20 - nx, nd1 = createCirclePoints([ 0.0, 0.0, 0.0 ], [ radius, 0.0, 0.0 ], [ 0.0, radius, 0.0 ], sampleElementOut, startRadians = 0.0) - - # Sample half mesenteric zone into equally spaced nodes - radiansAroundMZ = mzWidth/radius - radiansPerElementAroundMZ = radiansAroundMZ / elementsCountAroundMZ - for n1 in range(int(elementsCountAroundMZ*0.5 + 1)): - radiansAround = radiansPerElementAroundMZ * n1 - cosRadiansAround = math.cos(radiansAround) - sinRadiansAround = math.sin(radiansAround) - x = [ radius*cosRadiansAround, - radius*sinRadiansAround, - 0.0 ] - d1 = [ -radius*sinRadiansAround*radiansPerElementAroundMZ, - radius*cosRadiansAround*radiansPerElementAroundMZ, - 0.0] - xMZ.append(x) - d1MZ.append(d1) - arcLengthPerMZ = radius*radiansPerElementAroundMZ - - # Sample half non-mesenteric zone into equally spaced nodes - halfCircumference = math.pi*radius - xNonMZ, d1NonMZ, arcLengthPerNonMZ, arcLengthPerTransition = sampleHaustrum(nx, nd1, xMZ[-1], d1MZ[-1], halfCircumference, mzWidth*0.5, elementsCountAroundNonMZ) - - xListBaseLayer = xList + xMZ + xNonMZ[1:] - d1ListBaseLayer = d1List + d1MZ + d1NonMZ[1:] - lengthHalfList = len(xListBaseLayer) - - # Reflect to get other half - for n in range(1,int(lengthHalfList)-1): - idx = -n + lengthHalfList - 1 - x = xListBaseLayer[idx] - d1 = d1ListBaseLayer[idx] - xReflect = [x[0], -x[1], x[2]] - d1Reflect = [-d1[0], d1[1], d1[2]] - xListBaseLayer.append(xReflect) - d1ListBaseLayer.append(d1Reflect) - - # Generate node along segment length - lengthPerElementAlong = segmentLength / elementsCountAlongSegment - d2 = [0.0, 0.0, lengthPerElementAlong] - for n2 in range(elementsCountAlongSegment + 1): - for n1 in range(elementsCountAroundNonMZ + elementsCountAroundMZ): - z = lengthPerElementAlong*n2 - x = [xListBaseLayer[n1][0], xListBaseLayer[n1][1], z] - d1 = d1ListBaseLayer[n1] - xList.append(x) - d1List.append(d1) - d2List.append(d2) - - # Calculate uList for elements on outer surface along mid-length of segment - uList, totalArcLengthOuterMidLength = getuListFromOuterMidLengthProfile(xListBaseLayer, d1ListBaseLayer, segmentAxis, wallThickness, transitElementList) - - return annotationGroups, annotationArray, transitElementList, uList, totalArcLengthOuterMidLength, xList, d1List, d2List, segmentAxis diff --git a/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentteniacoli1.py b/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentteniacoli1.py deleted file mode 100644 index 773fa9b1..00000000 --- a/scaffoldmaker/meshtypes/meshtype_3d_colonsegmentteniacoli1.py +++ /dev/null @@ -1,1291 +0,0 @@ -""" -Generates a single 3-D colon segment mesh along a central -line, with variable numbers of elements around, along and -through wall, with variable radius and thickness along. -""" - -import math -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup -from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear -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 -from scaffoldmaker.utils import zinc_utils -from opencmiss.zinc.element import Element, Elementbasis -from opencmiss.zinc.field import Field -from opencmiss.zinc.node import Node - -class MeshType_3d_colonsegmentteniacoli1(Scaffold_base): - ''' - Generates a single 3-D colon segment mesh with variable - numbers of tenia coli, elements around, along the central - line, and through wall. The cross-section profile of the colon - segment varies with species and is dependent on the number - of tenia coli. - Pig: 2 tenia coli, bow tie profile - Human (Default): 3 tenia coli, triangular profile with rounded - corners at the inter-haustral septa, and a clover - profile in the intra-haustral region. - ''' - @staticmethod - def getName(): - return '3D Colon Segment Tenia Coli 1' - - @staticmethod - def getParameterSetNames(): - return [ - 'Default', - 'Human 1', - 'Pig 1'] - - @staticmethod - def getDefaultOptions(parameterSetName='Default'): - options = { - 'Number of elements around tenia coli' : 2, - 'Number of elements around haustrum' : 8, - 'Number of elements along segment' : 4, - 'Number of elements through wall' : 1, - 'Inner radius': 1.0, - 'Corner inner radius factor': 0.5, - 'Haustrum inner radius factor': 0.5, - 'Segment length end derivative factor': 0.5, - 'Segment length mid derivative factor': 2.0, - 'Segment length': 1.5, - 'Number of tenia coli': 3, - 'Tenia coli width': 0.2, - 'Tenia coli thickness': 0.03, - 'Wall thickness': 0.02, - 'Use cross derivatives' : False, - 'Use linear through wall' : True, - 'Refine' : False, - 'Refine number of elements around' : 1, - 'Refine number of elements along segment' : 1, - 'Refine number of elements through wall' : 1 - } - if 'Pig' in parameterSetName: - options['Number of elements along segment'] = 4 - options['Inner radius'] = 0.9 - options['Corner inner radius factor'] = 0.0 - options['Haustrum inner radius factor'] = 0.2 - options['Number of tenia coli'] = 2 - options['Tenia coli width'] = 0.5 - options['Tenia coli thickness'] = 0.05 - options['Wall thickness'] = 0.2 - return options - - @staticmethod - def getOrderedOptionNames(): - return [ - 'Number of elements around tenia coli', - 'Number of elements around haustrum', - 'Number of elements along segment', - 'Number of elements through wall', - 'Inner radius', - 'Corner inner radius factor', - 'Haustrum inner radius factor', - 'Segment length end derivative factor', - 'Segment length mid derivative factor', - 'Segment length', - 'Number of tenia coli', - 'Tenia coli width', - 'Tenia coli thickness', - 'Wall thickness', - 'Use cross derivatives', - 'Use linear through wall', - 'Refine', - 'Refine number of elements around', - 'Refine number of elements along segment', - 'Refine number of elements through wall' - ] - - @staticmethod - def checkOptions(options): - for key in [ - 'Number of elements through wall', - 'Refine number of elements around', - 'Refine number of elements along segment', - 'Refine number of elements through wall']: - if options[key] < 1: - options[key] = 1 - for key in [ - 'Number of elements around tenia coli', - 'Number of elements along segment']: - if options[key] < 2: - options[key] = 2 - if options['Number of elements around haustrum'] < 4: - options['Number of elements around haustrum'] = 4 - for key in [ - 'Number of elements around tenia coli', - 'Number of elements around haustrum']: - if options[key] % 2 > 0: - options[key] = options[key] + 1 - for key in [ - 'Inner radius', - 'Haustrum inner radius factor', - 'Segment length end derivative factor', - 'Segment length mid derivative factor', - 'Segment length', - 'Tenia coli thickness', - 'Wall thickness']: - if options[key] < 0.0: - options[key] = 0.0 - if options['Corner inner radius factor'] < 0.1: - options['Corner inner radius factor'] = 0.1 - for key in [ - 'Corner inner radius factor', - 'Segment length end derivative factor']: - if options[key] > 1.0: - options[key] = 1.0 - if options['Number of tenia coli'] < 2: - options['Number of tenia coli'] = 2 - elif options['Number of tenia coli'] > 3: - options['Number of tenia coli'] = 3 - if options['Tenia coli width'] < 0.2*options['Inner radius']: - options['Tenia coli width'] = round(0.2*options['Inner radius'], 2) - if options['Tenia coli width'] > round(math.sqrt(3)*0.5*options['Inner radius'],2): - options['Tenia coli width'] = round(math.sqrt(3)*0.5*options['Inner radius'],2) - - @staticmethod - def generateBaseMesh(region, options): - """ - Generate the base tricubic Hermite mesh. See also generateMesh(). - :param region: Zinc region to define model in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - :return: None - """ - elementsCountAroundTC = options['Number of elements around tenia coli'] - elementsCountAroundHaustrum = options['Number of elements around haustrum'] - elementsCountAlongSegment = options['Number of elements along segment'] - elementsCountThroughWall = options['Number of elements through wall'] - radius = options['Inner radius'] - cornerInnerRadiusFactor = options['Corner inner radius factor'] - haustrumInnerRadiusFactor = options['Haustrum inner radius factor'] - segmentLengthEndDerivativeFactor = options['Segment length end derivative factor'] - segmentLengthMidDerivativeFactor = options['Segment length mid derivative factor'] - segmentLength = options['Segment length'] - tcCount = options['Number of tenia coli'] - tcWidth = options['Tenia coli width'] - tcThickness = options['Tenia coli thickness'] - wallThickness = options['Wall thickness'] - useCrossDerivatives = options['Use cross derivatives'] - useCubicHermiteThroughWall = not(options['Use linear through wall']) - haustraSegmentCount = 1 - elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount - - cx = [ [ 0.0, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] - cd1 = [ [ segmentLength, 0.0, 0.0 ], [ segmentLength, 0.0, 0.0 ] ] - cd2 = [ [ 0.0, 1.0, 0.0 ], [ 0.0, 1.0, 0.0 ] ] - cd12 = [ [0.0, 0.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ] - - # Generate inner surface of a colon segment - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength, xHaustraInner, d1HaustraInner, d2HaustraInner, haustraSegmentAxis = getColonSegmentInnerPointsTeniaColi(region, elementsCountAroundTC, - elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, tcWidth, radius, cornerInnerRadiusFactor, haustrumInnerRadiusFactor, - segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, segmentLength, wallThickness) - - # Generate tube mesh - annotationGroups, nextNodeIdentifier, nextElementIdentifier, xList, d1List, d2List, d3List, sx, curvatureAlong, factorList = tubemesh.generatetubemesh(region, - elementsCountAround, elementsCountAlongSegment, elementsCountThroughWall, haustraSegmentCount, cx, cd1, cd2, cd12, - xHaustraInner, d1HaustraInner, d2HaustraInner, wallThickness, haustraSegmentAxis, segmentLength, - useCrossDerivatives, useCubicHermiteThroughWall, annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength) - - # Generate tenia coli - annotationGroupsTC, nextNodeIdentifier, nextElementIdentifier = getTeniaColi(region, nextNodeIdentifier, nextElementIdentifier, - useCrossDerivatives, useCubicHermiteThroughWall, xList, d1List, d2List, d3List, - elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, elementsCountThroughWall, wallThickness, - tcWidth, tcThickness, sx, curvatureAlong, factorList, uList, arcLengthOuterMidLength, segmentLength, tcCount) - - annotationGroups += annotationGroupsTC - - return annotationGroups - - @classmethod - def generateMesh(cls, region, options): - """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - """ - if not options['Refine']: - return cls.generateBaseMesh(region, options) - - refineElementsCountAround = options['Refine number of elements around'] - refineElementsCountAlong = options['Refine number of elements along segment'] - refineElementsCountThroughWall = options['Refine number of elements through wall'] - - baseRegion = region.createRegion() - baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) - - meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) - meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) - return meshrefinement.getAnnotationGroups() - -def getColonSegmentInnerPointsTeniaColi(region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment, tcCount, tcWidth, radius, - cornerInnerRadiusFactor, haustrumInnerRadiusFactor, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor, segmentLength, wallThickness): - """ - Generates a 3-D colon segment mesh with two or three tenia coli with variable - numbers of elements around, along the central line, and through wall. - Colon segment with two tenia coli (pig) has a circular profile at the - inter-haustral septa, and a bowtie profile in the intra-haustral region. - Colon segment with three tenia coli (human) has a triangular profile - with rounded corners at the inter-haustral septa, and a clover profile - in the intra-haustral region. - :param elementsCountAroundTC: Number of elements around each tenia coli. - :param elementsCountAroundHaustrum: Number of elements around haustrum. - :param elementsCountAlongSegment: Number of elements along colon segment. - :param tcCount: Number of tenia coli. - :param tcWidth: Width of tenia coli. - :param radius: Inner radius defined from center of triangular - profile to vertex of the triangle. - :param cornerInnerRadiusFactor: Roundness of triangular corners of - 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 - radius to obtain radius of intersecting circles in the middle cross-section - along a haustra segment. - :param segmentLengthEndDerivativeFactor: Factor is multiplied by segment - length to scale derivative along the end of a segment length. - :param segmentLengthMidDerivativeFactor: Factor is multiplied by segment - length to scale derivative along the mid length of the segment. - :param segmentLength: Length of a colon segment. - :param wallThickness: Thickness of wall. - :return annotationGroups, annotationArray: annotationArray stores annotation - names of elements around - :return transitElementList: stores true if element around is an element that - transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. - :return uList: List of xi for each node around mid-length haustra. - : return totalArcLengthOuterMidLengthHaustra: Total arclength around outer - surface on mid-length haustra. - :return coordinates, derivatives on inner surface of a colon segment. - """ - annotationGroups = [] - annotationArray = [] - - transitElementListHaustrum = [0]*int(elementsCountAroundTC*0.5) + [1] + [0]*int(elementsCountAroundHaustrum - 2) + [1] + [0]*int(elementsCountAroundTC*0.5) - transitElementList = transitElementListHaustrum * tcCount - - # create nodes - x = [ 0.0, 0.0, 0.0 ] - dx_ds1 = [ 0.0, 0.0, 0.0 ] - dx_ds2 = [ 0.0, 0.0, 0.0 ] - dx_ds3 = [ 0.0, 0.0, 0.0 ] - sampleElementOut = 20 - segmentAxis = [0.0, 0.0, 1.0] - haustrumRadius = (haustrumInnerRadiusFactor + 1)*radius - - xAround = [] - d1Around = [] - xHalfSetInterHaustra = [] - d1HalfSetInterHaustra = [] - nxHaustrum = [] - nd1Haustrum = [] - xHalfSetIntraHaustra = [] - d1HalfSetIntraHaustra = [] - xInnerRaw = [] - dx_ds2InnerRaw = [] - xInnerList = [] - dx_ds2InnerList = [] - xFinal = [] - d1Final = [] - d2Final = [] - uList = [] - - # Inter-haustral segment - # Set up profile - if tcCount == 2: # Circular profile - xLoop, d1Loop = createCirclePoints([ 0.0, 0.0, 0.0 ], [ radius, 0.0, 0.0 ], [ 0.0, radius, 0.0 ], sampleElementOut, startRadians = 0.0) - - else: # tcCount == 3, Triangular profile - cornerRC = cornerInnerRadiusFactor*radius - radiansRangeRC = [7*math.pi/4, 0.0, math.pi/4] - - for n1 in range(3): - radiansAround = n1*2.0*math.pi / 3.0 - cosRadiansAround = math.cos(radiansAround) - sinRadiansAround = math.sin(radiansAround) - xc = [(radius - cornerRC) * cosRadiansAround, (radius - cornerRC) * sinRadiansAround, 0.0] - - for n in range(3): - radiansRC = radiansAround + radiansRangeRC[n] - cosRadiansRC = math.cos(radiansRC) - sinRadiansRC = math.sin(radiansRC) - x = [xc[0] + cornerRC*cosRadiansRC, xc[1] + cornerRC*sinRadiansRC, 0.0] - xAround.append(x) - d1 = [ cornerRC*math.pi/4.0 * -sinRadiansRC, cornerRC*math.pi/4.0 * cosRadiansRC, 0.0] - d1Around.append(d1) - - xSample = xAround[1:9] +[xAround[0], xAround[1]] - d1Sample = d1Around[1:9] +[d1Around[0], d1Around[1]] - sx, sd1, se, sxi, _= interp.sampleCubicHermiteCurves(xSample, d1Sample, sampleElementOut) - xLoop = sx[:-1] - d1Loop = interp.smoothCubicHermiteDerivativesLoop(sx[:-1], sd1[:-1]) - - # Calculate arc length - arcLength = 0.0 - arcDistance = interp.getCubicHermiteArcLength(xLoop[0], d1Loop[0], xLoop[1], d1Loop[1]) - arcLength = arcDistance * sampleElementOut - arcStart = 0.0 - numberOfHalves = tcCount * 2 - arcEnd = arcLength/numberOfHalves - - # Find edge of TC - arcDistanceTCEdge = findEdgeOfTeniaColi(xLoop, d1Loop, tcWidth, arcStart, arcEnd) - - # Sample TC into equally sized elements - xTC, d1TC, _ = sampleTeniaColi(xLoop, d1Loop, arcDistanceTCEdge, elementsCountAroundTC) - - # Sample haustrum into equally sized elements - xHaustrum, d1Haustrum, _, _ = sampleHaustrum(xLoop, d1Loop, xTC[-1], d1TC[-1], arcLength/numberOfHalves, arcDistanceTCEdge, elementsCountAroundHaustrum) - - xHalfSetInterHaustra = xHalfSetInterHaustra + xTC + xHaustrum[1:] - d1HalfSetInterHaustra = d1HalfSetInterHaustra + d1TC + d1Haustrum[1:] - - # Intra-haustral segment - # Set up profile - if tcCount == 2: # Bow-tie profile - originRC = (radius*radius - haustrumRadius*haustrumRadius)/(-2.0*haustrumRadius) - RC = haustrumRadius - originRC - thetaStart = math.asin(-originRC/ RC) - thetaEnd = math.pi - thetaStart - rotOriginRC = [ 0.0, -originRC, 0.0 ] - - else: # tcCount == 3, Clover profile - xc = [(radius - cornerRC)* math.cos(0.0), (radius - cornerRC)*math.sin(0.0), 0.0] - pt1 = [xc[0] + cornerRC*math.cos(0.0), xc[1] + cornerRC*math.sin(0.0), 0.0] - xTC2 = radius* math.cos(2.0*math.pi/3.0) - yTC2 = radius* math.sin(2.0*math.pi/3.0) - originRC = (xTC2*xTC2 + yTC2*yTC2 - haustrumRadius*haustrumRadius) / (2*(-xTC2 - haustrumRadius)) - RC = haustrumRadius - originRC - - # Rotate to find originRC of 1st haustrum - yTC1 = pt1[1] - rotOriginRC = [ originRC*math.cos(-2.0/3.0*math.pi), originRC*math.sin(-2.0/3.0*math.pi), 0.0] - - thetaStart = math.asin((yTC1 + rotOriginRC[1]) / RC) - thetaEnd = math.pi - math.asin((yTC2 + rotOriginRC[1])/ RC) - - thetaHalfHaustrum = (thetaEnd + thetaStart)*0.5 - concaveFactor = 0.15 - thetaConcave = (thetaEnd - thetaStart)*concaveFactor + thetaStart - thetaCircularStart = (thetaEnd - thetaStart)*(concaveFactor + 0.5)*0.5 + thetaStart - - thetaSet = [thetaStart, thetaConcave] - xConcave, _ = getCircleXandD1FromRadians(thetaSet, RC, rotOriginRC) - d1 = [0.0, xConcave[1][1], 0.0] - nxHaustrum = nxHaustrum + xConcave - nd1Haustrum = nd1Haustrum + [d1, d1] - thetaCircularSet = [thetaCircularStart, thetaHalfHaustrum] - xCircular, d1Circular = getCircleXandD1FromRadians(thetaCircularSet, RC, rotOriginRC) - nxHaustrum = nxHaustrum + xCircular - nd1Haustrum = nd1Haustrum + d1Circular - smoothd1 = interp.smoothCubicHermiteDerivativesLine(nxHaustrum, nd1Haustrum, fixStartDirection = True, fixEndDirection = True) - - # Find max x along path - sxHaustrum, sd1Haustrum, _, _ , _ = interp.sampleCubicHermiteCurves(nxHaustrum, smoothd1, sampleElementOut) - xVal = [] - for i in range(len(sxHaustrum)): - xVal.append(sxHaustrum[i][0]) - maxValue = max(xVal) - maxIndex = xVal.index(maxValue) - yAtMaxValue = sxHaustrum[maxIndex][1] - - arcLength = 0.0 - for e in range(len(nxHaustrum)-1): - arcDistance = interp.getCubicHermiteArcLength(nxHaustrum[e], smoothd1[e], nxHaustrum[e+1], smoothd1[e+1]) - arcLength = arcLength + arcDistance - - arcDistanceAtMaxValue = arcLength / sampleElementOut * maxIndex - arcStart = 0.0 if yAtMaxValue > tcWidth*0.5 else arcDistanceAtMaxValue - arcEnd = arcDistanceAtMaxValue if yAtMaxValue > tcWidth*0.5 else arcLength - - # Find edge of TC - arcDistanceTCEdge = findEdgeOfTeniaColi(nxHaustrum, smoothd1, tcWidth, arcStart, arcEnd) - - # Sample TC into equally sized elements - xTC, d1TC, arcLengthPerTC = sampleTeniaColi(nxHaustrum, smoothd1, arcDistanceTCEdge, elementsCountAroundTC) - - # Sample haustrum into equally sized elements - xHaustrum, d1Haustrum, arcLengthPerHaustrum, arcLengthPerTransition = sampleHaustrum(nxHaustrum, smoothd1, xTC[-1], d1TC[-1], arcLength, arcDistanceTCEdge, elementsCountAroundHaustrum) - - xHalfSetIntraHaustra = xHalfSetIntraHaustra + xTC + xHaustrum[1:] - d1HalfSetIntraHaustra = d1HalfSetIntraHaustra + d1TC + d1Haustrum[1:] - - # Calculate uList for elements lying around outer surface along mid-length haustra - d2 = [] - for n in range(len(xHalfSetIntraHaustra)): - d2.append([0.0, 0.0, 0.0]) - xInnerMidLengthHaustra, d1InnerMidLengthHaustra, _ = getFullProfileFromHalfHaustrum(xHalfSetIntraHaustra, d1HalfSetIntraHaustra, d2, tcCount) - uList, totalArcLengthOuterMidLengthHaustra = getuListFromOuterMidLengthProfile(xInnerMidLengthHaustra, d1InnerMidLengthHaustra, segmentAxis, wallThickness, transitElementList) - - # Sample arclength of haustra segment - elementsCountAroundHalfHaustrum = int((elementsCountAroundTC + elementsCountAroundHaustrum)*0.5) - - for n1 in range(elementsCountAroundHalfHaustrum + 1): - v1 = [xHalfSetInterHaustra[n1][0], xHalfSetInterHaustra[n1][1], 0.0] - startArcLength = segmentLengthEndDerivativeFactor * segmentLength - d1 = [ c*startArcLength for c in segmentAxis] - v2 = [xHalfSetIntraHaustra[n1][0], xHalfSetIntraHaustra[n1][1], segmentLength/2] - midArcLength = segmentLengthMidDerivativeFactor * segmentLength - d2 = [c*midArcLength for c in segmentAxis] - v3 = [xHalfSetInterHaustra[n1][0], xHalfSetInterHaustra[n1][1], segmentLength] - d3 = [ c*startArcLength for c in segmentAxis] - nx = [v1, v2, v3] - nd1 = [d1, d2, d3] - sx, sd1, se, sxi, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAlongSegment) - xInnerRaw.append(sx) - dx_ds2InnerRaw.append(sd1) - - # Re-arrange sample order & calculate dx_ds1 and dx_ds3 from dx_ds2 - for n2 in range(elementsCountAlongSegment + 1): - xAround = [] - unitdx_ds1Around = [] - d2Around = [] - - for n1 in range(elementsCountAroundHalfHaustrum+1): - x = xInnerRaw[n1][n2] - xInnerList.append(x) - dx_ds2 = dx_ds2InnerRaw[n1][n2] - dx_ds2InnerList.append(dx_ds2) - unitTangent = vector.normalise(dx_ds2) - # Intra-Haustra segments - if n1 == 0: - unitdx_ds1 = vector.normalise(d1HalfSetIntraHaustra[n1]) - else: # points on clover - if n2 <= int(elementsCountAlongSegment/2): # first half of segmentLength - axisRot = vector.crossproduct3(segmentAxis, unitTangent) - elif n2 > int(elementsCountAlongSegment/2): # second half of segmentLength - axisRot = vector.crossproduct3(unitTangent, segmentAxis) - rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, math.pi/2) - rotNormal = [rotFrame[j][0]*unitTangent[0] + rotFrame[j][1]*unitTangent[1] + rotFrame[j][2]*unitTangent[2] for j in range(3)] - unitdx_ds3 = vector.normalise(rotNormal) - unitdx_ds1 = vector.crossproduct3(unitTangent, unitdx_ds3) - xAround.append(x) - d2Around.append(dx_ds2) - unitdx_ds1Around.append(unitdx_ds1) - - if n2 > 0 and n2 < elementsCountAlongSegment: - dx_ds1InnerAroundList = [] - if elementsCountAlongSegment%2 == 0 and n2 == int(elementsCountAlongSegment*0.5): - dx_ds1InnerAroundList = dx_ds1InnerAroundList + d1HalfSetIntraHaustra - else: - for n1 in range(elementsCountAroundHalfHaustrum): - v1 = xAround[n1] - d1 = unitdx_ds1Around[n1] - v2 = xAround[n1+1] - d2 = unitdx_ds1Around[n1+1] - arcLengthAround = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) - dx_ds1 = [c*arcLengthAround for c in d1] - dx_ds1InnerAroundList.append(dx_ds1) - # Account for d1 of node sitting on half haustrum - dx_ds1 = [c*arcLengthAround for c in unitdx_ds1Around[elementsCountAroundHalfHaustrum]] - dx_ds1InnerAroundList.append(dx_ds1) - d1Smoothed = interp.smoothCubicHermiteDerivativesLine(xAround, dx_ds1InnerAroundList, fixStartDerivative = True) - d1TCEdge = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC*0.5)], vector.magnitude(d1Smoothed[int(elementsCountAroundTC*0.5 - 1)])) - d1Transition = vector.setMagnitude(d1Smoothed[int(elementsCountAroundTC*0.5 + 1)], vector.magnitude(d1Smoothed[int(elementsCountAroundTC*0.5 + 2)])) - d1Corrected = [] - d1Corrected = d1Corrected + d1Smoothed[:int(elementsCountAroundTC*0.5)] - d1Corrected.append(d1TCEdge) - d1Corrected.append(d1Transition) - d1Corrected = d1Corrected + d1Smoothed[int(elementsCountAroundTC*0.5 + 2):] - else: - d1Corrected = d1HalfSetInterHaustra - - xAlongList, d1AlongList, d2AlongList = getFullProfileFromHalfHaustrum(xAround, d1Corrected, d2Around, tcCount) - xFinal = xFinal + xAlongList - d1Final = d1Final + d1AlongList - d2Final = d2Final + d2AlongList - - return annotationGroups, annotationArray, transitElementList, uList, totalArcLengthOuterMidLengthHaustra, xFinal, d1Final, d2Final, segmentAxis - -def findEdgeOfTeniaColi(nx, nd1, tcWidth, arcStart, arcEnd): - """ - Locate edge of tenia coli on a cubic hermite interpolated - curve defined by nodes nx and derivatives nd1. - :param nx: Coordinates of nodes along curve. - :param nd1: Derivatives of nodes along curve. - :param tcWidth: Width of tenia coli. - :param arcStart: Lower limit of arc distance to search for - edge of tenia coli. - :param arcEnd: Upper limit of arc distance to search for - edge of tenia coli. - :return: arc distance covered by tenia coli. - """ - xTol = 1.0E-6 - for iter in range(100): - arcDistance = (arcStart + arcEnd)*0.5 - x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) - diff = x[1] - tcWidth*0.5 - if abs(diff) > xTol: - if diff < 0.0: - arcStart = arcDistance - else: - arcEnd = arcDistance - else: - arcDistanceTCEdge = arcDistance - break - if iter > 99: - print('Search for TC boundary - Max iters reached:',iter) - - return arcDistanceTCEdge - -def sampleTeniaColi(nx, nd1, arcDistanceTCEdge, elementsCountAroundTC): - """ - Get systematically spaced points and derivatives over the width of - tenia coli over a cubic hermite interpolated curves defined by nodes - nx and derivatives nd1 lying along half of a haustrum. - :param nx: Coordinates of nodes along curve. - :param nd1: Derivatives of nodes along curve. - :param arcDistanceTCEdge: Arc distance covered by tenia coli. - :param elementsCountAroundTC: Number of elements around tenia coli. - :return: coordinates, derivatives of tenia coli, and arclength of - each tenia coli element. - """ - xTC = [] - d1TC = [] - arcDistancePerElementTC = arcDistanceTCEdge / (elementsCountAroundTC*0.5) - for e in range(int(elementsCountAroundTC*0.5)+1): - arcDistance = arcDistancePerElementTC * e - x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) - d1Scaled = vector.setMagnitude(d1, arcDistancePerElementTC) - xTC.append(x) - d1TC.append(d1Scaled) - - return xTC, d1TC, arcDistancePerElementTC - -def sampleHaustrum(nx, nd1, xTCLast, d1TCLast, arcLength, arcDistanceTCEdge, elementsCountAroundHaustrum): - """ - Get systematically spaced points and derivatives over the width of - haustrum over a cubic hermite interpolated curves defined by nodes - nx and derivatives nd1 lying along half of a haustrum. - :param nx: Coordinates of nodes along curve. - :param nd1: Derivatives of nodes along curve. - :param xTCLast: Coordinates of node lying on edge of tenia coli. - :param d1TCLast: Derivative of node lying on edge of tenia coli. - :param arcLength: Arc length of curve. - :param arcDistanceTCEdge: Arc distance covered by tenia coli. - :param elementsCountAroundHaustrum: Number of elements around haustrum. - :return: coordinates, derivatives on haustrum. - :return: arclength of haustrum element and transition elements. - """ - xHaustrum = [] - d1Haustrum = [] - elementLengths = [] - length = arcLength - arcDistanceTCEdge - elementsCountOut = int(elementsCountAroundHaustrum * 0.5) - addLengthStart = 0.5 * vector.magnitude(d1TCLast) - lengthFractionStart = 0.5 - proportionStart = 1.0 - elementLengthMid = (length - addLengthStart) / (elementsCountOut - 1.0 + proportionStart*lengthFractionStart ) - elementLengthProportionStart = proportionStart*lengthFractionStart*elementLengthMid - elementLengthProportionEnd = elementLengthMid - - for e in range(elementsCountOut): - xi = e/(elementsCountOut - 1) - elementLengths.append(elementLengthMid) - elementLengths[0] = addLengthStart + elementLengthProportionStart - elementLengths[-1] = 0.0 + elementLengthProportionEnd - haustrumElementLength = elementLengthMid - transitionElementLength = elementLengths[0] - - arcDistance = arcDistanceTCEdge - xHaustrum.append(xTCLast) - d1Scaled = vector.setMagnitude(d1TCLast, elementLengths[0]) - d1Haustrum.append(d1Scaled) - - for e in range(elementsCountOut): - arcDistance = arcDistance + elementLengths[e] - x, d1, _, _ = interp.getCubicHermiteCurvesPointAtArcDistance(nx, nd1, arcDistance) - d1Scaled = vector.setMagnitude(d1, elementLengths[e] if e > 0 else elementLengths[e+1]) - xHaustrum.append(x) - d1Haustrum.append(d1Scaled) - - return xHaustrum, d1Haustrum, haustrumElementLength, transitionElementLength - -def getCircleXandD1FromRadians(thetaSet, radius, origin): - """ - Gets the coordinates and derivatives along a circular path - based on an angular range. - :param thetaSet: Lower and upper limit of theta. - :param radius: Radius of circle. - :param origin: Origin of circle. - :return: coordinates, derivatives on lower and upper - limit of thetaSet. - """ - nx = [] - nd1 = [] - dTheta = thetaSet[1] - thetaSet[0] - for n in range(2): - theta = thetaSet[n] - x = [radius*math.cos(theta) - origin[0], - radius*math.sin(theta) - origin[1], - 0.0] - d1 = [-radius*math.sin(theta)*dTheta, - radius*math.cos(theta)*dTheta, - 0.0] - nx.append(x) - nd1.append(d1) - - return nx, nd1 - -def getFullProfileFromHalfHaustrum(xHaustrumHalfSet, d1HaustrumHalfSet, d2HaustrumHalfSet, tcCount): - """ - Gets the coordinates and derivatives of the entire profile - using points from first half of the first sector. The first - sector starts from the x-axis. The first half set of points - are reflected across the x-axis followed by rotation to get - the points in the second half of the first sector. The full - set of points in the first sector are then rotated to obtain - points in the other two sectors. - :param xHaustrumHalfSet: Coordinates of points in first - half of the first sector. - :param d1HaustrumHalfSet: Derivatives of points around first - half of the first sector. - :param d2HaustrumHalfSet: Derivatives of points along first - half of the first sector. - :param tcCount: Number of tenia coli. - :return: coordinates, derivatives of points over entire profile. - """ - xHaustrumHalfSet2 = [] - d1HaustrumHalfSet2 = [] - d2HaustrumHalfSet2 = [] - xHaustrum = [] - d1Haustrum = [] - d2Haustrum = [] - xHaustra = [] - d1Haustra = [] - d2Haustra = [] - rotAng = 2/3*math.pi if tcCount == 3 else math.pi - - for n in range(1,len(xHaustrumHalfSet)): - idx = -n + len(xHaustrumHalfSet) - 1 - x = xHaustrumHalfSet[idx] - d1 = d1HaustrumHalfSet[idx] - xReflect = [x[0], -x[1], x[2]] - d1Reflect = [d1[0], -d1[1], d1[2]] - xRot = [xReflect[0]*math.cos(rotAng) - xReflect[1]*math.sin(rotAng), - xReflect[0]*math.sin(rotAng) + xReflect[1]*math.cos(rotAng), - xReflect[2]] - d1Rot = [-(d1Reflect[0]*math.cos(rotAng) - d1Reflect[1]*math.sin(rotAng)), - -(d1Reflect[0]*math.sin(rotAng) + d1Reflect[1]*math.cos(rotAng)), - -d1Reflect[2]] - d2 = d2HaustrumHalfSet[idx] - d2Reflect = [d2[0], -d2[1], d2[2]] - d2Rot = [(d2Reflect[0]*math.cos(rotAng) - d2Reflect[1]*math.sin(rotAng)), - (d2Reflect[0]*math.sin(rotAng) + d2Reflect[1]*math.cos(rotAng)), - d2Reflect[2]] - xHaustrumHalfSet2.append(xRot) - d1HaustrumHalfSet2.append(d1Rot) - d2HaustrumHalfSet2.append(d2Rot) - - xHaustrum = xHaustrumHalfSet + xHaustrumHalfSet2 - d1Haustrum = d1HaustrumHalfSet + d1HaustrumHalfSet2 - d2Haustrum = d2HaustrumHalfSet + d2HaustrumHalfSet2 - - # Rotate to get all 3 sectors - xHaustra = xHaustra + xHaustrum[:-1] - d1Haustra = d1Haustra + d1Haustrum[:-1] - d2Haustra = d2Haustra + d2Haustrum[:-1] - - ang = [ 2/3*math.pi, -2/3*math.pi] if tcCount == 3 else [math.pi] - for i in range(tcCount - 1): - rotAng = ang[i] - cosRotAng = math.cos(rotAng) - sinRotAng = math.sin(rotAng) - for n in range(len(xHaustrum)- 1): - x = xHaustrum[n] - d1 = d1Haustrum[n] - d2 = d2Haustrum[n] - x = [ x[0]*cosRotAng - x[1]*sinRotAng, x[0]*sinRotAng + x[1]*cosRotAng, x[2]] - xHaustra.append(x) - dx_ds1 = [ d1[0]*cosRotAng - d1[1]*sinRotAng, d1[0]*sinRotAng + d1[1]*cosRotAng, d1[2]] - d1Haustra.append(dx_ds1) - dx_ds2 = [ d2[0]*cosRotAng - d2[1]*sinRotAng, d2[0]*sinRotAng + d2[1]*cosRotAng, d2[2]] - d2Haustra.append(dx_ds2) - - return xHaustra, d1Haustra, d2Haustra - -def getuListFromOuterMidLengthProfile(xInnerMidLength, d1InnerMidLength, - segmentAxis, wallThickness, transitElementList): - """ - Gets a list of xi for flat and texture coordinates calculated - from outer arclength of elements around mid-length of a segment. - :param xInnerMidLength: Coordinates of points on inner surface around - mid-length of segment. - :param d1InnerMidLength: Derivatives of points on inner surface around - mid-length of segment. - :param segmentAxis: Axis of segment. - :param wallThickness: Thickness of wall. - :param transitElementList: stores true if element around is an element that - transits from tenia coli / mesenteric zone to haustrum / non-mesenteric zone. - :return uList: List containing xi for each point on outer surface extruded by - wall thickness from points on inner surface. - :return totalArcLengthOuterMidLength: Total arclength around outer surface on - mid-length section of segment. - """ - unitNormList = [] - d1OuterMidLength = [] - - for n in range(len(xInnerMidLength)): - unitNormList.append(vector.normalise(vector.crossproduct3(d1InnerMidLength[n], segmentAxis))) - - xOuterMidLength, curvatureInnerMidLength = tubemesh.getOuterCoordinatesAndCurvatureFromInner(xInnerMidLength, d1InnerMidLength, unitNormList, wallThickness, 0, len(xInnerMidLength), transitElementList) - - for n1 in range(len(xOuterMidLength)): - factor = 1.0 + wallThickness * curvatureInnerMidLength[n1] - dx_ds1 = [ factor*c for c in d1InnerMidLength[n1]] - d1OuterMidLength.append(dx_ds1) - - arcLengthList = [] - for n1 in range(len(xOuterMidLength)): - arcLengthPerElement = interp.getCubicHermiteArcLength(xOuterMidLength[n1], d1OuterMidLength[n1], xOuterMidLength[(n1+1) % len(xOuterMidLength)], d1OuterMidLength[(n1+1)% len(xOuterMidLength)]) - arcLengthList.append(arcLengthPerElement) - - # Total arcLength - totalArcLengthOuterMidLength = 0.0 - for n1 in range(len(arcLengthList)): - totalArcLengthOuterMidLength += arcLengthList[n1] - - uList = [0.0] - arcDistance = 0 - for n in range(len(arcLengthList)): - arcDistance = arcDistance + arcLengthList[n] - xi = arcDistance / totalArcLengthOuterMidLength - uList.append(xi) - - return uList, totalArcLengthOuterMidLength - -def getTeniaColi(region, nodeIdentifier, elementIdentifier, useCrossDerivatives, - useCubicHermiteThroughWall, xList, d1List, d2List, d3List, - elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlong, elementsCountThroughWall, - wallThickness, tcWidth, tcThickness, sxCentralLine, curvatureAlong, factorList, uList, - arcLengthOuterMidLength, totalLengthAlong, tcCount): - """ - Create equally spaced nodes and elements for tenia coli over the outer - surface of the haustra. Nodes of the tenia coli is sampled from a cubic - hermite curve running through the left edge of tenia coli boundary on the - outer surface of the haustra, a midpoint lying at a distance of tenia - coli thickness above the midpoint of the boundary of tenia coli, and the - right edge of tenia coli boundary. - :param nodeIdentifier, elementIdentifier: First node and element identifier to - use for tenia coli. - :param xList, d1List, d2List, d3List: Coordinates and derivatives of nodes on haustra. - :param elementsCountAroundTC: Number of elements around tenia coli. - :param elementsCountAroundHaustrum: Number of elements around haustrum. - :param elementsCountAlong: Number of elements along scaffold. - :param elementsCountThroughWall: Number of elements through wall. - :param wallThickness: Thickness of wall. - :param tcWidth: Width of tenia coli. - :param tcThickness: Thickness of tenia coli at its thickest part. - :param sxCentralLine: Coordinates sampled from central line. - :param curvatureAlong: Curvatures along the colon for nodes on inner surface of colon. - :param factorList: Factors used for scaling d2 to account for curvature along colon. - :param uList: List of xi for each node around mid-length haustra. - :param arcLengthOuterMidLength: Total arclength of elements around outer surface along - mid-length of segment. - :param totalLengthAlong: Total length of colon along center line. - :param tcCount: Number of tenia coli. - :return: annotationGroups, nodeIdentifier, elementIdentifier - """ - - fm = region.getFieldmodule() - fm.beginChange() - cache = fm.createFieldcache() - coordinates = zinc_utils.getOrCreateCoordinateField(fm) - - if tcCount == 3: - tlGroup = AnnotationGroup(region, 'tenia libera', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - tmGroup = AnnotationGroup(region, 'tenia mesocolica', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - toGroup = AnnotationGroup(region, 'tenia omentalis', FMANumber = 'FMANumber unknown', lyphID = 'Lyph ID unknown') - annotationGroups = [tlGroup, tmGroup, toGroup] - else: - annotationGroups = [ ] - - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - nodetemplate = nodes.createNodetemplate() - nodetemplate.defineField(coordinates) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - if useCubicHermiteThroughWall: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) - - mesh = fm.findMeshByDimension(3) - - if useCubicHermiteThroughWall: - eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives) - else: - eftfactory = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) - eft = eftfactory.createEftBasic() - - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplate.defineField(coordinates, -1, eft) - - # Tenia coli edge elements - elementtemplate1 = mesh.createElementtemplate() - elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - eft1 = eftfactory.createEftWedgeXi1One() - elementtemplate1.defineField(coordinates, -1, eft1) - - elementtemplate2 = mesh.createElementtemplate() - elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - eft2 = eftfactory.createEftWedgeXi1Zero() - elementtemplate2.defineField(coordinates, -1, eft2) - - # create nodes - elementsCountAround = (elementsCountAroundTC + elementsCountAroundHaustrum)*tcCount - TCEdgeFactor = 1.5 - zero = [0.0, 0.0, 0.0] - prevNodeIdentifier = nodeIdentifier - 1 - prevElementIdentifier = elementIdentifier - 1 - - xTCOuterList = [] - d1TCOuterList = [] - d2TCOuterList = [] - d3TCOuterList = [] - bniList = [] - bniTCList = [] - idxList = [] - wList = [] - dwList = [] - - set1 = list(range(int(elementsCountAroundTC/2)+1)) - set2 = list(range(set1[-1] + elementsCountAroundHaustrum, set1[-1] + elementsCountAroundHaustrum + elementsCountAroundTC + 1)) - if tcCount == 3: - set3 = list(range(set2[-1] + elementsCountAroundHaustrum, set2[-1] + elementsCountAroundHaustrum + elementsCountAroundTC +1)) - set4 = list(range(set3[-1] + elementsCountAroundHaustrum, set3[-1] + elementsCountAroundHaustrum + int(elementsCountAroundTC/2))) - setTCIdx = set1 + set2 + set3 + set4 - else: - set3 = list(range(set2[-1] + elementsCountAroundHaustrum, set2[-1] + elementsCountAroundHaustrum + int(elementsCountAroundTC/2))) - setTCIdx = set1 + set2 + set3 - - for n2 in range(elementsCountAlong + 1): - xTCRaw = [] - d1TCRaw = [] - d2TCRaw = [] - d3TCRaw = [] - for N in range(tcCount): - idxTCMid = elementsCountThroughWall*(elementsCountAlong+1)*elementsCountAround + n2*elementsCountAround + 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) - v1 = xList[TCStartIdx] - v2 = xMid - d1MidScaled = [c*tcWidth*TCEdgeFactor for c in vector.normalise(d1Mid)] - v3 = xList[TCEndIdx] - nx = [v1, v2, v3] - nd1 = [d1List[TCStartIdx], d1MidScaled, d1List[TCEndIdx]] - sx, sd1, se, sxi, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) - xTCRaw = xTCRaw + sx[1:-1] - if elementsCountAroundTC == 2: - p = [v2[i] - v1[i] for i in range(3)] - A = vector.dotproduct(unitNorm, p) # A<0 if v2 is higher than v1 - d1 = [c*tcWidth*0.5 for c in vector.normalise(d1Mid)] if A < 0 else d1MidScaled - d1TCRaw = d1TCRaw + sd1[1:-1] if elementsCountAroundTC > 2 else d1TCRaw + [d1] - xTCInnerSet = list(range(TCStartIdx+1, TCEndIdx)) if N > 0 else list(range(TCStartIdx + 1, TCStartIdx + int(elementsCountAroundTC * 0.5))) + list(range(idxTCMid, idxTCMid + int(elementsCountAroundTC * 0.5))) - - for n in range(elementsCountAroundTC - 1): - xTCInner = xList[xTCInnerSet[n]] - xTCOuter = sx[n + 1] - d3 = [xTCOuter[i] - xTCInner[i] for i in range(3)] - d3TCRaw.append(d3) - node = nodes.findNodeByIdentifier(xTCInnerSet[n]+1) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, d3) - - innerIdx = xTCInnerSet[n] - elementsCountThroughWall*(elementsCountAlong+1)*elementsCountAround - d2 = d2List[xTCInnerSet[n]] - factor = factorList[innerIdx] - d2Unscaled = [ 1.0/factor*c for c in d2] - curvature = curvatureAlong[innerIdx] - distance = vector.magnitude([xTCOuter[i] - sxCentralLine[n2][i] for i in range(3)]) - newFactor = 1.0 + curvature*distance - dx_ds2 = [ newFactor*c for c in d2Unscaled] - d2TCRaw.append(dx_ds2) - - xTCOuterList = xTCOuterList + xTCRaw[int((elementsCountAroundTC-2)*0.5):] + xTCRaw[:int((elementsCountAroundTC-2)*0.5)] - d1TCOuterList = d1TCOuterList + d1TCRaw[int((elementsCountAroundTC-2)*0.5):] + d1TCRaw[:int((elementsCountAroundTC-2)*0.5)] - d2TCOuterList = d2TCOuterList + d2TCRaw[int((elementsCountAroundTC-2)*0.5):] + d2TCRaw[:int((elementsCountAroundTC-2)*0.5)] - d3TCOuterList = d3TCOuterList + d3TCRaw[int((elementsCountAroundTC-2)*0.5):] + d3TCRaw[:int((elementsCountAroundTC-2)*0.5)] - - for n in range(len(xTCOuterList)): - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xTCOuterList[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1TCOuterList[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2TCOuterList[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, d3TCOuterList[n]) - if useCrossDerivatives: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) - nodeIdentifier = nodeIdentifier + 1 - - # create elements - if tcCount == 3: - tlMeshGroup = tlGroup.getMeshGroup(mesh) - tmMeshGroup = tmGroup.getMeshGroup(mesh) - toMeshGroup = toGroup.getMeshGroup(mesh) - - for N in range(tcCount + 1): - if N == 0: - for e1 in range(int(elementsCountAroundTC*0.5)): - bni = e1 + 1 - bniTC = e1 + 1 - bniList.append(bni) - bniTCList.append(bniTC) - elif N > 0 and N < tcCount: - for e1 in range(elementsCountAroundTC): - bni = e1 + int(elementsCountAroundTC*0.5) + (N-1)*(elementsCountAroundTC+1) + 2 - bniList.append(bni) - bniTCList.append(bniTC + 1) - for e1 in range(elementsCountAroundTC - 1): - bniTC = bniTC + 1 - bniTCList.append(bniTC) - else: - for e1 in range(int(elementsCountAroundTC*0.5)): - bni = e1 + int(elementsCountAroundTC*0.5) + (N-1)*(elementsCountAroundTC+1) + 2 - bniList.append(bni) - if elementsCountAroundTC > 2: - bniTCList.append(bniTC + 1) - for e1 in range(int(elementsCountAroundTC*0.5 - 1)): - bniTC = bniTC + 1 - bniTCList.append(bniTC) - else: - bniTCList.append(1) - - for e1 in range((elementsCountAroundTC+1)*tcCount): - idxTC = elementsCountThroughWall*(elementsCountAlong+1)*elementsCountAround + setTCIdx[e1] +1 - idxList.append(idxTC) - - for e2 in range(elementsCountAlong): - e1 = -1 - A = e2*(elementsCountAroundTC-1)*tcCount + prevNodeIdentifier - for n1 in range(int(elementsCountAroundTC*0.5)): - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - nodeIdentifiers = regNodeIdentifiers if n1 < int(elementsCountAroundTC*0.5) - 1 else [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.createElement(elementIdentifier, elementtemplate if n1 < int(elementsCountAroundTC*0.5) - 1 else elementtemplate1) - result = element.setNodesByIdentifier(eft if n1 < int(elementsCountAroundTC*0.5) - 1 else eft1, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - toMeshGroup.addElement(element) - - for N in range(tcCount - 1): - for n1 in range(elementsCountAroundTC): - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - if n1 == 0: - nodeIdentifiers = [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.createElement(elementIdentifier, elementtemplate2) - result = element.setNodesByIdentifier(eft2, nodeIdentifiers) - elif n1 > 0 and n1 < elementsCountAroundTC - 1: - nodeIdentifiers = regNodeIdentifiers - element = mesh.createElement(elementIdentifier, elementtemplate) - result = element.setNodesByIdentifier(eft, nodeIdentifiers) - else: - nodeIdentifiers = [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.createElement(elementIdentifier, elementtemplate1) - result = element.setNodesByIdentifier(eft1, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - tlMeshGroup.addElement(element) if N == 0 else tmMeshGroup.addElement(element) - - for n1 in range(int(elementsCountAroundTC*0.5)): - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - nodeIdentifiers = regNodeIdentifiers if n1 > 0 else [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.createElement(elementIdentifier, elementtemplate if n1 > 0 else elementtemplate2) - result = element.setNodesByIdentifier(eft if n1 > 0 else eft2, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - if tcCount == 3: - toMeshGroup.addElement(element) - - # Define texture coordinates field for tenia coli - textureCoordinates = zinc_utils.getOrCreateTextureCoordinateField(fm) - textureNodetemplate1 = nodes.createNodetemplate() - textureNodetemplate1.defineField(textureCoordinates) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - - textureNodetemplate2 = nodes.createNodetemplate() - textureNodetemplate2.defineField(textureCoordinates) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) - if useCrossDerivatives: - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) - - bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) - eftTexture1 = bicubichermitelinear.createEftBasic() - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate.defineField(textureCoordinates, -1, eftTexture1) - - eftTexture2 = bicubichermitelinear.createEftWedgeXi1One() - elementtemplate1.defineField(textureCoordinates, -1, eftTexture2) - - eftTexture3 = bicubichermitelinear.createEftWedgeXi1Zero() - elementtemplate2.defineField(textureCoordinates, -1, eftTexture3) - - eftTexture4 = bicubichermitelinear.createEftOpenTube() - elementtemplate3 = mesh.createElementtemplate() - elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate3.defineField(textureCoordinates, -1, eftTexture4) - - eftTexture5 = bicubichermitelinear.createEftWedgeXi1ZeroOpenTube() - elementtemplate4 = mesh.createElementtemplate() - elementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate4.defineField(textureCoordinates, -1, eftTexture5) - - # Calculate texture coordinates and derivatives - nodeIdentifier = prevNodeIdentifier + 1 - - for N in range(tcCount + 1): # 4 - idxTCMid = N*(elementsCountAroundTC + elementsCountAroundHaustrum) - TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) - TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) - dTC = uList[idxTCMid] - uList[TCStartIdx] if N > 0 else uList[TCEndIdx] - uList[idxTCMid] - v1 = [uList[TCStartIdx], 0.0, 1.0] if N > 0 else [-dTC, 0.0, 1.0] - v2 = [ uList[idxTCMid], 0.0, 2.0] - v3 = [ uList[TCEndIdx], 0.0, 1.0] if N < tcCount else [ 1 + dTC , 0.0, 1.0] - d1 = d2 = d3 = [dTC, 0.0, 0.0] - nx = [v1, v2, v3] - nd1 = [d1, d2, d3] - sx, sd1, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) - if N > 0 and N < tcCount: - w = sx[1:-1] - dw = sd1[1:-1] if elementsCountAroundTC > 2 else nd1[1:-1] - elif N == 0: - w = sx[int(elementsCountAroundTC*0.5):-1] - dw = sd1[int(elementsCountAroundTC*0.5):-1] if elementsCountAroundTC > 2 else nd1[int(elementsCountAroundTC*0.5):-1] - else: - w = sx[1:int(elementsCountAroundTC*0.5)+1] - dw = sd1[1:int(elementsCountAroundTC*0.5)+1] if elementsCountAroundTC > 2 else nd1[1:int(elementsCountAroundTC*0.5)+1] - wList = wList + w - dwList = dwList + dw - - d2 = [0.0, 1.0 / elementsCountAlong, 0.0] - for n2 in range(elementsCountAlong + 1): - for n1 in range((elementsCountAroundTC-1) * tcCount): - u = [ wList[n1][0], - 1.0 / elementsCountAlong * n2, - wList[n1][2]] - node = nodes.findNodeByIdentifier(nodeIdentifier) - node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) - cache.setNode(node) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, u) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dwList[n1]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) - if useCrossDerivatives: - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - if n1 == 0: - u = [ wList[-1][0], - 1.0 / elementsCountAlong * n2, - wList[-1][2]] - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, u) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, dwList[n1]) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2) - if useCrossDerivatives: - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - nodeIdentifier = nodeIdentifier + 1 - - # Define flat coordinates field for tenia coli - flatCoordinates = zinc_utils.getOrCreateFlatCoordinateField(fm) - 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) - - flatElementtemplate = mesh.createElementtemplate() - flatElementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate.defineField(flatCoordinates, -1, eftTexture1) - - flatElementtemplate1 = mesh.createElementtemplate() - flatElementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate1.defineField(flatCoordinates, -1, eftTexture2) - - flatElementtemplate2 = mesh.createElementtemplate() - flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture3) - - flatElementtemplate3 = mesh.createElementtemplate() - flatElementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate3.defineField(flatCoordinates, -1, eftTexture4) - - flatElementtemplate4 = mesh.createElementtemplate() - flatElementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) - flatElementtemplate4.defineField(flatCoordinates, -1, eftTexture5) - - # Calculate flat coordinates and derivatives - nodeIdentifier = prevNodeIdentifier + 1 - - wList = [] - dwList = [] - factor = 3.0 if tcCount == 3 else 2.0 - - for N in range(tcCount + 1): - idxTCMid = N*(elementsCountAroundTC + elementsCountAroundHaustrum) - TCStartIdx = idxTCMid - int(elementsCountAroundTC*0.5) - TCEndIdx = idxTCMid + int(elementsCountAroundTC*0.5) - dTC = (uList[idxTCMid] - uList[TCStartIdx])*arcLengthOuterMidLength if N > 0 else (uList[TCEndIdx] - uList[idxTCMid])*arcLengthOuterMidLength - v1 = [uList[TCStartIdx]*arcLengthOuterMidLength, 0.0, wallThickness] if N > 0 else [-dTC, 0.0, wallThickness] - v2 = [ uList[idxTCMid]*arcLengthOuterMidLength, 0.0, wallThickness + tcThickness] - v3 = [ uList[TCEndIdx]*arcLengthOuterMidLength, 0.0, wallThickness] if N < tcCount else [ arcLengthOuterMidLength + dTC, 0.0, wallThickness] - d1 = d3 = [dTC, 0.0, 0.0] - d2 = [c*factor for c in [dTC, 0.0, 0.0]] - nx = [v1, v2, v3] - nd1 = [d1, d2, d3] - sx, sd1, _, _, _ = interp.sampleCubicHermiteCurves(nx, nd1, elementsCountAroundTC) - if N > 0 and N < tcCount: - w = sx[1:-1] - dw = sd1[1:-1] if elementsCountAroundTC > 2 else nd1[1:-1] - elif N == 0: - w = sx[int(elementsCountAroundTC*0.5):-1] - dw = sd1[int(elementsCountAroundTC*0.5):-1] if elementsCountAroundTC > 2 else nd1[int(elementsCountAroundTC*0.5):-1] - else: - w = sx[1:int(elementsCountAroundTC*0.5)+1] - dw = sd1[1:int(elementsCountAroundTC*0.5)+1] if elementsCountAroundTC > 2 else nd1[1:int(elementsCountAroundTC*0.5)+1] - wList = wList + w - dwList = dwList + dw - - dxds2 = [0.0, totalLengthAlong / elementsCountAlong, 0.0] - - for n2 in range(elementsCountAlong + 1): - for n1 in range((elementsCountAroundTC-1) * tcCount): - x = [ wList[n1][0], - totalLengthAlong / elementsCountAlong * n2, - wList[n1][2]] - node = nodes.findNodeByIdentifier(nodeIdentifier) - node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) - cache.setNode(node) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dwList[n1]) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dxds2) - if useCrossDerivatives: - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - if n1 == 0: - x = [ wList[-1][0], - totalLengthAlong / elementsCountAlong * n2, - wList[-1][2]] - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, x) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, dwList[n1]) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, dxds2) - if useCrossDerivatives: - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - nodeIdentifier = nodeIdentifier + 1 - - # Define texture coordinates field & flat coordinates field over elements - elementIdentifier = prevElementIdentifier + 1 - for e2 in range(elementsCountAlong): - e1 = -1 - A = e2*(elementsCountAroundTC-1)*tcCount + prevNodeIdentifier - - for n1 in range(int(elementsCountAroundTC*0.5)): - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - nodeIdentifiers = regNodeIdentifiers if n1 < int(elementsCountAroundTC*0.5) - 1 else [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.findElementByIdentifier(elementIdentifier) - element.merge(elementtemplate if n1 < int(elementsCountAroundTC*0.5) - 1 else elementtemplate1) - element.merge(flatElementtemplate if n1 < int(elementsCountAroundTC*0.5) - 1 else flatElementtemplate1) - element.setNodesByIdentifier(eftTexture1 if n1 < int(elementsCountAroundTC*0.5) - 1 else eftTexture2, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - - for N in range(tcCount -1): - for n1 in range(elementsCountAroundTC): - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - element = mesh.findElementByIdentifier(elementIdentifier) - if n1 == 0: - nodeIdentifiers = [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element.merge(elementtemplate2) - element.merge(flatElementtemplate2) - element.setNodesByIdentifier(eftTexture3, nodeIdentifiers) - elif n1 > 0 and n1 < elementsCountAroundTC - 1: - nodeIdentifiers = regNodeIdentifiers - element.merge(elementtemplate) - element.merge(flatElementtemplate) - element.setNodesByIdentifier(eftTexture1, nodeIdentifiers) - else: - nodeIdentifiers = [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element.merge(elementtemplate1) - element.merge(flatElementtemplate1) - element.setNodesByIdentifier(eftTexture2, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - - for n1 in range(int(elementsCountAroundTC*0.5)): - onOpening = (n1 == int(elementsCountAroundTC*0.5 - 1)) - e1 = e1 + 1 - regNodeIdentifiers = getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount) - nodeIdentifiers = regNodeIdentifiers if n1 > 0 else [regNodeIdentifiers[i] for i in range(4)] + [regNodeIdentifiers[4], regNodeIdentifiers[6]] - element = mesh.findElementByIdentifier(elementIdentifier) - if n1 > 0: - element.merge(elementtemplate3 if onOpening else elementtemplate) - element.merge(flatElementtemplate3 if onOpening else flatElementtemplate) - element.setNodesByIdentifier(eftTexture4 if onOpening else eftTexture1, nodeIdentifiers) - else: - element.merge(elementtemplate4 if onOpening else elementtemplate2) - element.merge(flatElementtemplate4 if onOpening else flatElementtemplate2) - result = element.setNodesByIdentifier(eftTexture5 if onOpening else eftTexture3, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - - fm.endChange() - - return annotationGroups, nodeIdentifier, elementIdentifier - -def getNodeIdentifierForRegularElement(e1, e2, A, elementsCountAroundTC, elementsCountAround, bniList, bniTCList, idxList, tcCount): - """ - Get node identifiers to create regular elements for tenia coli. - :param e1: Element count iterator around tenia coli - :param e2: Element count iterator along tenia coli - :param A: Element offset - :param elementsCountAroundTC: Number of elements around tenia coli - :param elementsCountAround: Number of elements around scaffold - :param bniList: Base node indices for nodes lying on the outer - surface of haustra over the boundary of tenia coli - :param bniTCList: Base node index for tenia coli nodes - :param idxList: List of global numbering of nodes lying on the outer - surface of haustra over the boundary of tenia coli - :param tcCount: Number of tenia coli - :return: nodeIdentifiers - """ - - bni111 = idxList[bniList[e1]-1] + e2*elementsCountAround - bni121 = idxList[bniList[e1]%((elementsCountAroundTC+1)*tcCount)] + e2*elementsCountAround - bni211 = bni111 + elementsCountAround - bni221 = bni121 + elementsCountAround - bni112 = bniTCList[e1] + A - bni122 = (bniTCList[e1]+1)%((elementsCountAroundTC - 1)*tcCount) + A if (bniTCList[e1]+1)%((elementsCountAroundTC - 1)*tcCount)>0 else bniTCList[e1]+1 + A - bni212 = bni112 + (elementsCountAroundTC-1)*tcCount - bni222 = bni122 + (elementsCountAroundTC-1)*tcCount - nodeIdentifiers = [ bni111, bni121, bni211, bni221, bni112, bni122, bni212, bni222] - - return nodeIdentifiers diff --git a/scaffoldmaker/scaffolds.py b/scaffoldmaker/scaffolds.py index 7739ade0..852721ce 100644 --- a/scaffoldmaker/scaffolds.py +++ b/scaffoldmaker/scaffolds.py @@ -11,8 +11,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1 from scaffoldmaker.meshtypes.meshtype_3d_boxhole1 import MeshType_3d_boxhole1 from scaffoldmaker.meshtypes.meshtype_3d_colon1 import MeshType_3d_colon1 -from scaffoldmaker.meshtypes.meshtype_3d_colonsegmentsimplemesentery1 import MeshType_3d_colonsegmentsimplemesentery1 -from scaffoldmaker.meshtypes.meshtype_3d_colonsegmentteniacoli1 import MeshType_3d_colonsegmentteniacoli1 +from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1 from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1 from scaffoldmaker.meshtypes.meshtype_3d_heart2 import MeshType_3d_heart2 from scaffoldmaker.meshtypes.meshtype_3d_heartarterialroot1 import MeshType_3d_heartarterialroot1 @@ -45,8 +44,7 @@ def __init__(self): MeshType_3d_box1, MeshType_3d_boxhole1, MeshType_3d_colon1, - MeshType_3d_colonsegmentsimplemesentery1, - MeshType_3d_colonsegmentteniacoli1, + MeshType_3d_colonsegment1, MeshType_3d_heart1, MeshType_3d_heart2, MeshType_3d_heartarterialroot1, diff --git a/scaffoldmaker/utils/interpolation.py b/scaffoldmaker/utils/interpolation.py index 24086e05..d3c38a05 100644 --- a/scaffoldmaker/utils/interpolation.py +++ b/scaffoldmaker/utils/interpolation.py @@ -787,3 +787,78 @@ def getDoubleCubicHermiteCurvesMidDerivative(ax, ad1, mx, bx, bd1): magb = vector.magnitude(bd1) magm = arcLengtha + arcLengthb - 0.5*(maga + magb) return vector.setMagnitude(md1, magm) + +def sampleParameterAlongLine(lengthList, paramList, elementsCountOut): + """ + Smooths derivative of parameter with linearly varying lengths, and + samples smoothed parameter at equal distances along the lengths + without including parameter as a component of coordinates. + :param lengthList: List of length locations along a line. + :param paramList: List of parameter values at length locations specified + in lengthList. + :param elementsCountOut: Number of output elements along length. + :return sP, sdP: Parameter values and rate of change at each sampled point. + """ + assert len(lengthList) == len(paramList), 'sampleParameterAlongLine. Mismatched number of lengths and parameters' + nodesCount = len(lengthList) + + md1 = [] + sP = [] + sdP = [] + + # Find smoothed parameter derivatives + # Middle + for n in range(1, nodesCount - 1): + # get mean of directions from point n to points (n - 1) and (n + 1) + nm = n - 1 + np = n + 1 + dirm = paramList[n ] - paramList[nm] + dirp = paramList[np] - paramList[n ] + # mean weighted by fraction towards that end, equivalent to harmonic mean + arcLengthm = lengthList[n ] - lengthList[nm] + arcLengthn = lengthList[np] - lengthList[n ] + arcLengthmp = arcLengthm + arcLengthn + wm = arcLengthn/arcLengthmp + wp = arcLengthm/arcLengthmp + md1.append(wm*dirm + wp*dirp) + + # Start + md1Start = interpolateLagrangeHermiteDerivative([paramList[0]], [paramList[1]], [md1[0]], 0.0) + + # End + md1End = interpolateHermiteLagrangeDerivative([paramList[-2]], [md1[-1]], [paramList[-1]], 1.0) + md1All = md1Start + md1 + md1End + + # Sample into equally spaced elements along line + distance = 0.0 + e = 0 + totalLength = lengthList[-1] + lengthPerElementOut = totalLength / elementsCountOut + dLength = [] + + for n in range(1, nodesCount): + dLength.append(lengthList[n] - lengthList[n-1]) + dLength.append(dLength[-1]) + + for eOut in range(elementsCountOut): + while e < nodesCount - 1: + if distance < lengthList[e + 1]: + partDistance = distance - lengthList[e] + xi = partDistance/(lengthList[e+1] - lengthList[e]) + p = interpolateCubicHermite([paramList[e]], [md1All[e]], [paramList[e+1]], [md1All[e+1]], xi)[0] + dpdxi = interpolateCubicHermiteDerivative([paramList[e]], [md1All[e]], [paramList[e+1]], [md1All[e+1]], xi)[0] + dxdxi = interpolateCubicHermiteDerivative([lengthList[e]], [dLength[e]], [lengthList[e+1]], [dLength[e+1]], xi)[0] + dpdx = dpdxi*1.0/dxdxi + dp = dpdx*lengthPerElementOut + sP.append(p) + sdP.append(dp) + break + e += 1 + distance += lengthPerElementOut + + # Last node + sP.append(paramList[-1]) + dpdx = md1All[-1] * 1.0/dLength[-1] + sdP.append(dpdx*lengthPerElementOut) + + return sP, sdP diff --git a/scaffoldmaker/utils/matrix.py b/scaffoldmaker/utils/matrix.py index fdb88ad5..23e09d45 100644 --- a/scaffoldmaker/utils/matrix.py +++ b/scaffoldmaker/utils/matrix.py @@ -19,3 +19,19 @@ def getRotationMatrixFromAxisAngle(rotAxis, theta): [rotAxis[2]*rotAxis[0]*C - rotAxis[1]*sinTheta, rotAxis[2]*rotAxis[1]*C + rotAxis[0]*sinTheta, rotAxis[2]*rotAxis[2]*C + cosTheta]]) return rotMatrix + +def rotateAboutZAxis(x, theta): + """ + Rotates matrix about z-axis. + : param x: matrix to be rotated. + : param theta: angle of rotation. + :return rotated matrix + """ + cosTheta = math.cos(theta) + sinTheta = math.sin(theta) + + xRot = [x[0]*cosTheta - x[1]*sinTheta, + x[0]*sinTheta + x[1]*cosTheta, + x[2]] + + return xRot diff --git a/scaffoldmaker/utils/tubemesh.py b/scaffoldmaker/utils/tubemesh.py index 641942aa..acacb426 100644 --- a/scaffoldmaker/utils/tubemesh.py +++ b/scaffoldmaker/utils/tubemesh.py @@ -14,67 +14,357 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -def generatetubemesh(region, - elementsCountAround, - elementsCountAlongSegment, - elementsCountThroughWall, - segmentCountAlong, - cx, cd1, cd2, cd12, - xInner, d1Inner, d2Inner, wallThickness, - segmentAxis, segmentLength, - useCrossDerivatives, - useCubicHermiteThroughWall, # or Zinc Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE etc. - annotationGroups, annotationArray, transitElementList, uList, arcLengthOuterMidLength, - firstNodeIdentifier = 1, firstElementIdentifier = 1 - ): - ''' - Generates a 3-D tubular mesh with variable numbers of elements - around, along the central axis, and radially through wall. The - tubular mesh is created from a segment profile which is mapped onto - the central line and lateral axes data - :param elementsCountAround: number of elements around tube - :param elementsCountAlongSegment: number of elements along segment profile - :param elementsCountThroughWall: number of elements through wall thickness - :param segmentCountAlong: number of segments along the tube - :param cx: coordinates on central line - :param cd1: derivative along central line - :param cd2: derivative representing cross axis - :param cd12: rate of change of cd2 along cd1 - :param xInner: coordinates on inner surface of segment profile - :param d1Inner: derivatives around inner surface of segment profile - :param d2Inner: derivatives along inner surface of segment profile - :param wallThickness: thickness of wall - :param segmentAxis: axis of segment profile - :param segmentLength: length of segment profile + +def warpSegmentPoints(xList, d1List, d2List, segmentAxis, segmentLength, + sx, sd1, sd2, elementsCountAround, elementsCountAlongSegment, nSegment): + """ + Warps points in segment to account for bending and twisting + along central path defined by nodes sx and derivatives sd1 and sd2. + :param xList: coordinates of segment points. + :param d1List: derivatives around axis of segment. + :param d2List: derivatives along axis of segment. + :param segmentAxis: axis perpendicular to segment plane. + :param sx: coordinates of points on central path. + :param sd1: derivatives of points along central path. + :param sd2: derivatives representing cross axes. + :param elementsCountAround: Number of elements around segment. + :param elementsCountAlong: Number of elements along segment. + :param nSegment: Segment index along central path. + :return coordinates and derivatives of warped points. + """ + xWarpedList = [] + d1WarpedList = [] + d2WarpedList = [] + smoothd2WarpedList = [] + d3WarpedUnitList = [] + + for nAlongSegment in range(elementsCountAlongSegment + 1): + n2 = elementsCountAlongSegment * nSegment + nAlongSegment + xElementAlongSegment = xList[elementsCountAround*nAlongSegment: elementsCountAround*(nAlongSegment+1)] + d1ElementAlongSegment = d1List[elementsCountAround*nAlongSegment: elementsCountAround*(nAlongSegment+1)] + d2ElementAlongSegment = d2List[elementsCountAround*nAlongSegment: elementsCountAround*(nAlongSegment+1)] + xMid = [0.0, 0.0, segmentLength/elementsCountAlongSegment* nAlongSegment] + + # Rotate to align segment axis with tangent of central line + unitTangent = vector.normalise(sd1[n2]) + cp = vector.crossproduct3(segmentAxis, unitTangent) + dp = vector.dotproduct(segmentAxis, unitTangent) + if vector.magnitude(cp)> 0.0: # path tangent not parallel to segment axis + axisRot = vector.normalise(cp) + thetaRot = math.acos(vector.dotproduct(segmentAxis, unitTangent)) + rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, thetaRot) + midRot = [rotFrame[j][0]*xMid[0] + rotFrame[j][1]*xMid[1] + rotFrame[j][2]*xMid[2] for j in range(3)] + else: # path tangent parallel to segment axis (z-axis) + if dp == -1.0: # path tangent opposite direction to segment axis + thetaRot = math.pi + axisRot = [1.0, 0, 0] + rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, thetaRot) + midRot = [rotFrame[j][0]*xMid[0] + rotFrame[j][1]*xMid[1] + rotFrame[j][2]*xMid[2] for j in range(3)] + else: # segment axis in same direction as unit tangent + midRot = xMid + translateMatrix = [sx[n2][j] - midRot[j] for j in range(3)] + + for n1 in range(elementsCountAround): + x = xElementAlongSegment[n1] + d1 = d1ElementAlongSegment[n1] + d2 = d2ElementAlongSegment[n1] + + if vector.magnitude(cp)> 0.0: # path tangent not parallel to segment axis + xRot1 = [rotFrame[j][0]*x[0] + rotFrame[j][1]*x[1] + rotFrame[j][2]*x[2] for j in range(3)] + d1Rot1 = [rotFrame[j][0]*d1[0] + rotFrame[j][1]*d1[1] + rotFrame[j][2]*d1[2] for j in range(3)] + d2Rot1 = [rotFrame[j][0]*d2[0] + rotFrame[j][1]*d2[1] + rotFrame[j][2]*d2[2] for j in range(3)] + + if n1 == 0: + # Project sd2 onto plane normal to sd1 + v = sd2[n2] + pt = [midRot[j] + sd2[n2][j] for j in range(3)] + dist = vector.dotproduct(v, unitTangent) + ptOnPlane = [pt[j] - dist*unitTangent[j] for j in range(3)] + newVector = [ptOnPlane[j] - midRot[j] for j in range(3)] + # Rotate first point to align with planar projection of sd2 + firstVector = vector.normalise([xRot1[j] - midRot[j] for j in range(3)]) + thetaRot2 = math.acos(vector.dotproduct(vector.normalise(newVector), firstVector)) + cp2 = vector.crossproduct3(vector.normalise(newVector), firstVector) + if vector.magnitude(cp2) > 0.0: + cp2 = vector.normalise(cp2) + signThetaRot2 = vector.dotproduct(unitTangent, cp2) + axisRot2 = unitTangent + rotFrame2 = matrix.getRotationMatrixFromAxisAngle(axisRot2, -signThetaRot2*thetaRot2) + else: + rotFrame2 = [ [1, 0, 0], [0, 1, 0], [0, 0, 1]] + + else: # path tangent parallel to segment axis + xRot1 = [rotFrame[j][0]*x[0] + rotFrame[j][1]*x[1] + rotFrame[j][2]*x[2] for j in range(3)] if dp == -1.0 else x + d1Rot1 = [rotFrame[j][0]*d1[0] + rotFrame[j][1]*d1[1] + rotFrame[j][2]*d1[2] for j in range(3)] if dp == -1.0 else d1 + d2Rot1 = [rotFrame[j][0]*d2[0] + rotFrame[j][1]*d2[1] + rotFrame[j][2]*d2[2] for j in range(3)] if dp == -1.0 else d2 + + # Rotate to align start of elementsAround with sd2 + if n1 == 0: + v = vector.normalise(sd2[n2]) + startVector = vector.normalise([xRot1[j] - midRot[j] for j in range(3)]) + axisRot2 = unitTangent + thetaRot2 = dp*-math.acos(vector.dotproduct(v, startVector)) + rotFrame2 = matrix.getRotationMatrixFromAxisAngle(axisRot2, thetaRot2) + + 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)] + xTranslate = [xRot2[j] + translateMatrix[j] for j in range(3)] + + xWarpedList.append(xTranslate) + d1WarpedList.append(d1Rot2) + d2WarpedList.append(d2Rot2) + + # Smooth d2 for segment + smoothd2Raw = [] + for n1 in range(elementsCountAround): + nx = [] + nd2 = [] + for n2 in range(elementsCountAlongSegment + 1): + n = n2*elementsCountAround + n1 + nx.append(xWarpedList[n]) + nd2.append(d2WarpedList[n]) + smoothd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixStartDerivative = True, fixEndDerivative = True) + smoothd2Raw.append(smoothd2) + + # Re-arrange smoothd2 + for n2 in range(elementsCountAlongSegment + 1): + for n1 in range(elementsCountAround): + smoothd2WarpedList.append(smoothd2Raw[n1][n2]) + + # Calculate unit d3 + for n in range(len(xWarpedList)): + d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1WarpedList[n]), vector.normalise(smoothd2WarpedList[n]))) + d3WarpedUnitList.append(d3Unit) + + return xWarpedList, d1WarpedList, smoothd2WarpedList, d3WarpedUnitList + + +def getCoordinatesFromInner(xInner, d1Inner, d2Inner, d3Inner, + sx, wallThicknessList, elementsCountAround, + elementsCountAlong, elementsCountThroughWall, transitElementList): + """ + Generates coordinates from inner to outer surface using coordinates + and derivatives of inner surface. + :param xInner: Coordinates on inner surface + :param d1Inner: Derivatives on inner surface around tube + :param d2Inner: Derivatives on inner surface along tube + :param d3Inner: Derivatives on inner surface through wall + :param sx: coordinates of points on central path. + :param wallThicknessList: Wall thickness for each element along tube + :param elementsCountAround: Number of elements around tube + :param elementsCountAlong: Number of elements along tube + :param elementsCountThroughWall: Number of elements through tube wall + :param transitElementList: stores true if element around is a transition + element that is between a big and a small element. + return nodes and derivatives for mesh, and curvature along inner surface. + """ + + xOuter = [] + curvatureAroundInner = [] + curvatureAlong = [] + curvatureList = [] + xList = [] + d1List = [] + d2List = [] + d3List = [] + + for n2 in range(elementsCountAlong + 1): + wallThickness = wallThicknessList[n2] + for n1 in range(elementsCountAround): + n = n2*elementsCountAround + n1 + norm = d3Inner[n] + # Calculate outer coordinates + x = [xInner[n][i] + norm[i]*wallThickness for i in range(3)] + xOuter.append(x) + # Calculate curvature along elements around + prevIdx = n - 1 if (n1 != 0) else (n2 + 1)*elementsCountAround - 1 + nextIdx = n + 1 if (n1 < elementsCountAround - 1) else n2*elementsCountAround + kappam = interp.getCubicHermiteCurvatureSimple(xInner[prevIdx], d1Inner[prevIdx], xInner[n], d1Inner[n], 1.0) + kappap = interp.getCubicHermiteCurvatureSimple(xInner[n], d1Inner[n], xInner[nextIdx], d1Inner[nextIdx], 0.0) + if not transitElementList[n1] and not transitElementList[(n1-1)%elementsCountAround]: + curvatureAround = 0.5*(kappam + kappap) + elif transitElementList[n1]: + curvatureAround = kappam + elif transitElementList[(n1-1)%elementsCountAround]: + curvatureAround = kappap + curvatureAroundInner.append(curvatureAround) + + # Calculate curvature along + if n2 == 0: + curvature = abs(interp.getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[n + elementsCountAround], d2Inner[n + elementsCountAround], vector.normalise(d3Inner[n]), 0.0)) + elif n2 == elementsCountAlong: + curvature = abs(interp.getCubicHermiteCurvature(xInner[n - elementsCountAround], d2Inner[n - elementsCountAround], xInner[n], d2Inner[n], vector.normalise(d3Inner[n]), 1.0)) + else: + curvature = 0.5*( + abs(interp.getCubicHermiteCurvature(xInner[n - elementsCountAround], d2Inner[n - elementsCountAround], xInner[n], d2Inner[n], vector.normalise(d3Inner[n]), 1.0)) + + abs(interp.getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[n + elementsCountAround], d2Inner[n + elementsCountAround], vector.normalise(d3Inner[n]), 0.0))) + curvatureAlong.append(curvature) + + for n3 in range(elementsCountThroughWall + 1): + xi3 = 1/elementsCountThroughWall * n3 + for n1 in range(elementsCountAround): + n = n2*elementsCountAround + n1 + norm = d3Inner[n] + innerx = xInner[n] + outerx = xOuter[n] + dWall = [wallThickness*c for c in norm] + # x + x = interp.interpolateCubicHermite(innerx, dWall, outerx, dWall, xi3) + xList.append(x) + + # dx_ds1 + factor = 1.0 + wallThickness*xi3 * curvatureAroundInner[n] + d1 = [ factor*c for c in d1Inner[n]] + d1List.append(d1) + + # dx_ds2 + curvature = curvatureAlong[n] + distance = vector.magnitude([x[i] - xInner[n][i] for i in range(3)]) + factor = 1.0 - curvature*distance + d2 = [ factor*c for c in d2Inner[n]] + d2List.append(d2) + curvatureList.append(curvature) + + #dx_ds3 + d3 = [c * wallThickness/elementsCountThroughWall for c in norm] + d3List.append(d3) + + return xList, d1List, d2List, d3List, curvatureList + +def createFlatAndTextureCoordinates(xiList, lengthAroundList, + totalLengthAlong, wallThickness, elementsCountAround, + elementsCountAlong, elementsCountThroughWall, transitElementList): + """ + Calculates flat coordinates and texture coordinates + for a tube when it is opened into a flat preparation. + :param xiList: List containing xi for each point around + the outer surface of the tube. + :param lengthAroundList: List of total arclength around the + outer surface for each element along. + :param totalLengthAlong: Total length along tube. + :param wallThickness: Thickness of wall. + :param elementsCountAround: Number of elements around tube. + :param elementsCountAlong: Number of elements along tube. + :param elementsCountThroughWall: Number of elements through wall. + :param transitElementList: stores true if element around is a + transition element between a big and small element. + :return: coordinates and derivatives of flat and texture coordinates fields. + """ + + # Calculate flat coordinates and derivatives + xFlatList = [] + d1FlatList = [] + d2FlatList = [] + for n2 in range(elementsCountAlong + 1): + xiFace = xiList[n2] + lengthAround = lengthAroundList[n2] + d1List = [] + for n1 in range(len(xiFace)): + d1 = (xiFace[n1] - xiFace[n1-1]) if n1 > 0 else (xiFace[n1+1] - xiFace[n1]) + d1List.append(d1) + + # To modify derivative along transition elements + for i in range(len(transitElementList)): + if transitElementList[i]: + d1List[i+1] = d1List[i+2] + + xPad = (lengthAroundList[0] - lengthAround)*0.5 + for n3 in range(elementsCountThroughWall + 1): + z = wallThickness / elementsCountThroughWall * n3 + for n1 in range(elementsCountAround + 1): + xFlat = [xPad + xiFace[n1] * lengthAround, + totalLengthAlong / elementsCountAlong * n2, + z] + d1Flat = [ d1List[n1]*lengthAround, 0.0, 0.0 ] + xFlatList.append(xFlat) + d1FlatList.append(d1Flat) + + for n2 in range(elementsCountAlong): + for n3 in range(elementsCountThroughWall + 1): + for n1 in range(elementsCountAround + 1 ): + nodeIdx = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + n3*(elementsCountAround + 1) + n1 + nodeNextElementAlong = nodeIdx + (elementsCountAround+1)*(elementsCountThroughWall + 1) + # print(nodeIdx + 1, nodeNextElementAlong + 1) + v1 = xFlatList[nodeNextElementAlong] + v2 = xFlatList[nodeIdx] + d1 = d2 = [v1[i] - v2[i] for i in range(3)] + arclength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, True) + d2Flat = vector.setMagnitude(d1, arclength) + d2FlatList.append(d2Flat) + d2FlatList = d2FlatList + d2FlatList[-(elementsCountAround+1)*(elementsCountThroughWall+1):] + + # Calculate texture coordinates and derivatives + xTextureList = [] + d1Texture = [] + d1TextureList = [] + d2TextureList = [] + + d2 = [0.0, 1.0 / elementsCountAlong, 0.0] + xiTexture = xiList[0] + + for n1 in range(len(xiTexture)): + d1 = [xiTexture[n1] - xiTexture[n1-1] if n1 > 0 else xiTexture[n1+1] - xiTexture[n1], + 0.0, + 0.0] + d1Texture.append(d1) + + # To modify derivative along transition elements + for i in range(len(transitElementList)): + if transitElementList[i]: + d1Texture[i+1] = d1Texture[i+2] + + for n2 in range(elementsCountAlong + 1): + for n3 in range(elementsCountThroughWall + 1): + for n1 in range(elementsCountAround + 1): + u = [ xiTexture[n1], + 1.0 / elementsCountAlong * n2, + 1.0 / elementsCountThroughWall * n3] + d1 = d1Texture[n1] + xTextureList.append(u) + d1TextureList.append(d1) + d2TextureList.append(d2) + + return xFlatList, d1FlatList, d2FlatList, xTextureList, d1TextureList, d2TextureList + +def createNodesAndElements(region, + x, d1, d2, d3, + xFlat, d1Flat, d2Flat, + xTexture, d1Texture, d2Texture, + elementsCountAround, elementsCountAlong, elementsCountThroughWall, + annotationGroups, annotationArray, + firstNodeIdentifier, firstElementIdentifier, + useCubicHermiteThroughWall, useCrossDerivatives): + """ + Create nodes and elements for the coordinates, flat coordinates, + and texture coordinates field. + :param x, d1, d2, d3: coordinates and derivatives of coordinates field. + :param xFlat, d1Flat, d2Flat, d3Flat: coordinates and derivatives of + flat coordinates field. + :param xTexture, d1Texture, d2Texture, d3Texture: coordinates and derivatives + of texture coordinates field. + :param elementsCountAround: Number of elements around tube. + :param elementsCountAlong: Number of elements along tube. + :param elementsCountThroughWall: Number of elements through wall. + :param annotationGroups: stores information about annotation groups. + :param annotationArray: stores annotation names of elements around. + :param firstNodeIdentifier, firstElementIdentifier: first node and + element identifier to use. :param useCubicHermiteThroughWall: use linear when false - :param annotationGroups: Empty when not required - :param annotationArray: Array storing annotation group name for elements around - :param transitElementList: True false list of where transition elements are - located around - :param uList: List of xi for each node around representative cross-sectional - profile - :param arcLengthOuterMidLength: Total arclength of elements around outer surface - along mid-length of a segment - :return nodeIdentifier, elementIdentifier - :return xList, d1List, d2List, d3List: List of coordinates and derivatives - on tube - :return sx: List of coordinates sampled from central line - :return curvatureAlong, factorList: List of curvature and scale factor along mesh - for each node on inner surface of mesh. - ''' - - zero = [0.0, 0.0, 0.0] - elementsCountAlong = elementsCountAlongSegment*segmentCountAlong - - # Sample central line to get same number of elements as elementsCountAlong - sx, sd1, se, sxi, ssf = interp.sampleCubicHermiteCurves(cx, cd1, elementsCountAlong) - sd2, _ = interp.interpolateSampleCubicHermite(cd2, cd12, se, sxi, ssf) + :param useCrossDerivatives: use cross derivatives when true + :return nodeIdentifier, elementIdentifier, annotationGroups + """ + + nodeIdentifier = firstNodeIdentifier + elementIdentifier = firstElementIdentifier + zero = [ 0.0, 0.0, 0.0 ] fm = region.getFieldmodule() fm.beginChange() cache = fm.createFieldcache() - coordinates = zinc_utils.getOrCreateCoordinateField(fm) + # Coordinates field + coordinates = zinc_utils.getOrCreateCoordinateField(fm) nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) nodetemplate = nodes.createNodetemplate() nodetemplate.defineField(coordinates) @@ -102,238 +392,11 @@ def generatetubemesh(region, elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) result = elementtemplate.defineField(coordinates, -1, eft) - # create nodes - nodeIdentifier = firstNodeIdentifier - x = [ 0.0, 0.0, 0.0 ] - dx_ds1 = [ 0.0, 0.0, 0.0 ] - dx_ds2 = [ 0.0, 0.0, 0.0 ] - dx_ds3 = [ 0.0, 0.0, 0.0 ] - xInnerList = [] - d1InnerList = [] - d2InnerList = [] - d3InnerUnitList = [] - xList = [] - dx_ds1List = [] - dx_ds2List = [] - dx_ds3List = [] - curvatureAlong = [] - smoothd2Raw = [] - smoothd2InnerList = [] - d1List = [] - -# Map each face along segment profile to central line - for nSegment in range(segmentCountAlong): - for nAlongSegment in range(elementsCountAlongSegment + 1): - n2 = nSegment*elementsCountAlongSegment + nAlongSegment - if nSegment == 0 or (nSegment > 0 and nAlongSegment > 0): - # Rotate to align segment axis with tangent of central line - segmentMid = [0.0, 0.0, segmentLength/elementsCountAlongSegment* nAlongSegment] - unitTangent = vector.normalise(sd1[n2]) - cp = vector.crossproduct3(segmentAxis, unitTangent) - if vector.magnitude(cp)> 0.0: - axisRot = vector.normalise(cp) - thetaRot = math.acos(vector.dotproduct(segmentAxis, unitTangent)) - rotFrame = matrix.getRotationMatrixFromAxisAngle(axisRot, thetaRot) - midRot = [rotFrame[j][0]*segmentMid[0] + rotFrame[j][1]*segmentMid[1] + rotFrame[j][2]*segmentMid[2] for j in range(3)] - translateMatrix = [sx[n2][j] - midRot[j] for j in range(3)] - else: - midRot = segmentMid - - for n1 in range(elementsCountAround): - n = nAlongSegment*elementsCountAround + n1 - x = xInner[n] - d1 = d1Inner[n] - d2 = d2Inner[n] - if vector.magnitude(cp)> 0.0: - xRot1 = [rotFrame[j][0]*x[0] + rotFrame[j][1]*x[1] + rotFrame[j][2]*x[2] for j in range(3)] - d1Rot1 = [rotFrame[j][0]*d1[0] + rotFrame[j][1]*d1[1] + rotFrame[j][2]*d1[2] for j in range(3)] - d2Rot1 = [rotFrame[j][0]*d2[0] + rotFrame[j][1]*d2[1] + rotFrame[j][2]*d2[2] for j in range(3)] - if n1 == 0: - # Project sd2 onto plane normal to sd1 - v = sd2[n2] - pt = [midRot[j] + sd2[n2][j] for j in range(3)] - dist = vector.dotproduct(v, unitTangent) - ptOnPlane = [pt[j] - dist*unitTangent[j] for j in range(3)] - newVector = [ptOnPlane[j] - midRot[j] for j in range(3)] - # Rotate first point to align with planar projection of sd2 - firstVector = vector.normalise([xRot1[j] - midRot[j] for j in range(3)]) - thetaRot2 = math.acos(vector.dotproduct(vector.normalise(newVector), firstVector)) - cp2 = vector.crossproduct3(vector.normalise(newVector), firstVector) - if vector.magnitude(cp2) > 0.0: - cp2 = vector.normalise(cp2) - signThetaRot2 = vector.dotproduct(unitTangent, cp2) - axisRot2 = unitTangent - rotFrame2 = matrix.getRotationMatrixFromAxisAngle(axisRot2, -signThetaRot2*thetaRot2) - else: - rotFrame2 = [ [1, 0, 0], [0, 1, 0], [0, 0, 1]] - 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)] - else: - xRot2 = x - d1Rot2 = d1 - d2Rot2 = d2 - xTranslate = [xRot2[j] + translateMatrix[j] for j in range(3)] - xInnerList.append(xTranslate) - d1InnerList.append(d1Rot2) - d2InnerList.append(d2Rot2) - d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1Rot2), vector.normalise(d2Rot2))) - d3InnerUnitList.append(d3Unit) - - for n1 in range(elementsCountAround): - nx = [] - nd2 = [] - for n2 in range(elementsCountAlong + 1): - n = n2*elementsCountAround + n1 - nx.append(xInnerList[n]) - nd2.append(d2InnerList[n]) - smoothd2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2) - smoothd2Raw.append(smoothd2) - - for n2 in range(elementsCountAlong + 1): - for n1 in range(elementsCountAround): - smoothd2InnerList.append(smoothd2Raw[n1][n2]) - n = elementsCountAround * n2 + n1 - if n2 == 0: - curvature = abs(interp.getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], vector.normalise(sd2[n2]), 0.0)) - elif n2 == elementsCountAlong: - curvature = abs(interp.getCubicHermiteCurvature(sx[n2-1], sd1[n2-1], sx[n2], sd1[n2], vector.normalise(sd2[n2]), 1.0)) - else: - curvature = 0.5*( - abs(interp.getCubicHermiteCurvature(sx[n2-1], sd1[n2-1], sx[n2], sd1[n2], vector.normalise(sd2[n2]), 1.0)) + - abs(interp.getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], vector.normalise(sd2[n2]), 0.0))) - curvatureAlong.append(curvature) - - # Pre-calculate node locations and derivatives on outer boundary - xOuterList, curvatureInner = getOuterCoordinatesAndCurvatureFromInner(xInnerList, d1InnerList, d3InnerUnitList, wallThickness, elementsCountAlong, elementsCountAround, transitElementList) - - # Interpolate to get nodes through wall - for n3 in range(elementsCountThroughWall + 1): - xi3 = 1/elementsCountThroughWall * n3 - x, dx_ds1, dx_ds2, dx_ds3, factorList = interpolatefromInnerAndOuter( xInnerList, xOuterList, - wallThickness, xi3, sx, curvatureInner, curvatureAlong, d1InnerList, smoothd2InnerList, d3InnerUnitList, - elementsCountAround, elementsCountAlong, elementsCountThroughWall) - xList = xList + x - dx_ds1List = dx_ds1List + dx_ds1 - dx_ds2List = dx_ds2List + dx_ds2 - dx_ds3List = dx_ds3List + dx_ds3 - - for n in range(len(xList)): - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xList[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1List[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2List[n]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3List[n]) - if useCrossDerivatives: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) - # print('NodeIdentifier = ', nodeIdentifier, xList[n]) - nodeIdentifier = nodeIdentifier + 1 - - # # For debugging - Nodes along central line - # for pt in range(len(sx)): - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, sx[pt]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, sd1[pt]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, vector.normalise(sd2[pt])) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, vector.normalise(testRot[pt])) - # nodeIdentifier = nodeIdentifier + 1 - - # create elements - elementIdentifier = firstElementIdentifier - now = (elementsCountAlong + 1)*elementsCountAround - for e3 in range(elementsCountThroughWall): - for e2 in range(elementsCountAlong): - for e1 in range(elementsCountAround): - element = mesh.createElement(elementIdentifier, elementtemplate) - bni11 = e3*now + e2*elementsCountAround + e1 + 1 - bni12 = e3*now + e2*elementsCountAround + (e1 + 1) % elementsCountAround + 1 - bni21 = e3*now + (e2 + 1)*elementsCountAround + e1 + 1 - bni22 = e3*now + (e2 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1 - nodeIdentifiers = [ bni11, bni12, bni21, bni22, bni11 + now, bni12 + now, bni21 + now, bni22 + now ] - result = element.setNodesByIdentifier(eft, nodeIdentifiers) - elementIdentifier = elementIdentifier + 1 - if annotationGroups: - for annotationGroup in annotationGroups: - if annotationArray[e1] == annotationGroup._name: - meshGroup = annotationGroup.getMeshGroup(mesh) - meshGroup.addElement(element) - - # Define texture coordinates field - textureCoordinates = zinc_utils.getOrCreateTextureCoordinateField(fm) - textureNodetemplate1 = nodes.createNodetemplate() - textureNodetemplate1.defineField(textureCoordinates) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - textureNodetemplate1.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - - textureNodetemplate2 = nodes.createNodetemplate() - textureNodetemplate2.defineField(textureCoordinates) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_VALUE, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS1, 2) - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D_DS2, 2) - if useCrossDerivatives: - textureNodetemplate2.setValueNumberOfVersions(textureCoordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 2) - + # Flat coordinates field bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) eftTexture1 = bicubichermitelinear.createEftBasic() - - elementtemplate1 = mesh.createElementtemplate() - elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate1.defineField(textureCoordinates, -1, eftTexture1) - eftTexture2 = bicubichermitelinear.createEftOpenTube() - elementtemplate2 = mesh.createElementtemplate() - elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate2.defineField(textureCoordinates, -1, eftTexture2) - # Calculate texture coordinates and derivatives - d2 = [0.0, 1.0 / elementsCountAlong, 0.0] - - for n1 in range(len(uList)): - d1 = [uList[n1] - uList[n1-1] if n1 > 0 else uList[n1+1] - uList[n1], - 0.0, - 0.0] - d1List.append(d1) - - # To modify derivative along transition elements - for i in range(len(transitElementList)): - if transitElementList[i]: - d1List[i+1] = d1List[i+2] - - nodeIdentifier = firstNodeIdentifier - for n3 in range(elementsCountThroughWall + 1): - for n2 in range(elementsCountAlong + 1): - for n1 in range(elementsCountAround): - u = [ uList[n1], - 1.0 / elementsCountAlong * n2, - 1.0 / elementsCountThroughWall * n3] - d1 = d1List[n1] - node = nodes.findNodeByIdentifier(nodeIdentifier) - node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) - cache.setNode(node) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, u) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) - if useCrossDerivatives: - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - if n1 == 0: - u = [ 1.0, 1.0 / elementsCountAlong * n2, 1.0 / elementsCountThroughWall * n3] - d1 = d1List[-1] - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, u) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1) - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2) - if useCrossDerivatives: - textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) - nodeIdentifier = nodeIdentifier + 1 - - # Define flat coordinates field flatCoordinates = zinc_utils.getOrCreateFlatCoordinateField(fm) flatNodetemplate1 = nodes.createNodetemplate() flatNodetemplate1.defineField(flatCoordinates) @@ -359,157 +422,105 @@ def generatetubemesh(region, flatElementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) flatElementtemplate2.defineField(flatCoordinates, -1, eftTexture2) - # Calculate texture coordinates and derivatives - totalLengthAlong = segmentLength*segmentCountAlong - d2 = [0.0, totalLengthAlong/elementsCountAlong, 0.0] + # Texture coordinates field + textureCoordinates = zinc_utils.getOrCreateTextureCoordinateField(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) - d1List = [] - for n1 in range(len(uList)): - d1 = [(uList[n1] - uList[n1-1])*arcLengthOuterMidLength if n1 > 0 else (uList[n1+1] - uList[n1])*arcLengthOuterMidLength, - 0.0, - 0.0] - d1List.append(d1) + 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) - # To modify derivative along transition elements - for i in range(len(transitElementList)): - if transitElementList[i]: - d1List[i+1] = d1List[i+2] + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate1.defineField(textureCoordinates, -1, eftTexture1) + elementtemplate2 = mesh.createElementtemplate() + elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate2.defineField(textureCoordinates, -1, eftTexture2) + + # Create nodes + # Coordinates field + for n in range(len(x)): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, d3[n]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + # print('NodeIdentifier = ', nodeIdentifier, xList[n]) + nodeIdentifier = nodeIdentifier + 1 + + # Flat and texture coordinates fields nodeIdentifier = firstNodeIdentifier - for n3 in range(elementsCountThroughWall + 1): - for n2 in range(elementsCountAlong + 1): + for n2 in range(elementsCountAlong + 1): + for n3 in range(elementsCountThroughWall + 1): for n1 in range(elementsCountAround): - x = [ uList[n1]*arcLengthOuterMidLength, - totalLengthAlong / elementsCountAlong * n2, - wallThickness / elementsCountThroughWall * n3] - d1 = d1List[n1] + i = n2*(elementsCountAround + 1)*(elementsCountThroughWall + 1) + (elementsCountAround + 1)*n3 + n1 node = nodes.findNodeByIdentifier(nodeIdentifier) node.merge(flatNodetemplate2 if n1 == 0 else flatNodetemplate1) + node.merge(textureNodetemplate2 if n1 == 0 else textureNodetemplate1) cache.setNode(node) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) + # print('NodeIdentifier', nodeIdentifier, 'version 1, xList Index =', i+1) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xFlat[i]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Flat[i]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Flat[i]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xTexture[i]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1Texture[i]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Texture[i]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) if n1 == 0: - x = [ arcLengthOuterMidLength, totalLengthAlong / elementsCountAlong * n2, wallThickness / elementsCountThroughWall * n3] - d1 = d1List[-1] - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, x) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1) - flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2) + # print('NodeIdentifier', nodeIdentifier, 'version 2, xList Index =', i+elementsCountAround+1) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xFlat[i+elementsCountAround]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Flat[i+elementsCountAround]) + flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Flat[i+elementsCountAround]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 2, xTexture[i+elementsCountAround]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 2, d1Texture[i+elementsCountAround]) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 2, d2Texture[i+elementsCountAround]) if useCrossDerivatives: flatCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) + textureCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 2, zero) nodeIdentifier = nodeIdentifier + 1 - # Define flat coordinates field & texture coordinates field over elements - elementIdentifier = firstElementIdentifier - now = (elementsCountAlong + 1)*elementsCountAround - - for e3 in range(elementsCountThroughWall): - for e2 in range(elementsCountAlong): + # create elements + now = elementsCountAround*(elementsCountThroughWall+1) + for e2 in range(elementsCountAlong): + for e3 in range(elementsCountThroughWall): for e1 in range(elementsCountAround): + bni11 = e2*now + e3*elementsCountAround + e1 + 1 + bni12 = e2*now + e3*elementsCountAround + (e1 + 1) % elementsCountAround + 1 + bni21 = e2*now + (e3 + 1)*elementsCountAround + e1 + 1 + bni22 = e2*now + (e3 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1 + nodeIdentifiers = [ bni11, bni12, bni11 + now, bni12 + now, bni21, bni22, bni21 + now, bni22 + now ] onOpening = e1 > elementsCountAround - 2 - element = mesh.findElementByIdentifier(elementIdentifier) - element.merge(elementtemplate2 if onOpening else elementtemplate1) + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nodeIdentifiers) element.merge(flatElementtemplate2 if onOpening else flatElementtemplate1) - bni11 = e3*now + e2*elementsCountAround + e1 + 1 - bni12 = e3*now + e2*elementsCountAround + (e1 + 1) % elementsCountAround + 1 - bni21 = e3*now + (e2 + 1)*elementsCountAround + e1 + 1 - bni22 = e3*now + (e2 + 1)*elementsCountAround + (e1 + 1) % elementsCountAround + 1 - nodeIdentifiers = [ bni11, bni12, bni21, bni22, bni11 + now, bni12 + now, bni21 + now, bni22 + now ] + element.merge(elementtemplate2 if onOpening else elementtemplate1) element.setNodesByIdentifier(eftTexture2 if onOpening else eftTexture1, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 + if annotationGroups: + for annotationGroup in annotationGroups: + if annotationArray[e1] == annotationGroup._name: + meshGroup = annotationGroup.getMeshGroup(mesh) + meshGroup.addElement(element) fm.endChange() - return annotationGroups, nodeIdentifier, elementIdentifier, xList, dx_ds1List, dx_ds2List, dx_ds3List, sx, curvatureAlong, factorList - -def getOuterCoordinatesAndCurvatureFromInner(xInner, d1Inner, d3Inner, wallThickness, elementsCountAlong, elementsCountAround, transitElementList): - """ - Generates coordinates on outer surface and curvature of inner - surface from coordinates and derivatives of inner surface using - wall thickness and normals. - param xInner: Coordinates on inner surface - param d1Inner: Derivatives on inner surface around tube - param d3Inner: Derivatives on inner surface through wall - param wallThickness: Thickness of wall - param elementsCountAlong: Number of elements along tube - param elementsCountAround: Number of elements around tube - return xOuter: Coordinates on outer surface - return curvatureInner: Curvature of coordinates on inner surface - """ - xOuter = [] - curvatureInner = [] - for n2 in range(elementsCountAlong + 1): - for n1 in range(elementsCountAround): - n = n2*elementsCountAround + n1 - x = [xInner[n][i] + d3Inner[n][i]*wallThickness for i in range(3)] - prevIdx = n - 1 if (n1 != 0) else (n2 + 1)*elementsCountAround - 1 - nextIdx = n + 1 if (n1 < elementsCountAround - 1) else n2*elementsCountAround - norm = d3Inner[n] - kappam = interp.getCubicHermiteCurvatureSimple(xInner[prevIdx], d1Inner[prevIdx], xInner[n], d1Inner[n], 1.0) - kappap = interp.getCubicHermiteCurvatureSimple(xInner[n], d1Inner[n], xInner[nextIdx], d1Inner[nextIdx], 0.0) - if not transitElementList[n1] and not transitElementList[(n1-1)%elementsCountAround]: - curvatureAround = 0.5*(kappam + kappap) - elif transitElementList[n1]: - curvatureAround = kappam - elif transitElementList[(n1-1)%elementsCountAround]: - curvatureAround = kappap - xOuter.append(x) - curvatureInner.append(curvatureAround) - - return xOuter, curvatureInner - -def interpolatefromInnerAndOuter( xInner, xOuter, thickness, xi3, sx, curvatureInner, curvatureAlong, - d1Inner, d2Inner, d3InnerUnit, elementsCountAround, elementsCountAlong, elementsCountThroughWall): - """ - Generate coordinates and derivatives at xi3 by interpolating with - inner and outer coordinates and derivatives. - :param xInner: Coordinates on inner surface - :param xOuter: Coordinates on outer surface - :param thickness: Thickness of wall - :param sx: List of coordinates sampled from central line - :param curvatureInner: Curvature of coordinates on inner surface - :param curvatureAlong: Curvature of coordinates on inner surface along mesh - :param d1Inner: Derivatives on inner surface around tube - :param d2Inner: Derivatives on inner surface along tube - :param d3InnerUnit: Unit derivatives on inner surface through wall - :param elementsCountAround: Number of elements around tube - :param elementsCountAlong: Number of elements along tube - :param elementsCountThroughWall: Number of elements through wall - :return xList, dx_ds1List, dx_ds2List, dx_ds3List: Coordinates and derivatives on xi3 - :return factorList: List of factors used for scaling d2 - """ - xList = [] - dx_ds1List = [] - dx_ds2List = [] - dx_ds3List =[] - factorList = [] - - for n2 in range(elementsCountAlong + 1): - for n1 in range(elementsCountAround): - n = n2*elementsCountAround + n1 - norm = d3InnerUnit[n] - # x - innerx = xInner[n] - outerx = xOuter[n] - dWall = [thickness*c for c in norm] - x = interp.interpolateCubicHermite(innerx, dWall, outerx, dWall, xi3) - xList.append(x) - # dx_ds1 - factor = 1.0 + thickness*xi3 * curvatureInner[n] - dx_ds1 = [ factor*c for c in d1Inner[n]] - dx_ds1List.append(dx_ds1) - # dx_ds2 - curvature = curvatureAlong[n] - distance = vector.magnitude([x[i] - sx[n2][i] for i in range(3)]) - factor = 1.0 + curvature*distance - dx_ds2 = [ factor*c for c in d2Inner[n]] - dx_ds2List.append(dx_ds2) - factorList.append(factor) - - #dx_ds3 - dx_ds3 = [c * thickness/elementsCountThroughWall for c in norm] - dx_ds3List.append(dx_ds3) - - return xList, dx_ds1List, dx_ds2List, dx_ds3List, factorList + return nodeIdentifier, elementIdentifier, annotationGroups