diff --git a/src/scaffoldmaker/annotation/body_terms.py b/src/scaffoldmaker/annotation/body_terms.py index 89ae6c7b..cda4a8a8 100644 --- a/src/scaffoldmaker/annotation/body_terms.py +++ b/src/scaffoldmaker/annotation/body_terms.py @@ -17,6 +17,8 @@ ( "core", "" ), ( "non core", "" ), ( "core boundary", "" ), + ( "arm", "UBERON:0001460"), + ( "leg", "UBERON:0000978") ] def get_body_term(name : str): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py index b87b61df..293cb2ad 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py @@ -1,5 +1,5 @@ """ -Generates a 3-D Hermite bifurcating tube network. +Generates a 3-D Hermite bifurcating tube network (with optional solid core). """ from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base @@ -7,9 +7,81 @@ from scaffoldmaker.utils.tubenetworkmesh import TubeNetworkMeshBuilder, TubeNetworkMeshGenerateData +def calculateElementsCountAcrossMinor(cls, options): + # Calculate elementsCountAcrosMinor + elementsCountAcrossMinor = int(((options["Elements count around"] - 4) / 4 - + (options['Number of elements across core major'] / 2)) * 2 + 6) + + return elementsCountAcrossMinor + +def setParametersToDefault(cls, options): + options["Elements count around"] = 8 + options["Number of elements across core major"] = 4 + options["Number of elements across core transition"] = 1 + + return options + +def checkCoreParameters(cls, options, dependentChanges=False): + # Check elements count around are ok + if options["Elements count around"] < 8: + dependentChanges = True + options = setParametersToDefault(cls, options) + if options["Elements count around"] % 4: + options["Elements count around"] += 4 - options["Elements count around"] % 4 + + # Check elements count across major are ok + if options['Number of elements across core major'] < 4: + options['Number of elements across core major'] = 4 + if options['Number of elements across core major'] % 2: + options['Number of elements across core major'] += 1 + + # Check elements count across transition + if options['Number of elements across core transition'] < 1: + options['Number of elements across core transition'] = 1 + + # Calculate elementsCountAcrossMinor based on the current set of elementsCountAround and elementsCountAcrossMajor + elementsCountAcrossMinor = calculateElementsCountAcrossMinor(cls, options) + + # Rcrit check + Rcrit = max( + min(options['Number of elements across core major'] - 4, elementsCountAcrossMinor - 4) // 2 + 1, 0) + if Rcrit < options['Number of elements across core transition']: + dependentChanges = True + options['Number of elements across core transition'] = Rcrit + + # Number of elements around sanity check + eM = options["Number of elements across core major"] + em = elementsCountAcrossMinor + eC = (eM - 1) * (em - 1) - ((eM - 3) * (em - 3)) + if options["Elements count around"] != eC: + dependentChanges = True + # Reset parameter values + setParametersToDefault(cls, options) + + annotationElementsCountsAround = options["Annotation elements counts around"] + if len(annotationElementsCountsAround) == 0: + options["Annotation elements count around"] = [0] + else: + for i in range(len(annotationElementsCountsAround)): + if annotationElementsCountsAround[i] <= 0: + annotationElementsCountsAround[i] = 0 + elif annotationElementsCountsAround[i] < 4: + annotationElementsCountsAround[i] = 4 + + if options["Use linear through wall"]: + dependentChanges = True + options["Use linear through wall"] = False + + return options, dependentChanges + + class MeshType_3d_tubenetwork1(Scaffold_base): """ Generates a 3-D hermite tube network, with linear or cubic basis through wall. + Number of elements generated are primarily controlled by "elements count around" and "elements count across major". + The elements count around parameter determines the number of elements around a circular rim. + The elements count across major determines the number of elements across the major axis of an ellipse (y-axis in the scaffold). + The number of elements across minor axis (z-axis) is calculated internally based on the two elements count parameters. """ @classmethod @@ -31,7 +103,11 @@ def getDefaultOptions(cls, parameterSetName="Default"): "Annotation elements counts around": [0], "Target element density along longest segment": 4.0, "Use linear through wall": False, - "Show trim surfaces": False + "Show trim surfaces": False, + "Core": False, + 'Number of elements across core major': 4, + 'Number of elements across core transition': 1, + "Annotation elements counts across major": [0] } return options @@ -44,7 +120,11 @@ def getOrderedOptionNames(): "Annotation elements counts around", "Target element density along longest segment", "Use linear through wall", - "Show trim surfaces" + "Show trim surfaces", + "Core", + 'Number of elements across core major', + 'Number of elements across core transition', + "Annotation elements counts across major" ] @classmethod @@ -82,23 +162,33 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non def checkOptions(cls, options): if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"): options["Network layout"] = cls.getOptionScaffoldPackage("Network layout", MeshType_1d_network_layout1) - elementsCountAround = options["Elements count around"] - if options["Elements count around"] < 4: - options["Elements count around"] = 4 + dependentChanges = False + + # Parameters specific to the core + if options["Core"] == True: + options, dependentChanges = checkCoreParameters(cls, options, dependentChanges) + + # When core option is false + else: + if options["Elements count around"] < 4: + options["Elements count around"] = 4 + annotationElementsCountsAround = options["Annotation elements counts around"] + if len(annotationElementsCountsAround) == 0: + options["Annotation elements count around"] = [0] + else: + for i in range(len(annotationElementsCountsAround)): + if annotationElementsCountsAround[i] <= 0: + annotationElementsCountsAround[i] = 0 + elif annotationElementsCountsAround[i] < 4: + annotationElementsCountsAround[i] = 4 + + # Common parameters if options["Elements count through wall"] < 1: options["Elements count through wall"] = 1 - annotationElementsCountsAround = options["Annotation elements counts around"] - if len(annotationElementsCountsAround) == 0: - options["Annotation elements count around"] = [0] - else: - for i in range(len(annotationElementsCountsAround)): - if annotationElementsCountsAround[i] <= 0: - annotationElementsCountsAround[i] = 0 - elif annotationElementsCountsAround[i] < 4: - annotationElementsCountsAround[i] = 4 + if options["Target element density along longest segment"] < 1.0: options["Target element density along longest segment"] = 1.0 - dependentChanges = False + return dependentChanges @classmethod @@ -112,6 +202,7 @@ def generateBaseMesh(cls, region, options): layoutRegion = region.createRegion() networkLayout = options["Network layout"] networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters + layoutAnnotationGroups = networkLayout.getAnnotationGroups() networkMesh = networkLayout.getConstructionObject() tubeNetworkMeshBuilder = TubeNetworkMeshBuilder( @@ -119,8 +210,13 @@ def generateBaseMesh(cls, region, options): targetElementDensityAlongLongestSegment=options["Target element density along longest segment"], defaultElementsCountAround=options["Elements count around"], elementsCountThroughWall=options["Elements count through wall"], - layoutAnnotationGroups=networkLayout.getAnnotationGroups(), - annotationElementsCountsAround=options["Annotation elements counts around"]) + layoutAnnotationGroups=layoutAnnotationGroups, + annotationElementsCountsAround=options["Annotation elements counts around"], + defaultElementsCountAcrossMajor=options['Number of elements across core major'], + elementsCountTransition=options['Number of elements across core transition'], + annotationElementsCountsAcrossMajor=options["Annotation elements counts across major"], + isCore=options["Core"]) + tubeNetworkMeshBuilder.build() generateData = TubeNetworkMeshGenerateData( region, 3, diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py new file mode 100644 index 00000000..29cc017c --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py @@ -0,0 +1,418 @@ +""" +Generates a 3D body coordinates using tube network mesh. +""" + +from cmlibs.zinc.element import Element +from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1 +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.scaffoldpackage import ScaffoldPackage +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, \ + getAnnotationGroupForTerm +from scaffoldmaker.annotation.body_terms import get_body_term +from scaffoldmaker.utils.tubenetworkmesh import TubeNetworkMeshBuilder, TubeNetworkMeshGenerateData +from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters +from cmlibs.zinc.node import Node + + +def getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName): + assert parameterSetName in cls.getParameterSetNames() + if parameterSetName in ("Default", "Human 1"): + return ScaffoldPackage(MeshType_1d_network_layout1, { + "scaffoldSettings": { + "Structure": "1-2-3, 3-4,4.2-5,4.3-6,4.1-7,7-8,8-9.1,9.2-10,9.3-11, \ + 5-12-13-14,6-15-16-17,10-18-19-20, 11-21-22-23", + "Define inner coordinates": True + }, + "meshEdits": exnode_string_from_nodeset_field_parameters( + ["coordinates", "inner coordinates"], + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [[ + (1, [[0.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, -0.000]]), + (2, [[0.500, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, -0.000]]), + (3, [[1.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, 0.100, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, 0.100]]), + (4, [[1.500, 0.000, 0.000], [[0.643, 0.000, 0.000], [0.028, -1.246, 0.000], [0.028, 1.246, 0.000]], [[0.000, 0.600, 0.000], [0.396, 0.009, 0.000], [-0.395, 0.009, 0.000]], [[0.000, 0.100, 0.000], [-0.374, 0.001, 0.000], [0.376, 0.005, 0.000]], [[0.000, 0.000, 0.600], [-0.000, 0.000, 0.400], [0.000, -0.000, 0.400]], [[0.000, 0.000, 0.100], [0.000, 0.000, -0.270], [0.000, 0.000, -0.268]]]), + (5, [[1.788, -1.014, 0.000], [0.546, -0.709, 0.000], [0.121, 0.093, 0.000], [-0.189, 0.048, 0.000], [-0.000, 0.000, 0.200], [0.000, 0.000, -0.130]]), + (6, [[1.788, 1.014, 0.000], [0.547, 0.709, 0.000], [-0.120, 0.093, 0.000], [0.186, 0.045, 0.000], [0.000, -0.000, 0.201], [0.000, 0.000, -0.130]]), + (7, [[2.400, -0.000, 0.000], [0.788, 0.000, 0.000], [0.000, 0.600, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.600], [0.000, 0.000, -0.000]]), + (8, [[3.100, -0.000, 0.000], [0.747, 0.000, 0.000], [0.000, 0.600, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.600], [0.000, 0.000, -0.000]]), + (9, [[3.900, -0.000, 0.000], [[0.853, 0.000, 0.000], [0.731, -0.720, 0.000], [0.731, 0.720, 0.000]], [[0.000, 0.600, 0.000], [0.496, 0.504, 0.000], [-0.470, 0.478, 0.000]], [[-0.000, 0.000, -0.000], [-0.839, -0.173, -0.006], [0.775, -0.153, 0.006]], [[0.000, 0.000, 0.600], [-0.000, 0.000, 0.350], [0.000, -0.000, 0.350]], [[-0.000, -0.000, 0.000], [0.006, 0.004, -0.022], [0.006, -0.004, -0.022]]]), + (10, [[5.250, -0.600, 0.000], [1.915, -0.427, 0.000], [0.056, 0.249, 0.000], [-0.211, -0.182, 0.006], [-0.000, 0.000, 0.300], [-0.003, -0.005, -0.078]]), + (11, [[5.250, 0.600, 0.000], [1.915, 0.427, 0.000], [-0.056, 0.249, 0.000], [0.206, -0.165, -0.006], [0.000, -0.000, 0.300], [-0.003, 0.005, -0.078]]), + (12, [[2.439, -1.399, 0.000], [0.920, -0.382, 0.000], [0.057, 0.137, 0.000], [-0.009, -0.011, 0.001], [-0.000, 0.000, 0.140], [0.006, -0.002, -0.030]]), + (13, [[3.807, -1.755, 0.000], [1.181, -0.272, 0.078], [0.028, 0.121, -0.000], [-0.018, -0.063, 0.000], [-0.009, 0.002, 0.140], [-0.023, 0.005, -0.056]]), + (14, [[4.841, -1.961, 0.137], [0.882, -0.141, 0.194], [0.002, 0.015, -0.000], [-0.023, -0.149, 0.002], [-0.006, 0.001, 0.029], [0.054, -0.011, -0.158]]), + (15, [[2.439, 1.399, 0.000], [0.946, 0.385, 0.000], [-0.056, 0.138, 0.000], [0.012, -0.011, -0.001], [0.000, -0.000, 0.140], [0.006, 0.002, -0.031]]), + (16, [[3.911, 1.761, 0.000], [1.211, 0.255, 0.091], [-0.025, 0.121, 0.000], [0.018, -0.056, -0.000], [-0.010, -0.002, 0.140], [-0.028, -0.005, -0.042]]), + (17, [[4.932, 1.939, 0.153], [0.825, 0.101, 0.214], [-0.004, 0.029, 0.000], [0.014, -0.129, -0.002], [-0.015, -0.002, 0.058], [0.044, 0.008, -0.114]]), + (18, [[8.100, -0.665, 0.000], [0.775, -0.027, 0.275], [0.007, 0.250, 0.003], [0.032, -0.000, 0.006], [-0.069, -0.000, 0.194], [-0.224, -0.006, -0.056]]), + (19, [[8.459, -0.683, 0.298], [0.258, -0.010, 0.370], [0.004, 0.250, 0.004], [-0.008, 0.000, -0.013], [-0.181, 0.001, 0.126], [0.051, -0.005, -0.001]]), + (20, [[8.601, -0.685, 0.694], [0.024, 0.006, 0.399], [-0.001, 0.250, -0.004], [-0.002, -0.000, 0.005], [-0.141, -0.000, 0.009], [0.109, 0.005, -0.120]]), + (21, [[8.100, 0.665, 0.000], [0.775, 0.027, 0.275], [-0.007, 0.250, -0.003], [-0.032, -0.000, -0.009], [-0.069, 0.000, 0.194], [-0.224, 0.008, -0.055]]), + (22, [[8.459, 0.683, 0.298], [0.258, 0.010, 0.370], [0.000, 0.250, -0.007], [0.008, 0.000, 0.013], [-0.182, 0.004, 0.127], [0.049, 0.005, 0.002]]), + (23, [[8.601, 0.685, 0.694], [0.024, -0.006, 0.399], [0.001, 0.250, 0.004], [-0.007, 0.000, -0.005], [-0.148, 0.000, 0.009], [0.101, -0.015, -0.124]])], [ + (1, [[0.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, 0.000, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, 0.000]]), + (2, [[0.500, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, -0.000]]), + (3, [[1.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, 0.070, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, 0.070]]), + (4, [[1.500, 0.000, 0.000], [[0.643, 0.000, 0.000], [0.028, -1.246, 0.000], [0.028, 1.246, 0.000]], [[0.000, 0.420, 0.000], [0.277, 0.006, 0.000], [-0.277, 0.006, 0.000]], [[0.000, 0.070, 0.000], [-0.262, 0.001, 0.000], [0.263, 0.003, 0.000]], [[0.000, 0.000, 0.420], [-0.000, 0.000, 0.280], [0.000, -0.000, 0.280]], [[0.000, 0.000, 0.070], [0.000, 0.000, -0.189], [0.000, 0.000, -0.188]]]), + (5, [[1.788, -1.014, 0.000], [0.546, -0.709, 0.000], [0.085, 0.065, 0.000], [-0.133, 0.033, 0.000], [-0.000, 0.000, 0.140], [0.000, 0.000, -0.091]]), + (6, [[1.788, 1.014, 0.000], [0.547, 0.709, 0.000], [-0.084, 0.065, 0.000], [0.130, 0.032, 0.000], [0.000, -0.000, 0.141], [0.000, 0.000, -0.091]]), + (7, [[2.400, -0.000, 0.000], [0.788, 0.000, 0.000], [0.000, 0.420, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.420], [0.000, 0.000, 0.000]]), + (8, [[3.100, -0.000, 0.000], [0.747, 0.000, 0.000], [0.000, 0.420, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.420], [0.000, 0.000, -0.000]]), + (9, [[3.900, -0.000, 0.000], [[0.853, 0.000, 0.000], [0.731, -0.720, 0.000], [0.731, 0.720, 0.000]], [[0.000, 0.420, 0.000], [0.347, 0.353, 0.000], [-0.329, 0.334, 0.000]], [[-0.000, 0.000, -0.000], [-0.587, -0.121, -0.004], [0.543, -0.107, 0.004]], [[0.000, 0.000, 0.420], [-0.000, 0.000, 0.245], [0.000, -0.000, 0.245]], [[-0.000, -0.000, 0.000], [0.004, 0.003, -0.015], [0.004, -0.003, -0.015]]]), + (10, [[5.250, -0.600, 0.000], [1.915, -0.427, 0.000], [0.039, 0.174, 0.000], [-0.148, -0.127, 0.004], [-0.000, 0.000, 0.210], [-0.002, -0.004, -0.055]]), + (11, [[5.250, 0.600, 0.000], [1.915, 0.427, 0.000], [-0.039, 0.174, 0.000], [0.144, -0.115, -0.004], [0.000, -0.000, 0.210], [-0.002, 0.004, -0.055]]), + (12, [[2.439, -1.399, 0.000], [0.920, -0.382, 0.000], [0.040, 0.096, 0.000], [-0.006, -0.008, 0.001], [-0.000, 0.000, 0.098], [0.004, -0.002, -0.021]]), + (13, [[3.807, -1.755, 0.000], [1.181, -0.272, 0.078], [0.019, 0.084, -0.000], [-0.013, -0.044, 0.000], [-0.006, 0.001, 0.098], [-0.016, 0.003, -0.039]]), + (14, [[4.841, -1.961, 0.137], [0.882, -0.141, 0.194], [0.002, 0.011, -0.000], [-0.016, -0.104, 0.001], [-0.004, 0.001, 0.021], [0.038, -0.008, -0.110]]), + (15, [[2.439, 1.399, 0.000], [0.946, 0.385, 0.000], [-0.039, 0.096, 0.000], [0.008, -0.007, -0.001], [0.000, -0.000, 0.098], [0.004, 0.002, -0.022]]), + (16, [[3.911, 1.761, 0.000], [1.211, 0.255, 0.091], [-0.018, 0.085, 0.000], [0.013, -0.039, -0.000], [-0.007, -0.002, 0.098], [-0.019, -0.004, -0.029]]), + (17, [[4.932, 1.939, 0.153], [0.825, 0.101, 0.214], [-0.003, 0.020, 0.000], [0.010, -0.090, -0.002], [-0.010, -0.002, 0.041], [0.031, 0.006, -0.080]]), + (18, [[8.100, -0.665, 0.000], [0.775, -0.027, 0.275], [0.005, 0.175, 0.002], [0.022, -0.000, 0.004], [-0.048, -0.000, 0.136], [-0.157, -0.004, -0.039]]), + (19, [[8.459, -0.683, 0.298], [0.258, -0.010, 0.370], [0.003, 0.175, 0.003], [-0.005, 0.000, -0.009], [-0.127, 0.001, 0.088], [0.036, -0.004, -0.000]]), + (20, [[8.601, -0.685, 0.694], [0.024, 0.006, 0.399], [-0.001, 0.175, -0.002], [-0.001, -0.000, 0.003], [-0.099, -0.000, 0.006], [0.077, 0.004, -0.084]]), + (21, [[8.100, 0.665, 0.000], [0.775, 0.027, 0.275], [-0.005, 0.175, -0.002], [-0.022, -0.000, -0.006], [-0.048, 0.000, 0.136], [-0.157, 0.006, -0.039]]), + (22, [[8.459, 0.683, 0.298], [0.258, 0.010, 0.370], [0.000, 0.175, -0.005], [0.005, 0.000, 0.009], [-0.127, 0.002, 0.089], [0.034, 0.004, 0.001]]), + (23, [[8.601, 0.685, 0.694], [0.024, -0.006, 0.399], [0.001, 0.175, 0.002], [-0.005, 0.000, -0.003], [-0.104, 0.000, 0.006], [0.071, -0.010, -0.087]]) + ]]), + "userAnnotationGroups": [ + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '1-3', + 'name': get_body_term('head')[0], + 'ontId': get_body_term('head')[1] + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '6-8', + 'name': 'torso', + 'ontId': 'None' + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '4,11-13, 5,14-16', + 'name': get_body_term('arm')[0], + 'ontId': get_body_term('arm')[1] + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '9,17-18-19, 10,20-21-22', + 'name': get_body_term('leg')[0], + 'ontId': get_body_term('leg')[1] + } + ] + }) + elif parameterSetName in ("Human 2 Coarse", "Human 2 Medium", "Human 2 Fine"): + return ScaffoldPackage(MeshType_1d_network_layout1, { + "scaffoldSettings": { + "Structure": "1-2-3, 3-4,4.2-5,4.3-6,4.1-7,7-8,8-9.1,9.2-10,9.3-11, \ + 5-12-13-14,6-15-16-17,10-18-19-20, 11-21-22-23", + "Define inner coordinates": True + }, + "meshEdits": exnode_string_from_nodeset_field_parameters( + ["coordinates", "inner coordinates"], + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [[ + (1, [[0.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, -0.000]]), + (2, [[0.500, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, -0.000]]), + (3, [[1.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.400, 0.000], [0.000, 0.100, 0.000], [0.000, 0.000, 0.400], [0.000, 0.000, 0.100]]), + (4, [[1.500, 0.000, 0.000], [[0.643, 0.000, 0.000], [0.012, -1.338, 0.000], [0.054, 1.358, 0.000]], [[0.000, 0.600, 0.000], [0.277, 0.003, 0.000], [-0.277, 0.011, 0.000]], [[0.000, 0.100, 0.000], [-0.083, -0.206, 0.000], [0.082, -0.198, 0.000]], [[0.000, 0.000, 0.600], [-0.000, 0.000, 0.280], [0.000, -0.000, 0.280]], [[0.000, 0.000, 0.100], [0.000, 0.000, -0.195], [0.000, 0.000, -0.195]]]), + (5, [[1.806, -1.023, 0.000], [0.599, -0.599, 0.000], [0.158, 0.158, 0.000], [-0.253, 0.182, 0.000], [-0.000, 0.000, 0.150], [0.000, 0.000, -0.065]]), + (6, [[1.825, 1.023, 0.000], [0.592, 0.587, 0.000], [-0.156, 0.158, 0.000], [0.242, 0.173, 0.000], [0.000, -0.000, 0.150], [0.000, 0.000, -0.065]]), + (7, [[2.400, -0.000, 0.000], [0.788, 0.000, 0.000], [0.000, 0.600, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.600], [0.000, 0.000, -0.000]]), + (8, [[3.100, -0.000, 0.000], [0.747, 0.000, 0.000], [0.000, 0.600, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.600], [0.000, 0.000, -0.000]]), + (9, [[3.900, -0.000, 0.000], [[0.853, 0.000, 0.000], [0.731, -0.720, 0.000], [0.731, 0.720, 0.000]], [[0.000, 0.600, 0.000], [0.496, 0.504, 0.000], [-0.470, 0.478, 0.000]], [[-0.000, 0.000, -0.000], [-0.839, -0.173, -0.006], [0.776, -0.153, 0.006]], [[0.000, 0.000, 0.600], [-0.000, 0.000, 0.350], [0.000, -0.000, 0.350]], [[-0.000, -0.000, 0.000], [0.006, 0.004, -0.022], [0.006, -0.004, -0.022]]]), + (10, [[5.250, -0.600, 0.000], [1.915, -0.427, 0.000], [0.056, 0.249, 0.000], [-0.211, -0.182, 0.006], [-0.000, 0.000, 0.300], [-0.003, -0.005, -0.078]]), + (11, [[5.250, 0.600, 0.000], [1.915, 0.427, 0.000], [-0.056, 0.249, 0.000], [0.206, -0.165, -0.006], [0.000, -0.000, 0.300], [-0.003, 0.005, -0.078]]), + (12, [[2.450, -1.213, 0.000], [1.000, -0.145, 0.000], [0.022, 0.151, 0.000], [0.031, -0.007, 0.000], [-0.000, 0.000, 0.150], [0.000, 0.000, -0.015]]), + (13, [[4.327, -1.224, 0.000], [0.756, 0.003, 0.000], [-0.001, 0.200, 0.000], [-0.001, 0.075, -0.001], [0.000, -0.000, 0.120], [-0.000, 0.000, -0.025]]), + (14, [[4.800, -1.218, 0.000], [0.190, 0.009, 0.000], [-0.015, 0.300, -0.001], [-0.045, 0.124, -0.001], [-0.000, 0.000, 0.100], [-0.000, 0.000, -0.015]]), + (15, [[2.450, 1.213, 0.000], [0.979, 0.146, 0.000], [-0.023, 0.151, 0.000], [-0.029, -0.008, 0.000], [0.000, -0.000, 0.150], [0.000, 0.000, -0.015]]), + (16, [[4.327, 1.224, 0.000], [0.756, -0.003, 0.000], [0.001, 0.200, 0.000], [0.001, 0.075, 0.001], [-0.000, 0.000, 0.120], [-0.000, -0.000, -0.025]]), + (17, [[4.800, 1.218, 0.000], [0.190, -0.009, 0.000], [0.015, 0.300, 0.001], [0.045, 0.124, 0.001], [-0.000, -0.000, 0.100], [-0.000, -0.000, -0.015]]), + (18, [[8.100, -0.665, 0.000], [0.775, -0.027, 0.275], [0.007, 0.250, 0.003], [0.032, -0.000, 0.006], [-0.069, -0.000, 0.194], [-0.224, -0.006, -0.056]]), + (19, [[8.459, -0.683, 0.298], [0.258, -0.010, 0.370], [0.004, 0.250, 0.004], [-0.008, 0.000, -0.013], [-0.181, 0.001, 0.126], [0.051, -0.006, -0.001]]), + (20, [[8.601, -0.685, 0.694], [0.024, 0.006, 0.399], [-0.001, 0.250, -0.004], [-0.003, -0.000, 0.005], [-0.141, -0.000, 0.009], [0.111, 0.004, -0.120]]), + (21, [[8.100, 0.665, 0.000], [0.775, 0.027, 0.275], [-0.007, 0.250, -0.003], [-0.032, -0.000, -0.008], [-0.069, 0.000, 0.194], [-0.224, 0.008, -0.055]]), + (22, [[8.459, 0.683, 0.298], [0.258, 0.010, 0.370], [-0.000, 0.250, -0.007], [0.008, 0.000, 0.013], [-0.182, 0.003, 0.127], [0.049, 0.005, 0.002]]), + (23, [[8.601, 0.685, 0.694], [0.024, -0.006, 0.399], [0.001, 0.250, 0.004], [-0.006, 0.000, -0.005], [-0.148, 0.000, 0.009], [0.101, -0.013, -0.124]])], [ + (1, [[0.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, 0.000, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, 0.000]]), + (2, [[0.500, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, -0.000]]), + (3, [[1.000, 0.000, 0.000], [0.500, 0.000, 0.000], [0.000, 0.280, 0.000], [0.000, 0.070, 0.000], [0.000, 0.000, 0.280], [0.000, 0.000, 0.070]]), + (4, [[1.500, 0.000, 0.000], [[0.643, 0.000, 0.000], [0.012, -1.338, 0.000], [0.054, 1.358, 0.000]], [[0.000, 0.420, 0.000], [0.194, 0.002, 0.000], [-0.194, 0.008, 0.000]], [[0.000, 0.070, 0.000], [-0.058, -0.144, 0.000], [0.057, -0.139, 0.000]], [[0.000, 0.000, 0.420], [-0.000, 0.000, 0.196], [0.000, -0.000, 0.196]], [[0.000, 0.000, 0.070], [0.000, 0.000, -0.137], [0.000, 0.000, -0.137]]]), + (5, [[1.806, -1.023, 0.000], [0.599, -0.599, 0.000], [0.111, 0.111, 0.000], [-0.177, 0.128, 0.000], [-0.000, 0.000, 0.105], [0.000, 0.000, -0.045]]), + (6, [[1.825, 1.023, 0.000], [0.592, 0.587, 0.000], [-0.109, 0.110, 0.000], [0.169, 0.121, 0.000], [0.000, -0.000, 0.105], [0.000, 0.000, -0.045]]), + (7, [[2.400, -0.000, 0.000], [0.788, 0.000, 0.000], [0.000, 0.420, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.420], [0.000, 0.000, 0.000]]), + (8, [[3.100, -0.000, 0.000], [0.747, 0.000, 0.000], [0.000, 0.420, 0.000], [0.000, -0.000, 0.000], [0.000, 0.000, 0.420], [0.000, 0.000, -0.000]]), + (9, [[3.900, -0.000, 0.000], [[0.853, 0.000, 0.000], [0.731, -0.720, 0.000], [0.731, 0.720, 0.000]], [[0.000, 0.420, 0.000], [0.347, 0.353, 0.000], [-0.329, 0.335, 0.000]], [[-0.000, 0.000, -0.000], [-0.587, -0.121, -0.004], [0.543, -0.107, 0.004]], [[0.000, 0.000, 0.420], [-0.000, 0.000, 0.245], [0.000, -0.000, 0.245]], [[-0.000, -0.000, 0.000], [0.004, 0.003, -0.015], [0.004, -0.003, -0.015]]]), + (10, [[5.250, -0.600, 0.000], [1.915, -0.427, 0.000], [0.039, 0.174, 0.000], [-0.147, -0.127, 0.004], [-0.000, 0.000, 0.210], [-0.002, -0.004, -0.055]]), + (11, [[5.250, 0.600, 0.000], [1.915, 0.427, 0.000], [-0.039, 0.174, 0.000], [0.144, -0.115, -0.004], [0.000, -0.000, 0.210], [-0.002, 0.004, -0.055]]), + (12, [[2.450, -1.213, 0.000], [1.000, -0.145, 0.000], [0.015, 0.106, 0.000], [0.022, -0.005, 0.000], [-0.000, 0.000, 0.105], [0.000, 0.000, -0.011]]), + (13, [[4.327, -1.224, 0.000], [0.756, 0.003, 0.000], [-0.000, 0.140, 0.000], [-0.000, 0.052, -0.000], [0.000, -0.000, 0.084], [-0.000, 0.000, -0.018]]), + (14, [[4.800, -1.218, 0.000], [0.190, 0.009, 0.000], [-0.010, 0.210, -0.001], [-0.032, 0.087, -0.001], [-0.000, 0.000, 0.070], [-0.000, 0.000, -0.011]]), + (15, [[2.450, 1.213, 0.000], [0.979, 0.146, 0.000], [-0.016, 0.106, 0.000], [-0.020, -0.006, 0.000], [0.000, -0.000, 0.105], [0.000, 0.000, -0.011]]), + (16, [[4.327, 1.224, 0.000], [0.756, -0.003, 0.000], [0.000, 0.140, 0.000], [0.001, 0.052, 0.000], [-0.000, 0.000, 0.084], [-0.000, -0.000, -0.018]]), + (17, [[4.800, 1.218, 0.000], [0.190, -0.009, 0.000], [0.010, 0.210, 0.001], [0.031, 0.087, 0.001], [-0.000, -0.000, 0.070], [-0.000, -0.000, -0.011]]), + (18, [[8.100, -0.665, 0.000], [0.775, -0.027, 0.275], [0.005, 0.175, 0.002], [0.022, -0.000, 0.004], [-0.048, -0.000, 0.136], [-0.157, -0.004, -0.039]]), + (19, [[8.459, -0.683, 0.298], [0.258, -0.010, 0.370], [0.003, 0.175, 0.003], [-0.005, 0.000, -0.009], [-0.127, 0.001, 0.089], [0.036, -0.004, -0.000]]), + (20, [[8.601, -0.685, 0.694], [0.024, 0.006, 0.399], [-0.001, 0.175, -0.002], [-0.002, -0.000, 0.003], [-0.099, -0.000, 0.006], [0.077, 0.003, -0.084]]), + (21, [[8.100, 0.665, 0.000], [0.775, 0.027, 0.275], [-0.005, 0.175, -0.002], [-0.022, -0.000, -0.006], [-0.048, 0.000, 0.136], [-0.157, 0.006, -0.039]]), + (22, [[8.459, 0.683, 0.298], [0.258, 0.010, 0.370], [-0.000, 0.175, -0.005], [0.005, 0.000, 0.009], [-0.127, 0.002, 0.089], [0.034, 0.004, 0.001]]), + (23, [[8.601, 0.685, 0.694], [0.024, -0.006, 0.399], [0.001, 0.175, 0.002], [-0.004, 0.000, -0.003], [-0.104, 0.000, 0.006], [0.071, -0.009, -0.087]]) + ]]), + "userAnnotationGroups": [ + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '1-3', + 'name': get_body_term('head')[0], + 'ontId': get_body_term('head')[1] + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '6-8', + 'name': 'torso', + 'ontId': 'None' + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '4,11-13, 5,14-16', + 'name': get_body_term('arm')[0], + 'ontId': get_body_term('arm')[1] + }, + { + '_AnnotationGroup': True, + 'dimension': 1, + 'identifierRanges': '9,17-18-19, 10,20-21-22', + 'name': get_body_term('leg')[0], + 'ontId': get_body_term('leg')[1] + } + ] + }) + + +class MeshType_3d_wholebody2(Scaffold_base): + """ + Generates a 3-D hermite bifurcating tube network, with linear basis through wall. + """ + + @staticmethod + def getName(): + return "3D Whole Body 2" + + @classmethod + def getParameterSetNames(cls): + return [ + "Default", + "Human 1", + "Human 2 Coarse", + "Human 2 Medium", + "Human 2 Fine" + ] + + @classmethod + def getDefaultOptions(cls, parameterSetName="Default"): + + options = { + 'Base parameter set': parameterSetName, + "Network layout": getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName), + "Number of elements around head": 12, + "Number of elements around torso": 12, + "Number of elements around arm": 8, + "Number of elements around leg": 8, + "Number of elements through wall": 1, + "Target element density along longest segment": 5.0, + "Show trim surfaces": False, + "Use Core": True, + "Number of elements across head core major": 6, + "Number of elements across torso core major": 6, + "Number of elements across arm core major": 4, + "Number of elements across leg core major": 4, + "Number of elements across core transition": 1 + } + + if "Human 2 Medium" in parameterSetName: + options["Number of elements around head"] = 16 + options["Number of elements around torso"] = 16 + options["Number of elements around arm"] = 12 + options["Number of elements around leg"] = 12 + + options["Number of elements across head core major"] = 8 + options["Number of elements across torso core major"] = 8 + options["Number of elements across arm core major"] = 6 + options["Number of elements across leg core major"] = 6 + + elif "Human 2 Fine" in parameterSetName: + options["Number of elements around head"] = 24 + options["Number of elements around torso"] = 24 + options["Number of elements around arm"] = 20 + options["Number of elements around leg"] = 20 + + options["Number of elements across head core major"] = 10 + options["Number of elements across torso core major"] = 10 + options["Number of elements across arm core major"] = 8 + options["Number of elements across leg core major"] = 8 + + return options + + @classmethod + def getOrderedOptionNames(cls): + optionNames = [ + "Network layout", + "Number of elements around head", + "Number of elements around torso", + "Number of elements around arm", + "Number of elements around leg", + "Number of elements through wall", + "Target element density along longest segment", + "Show trim surfaces", + "Use Core", + "Number of elements across head core major", + "Number of elements across torso core major", + "Number of elements across arm core major", + "Number of elements across leg core major"] + return optionNames + + @classmethod + def getOptionValidScaffoldTypes(cls, optionName): + if optionName == "Network layout": + return [MeshType_1d_network_layout1] + return [] + + @classmethod + def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType): + if optionName == "Network layout": + return cls.getParameterSetNames() + assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), \ + cls.__name__ + ".getOptionScaffoldTypeParameterSetNames. " + \ + "Invalid option \"" + optionName + "\" scaffold type " + scaffoldType.getName() + return scaffoldType.getParameterSetNames() + + @classmethod + def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): + """ + :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. + :return: ScaffoldPackage. + """ + if parameterSetName: + assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ + "Invalid parameter set " + str(parameterSetName) + " for scaffold " + str(scaffoldType.getName()) + \ + " in option " + str(optionName) + " of scaffold " + cls.getName() + if optionName == "Network layout": + if not parameterSetName: + parameterSetName = "Default" + return getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName) + assert False, cls.__name__ + ".getOptionScaffoldPackage: Option " + optionName + " is not a scaffold" + + @classmethod + def checkOptions(cls, options): + if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"): + options["Network layout"] = cls.getOptionScaffoldPackage('Network layout', MeshType_1d_network_layout1) + for key in [ + "Number of elements around head", + "Number of elements around torso", + "Number of elements around arm", + "Number of elements around leg" + ]: + if options[key] < 8: + options[key] = 8 + + for key in [ + "Number of elements across head core major", + "Number of elements across torso core major", + "Number of elements across arm core major", + "Number of elements across leg core major" + ]: + if options[key] < 4: + options[key] = 4 + + if options["Number of elements through wall"] < 0: + options["Number of elements through wall"] = 1 + + if options["Target element density along longest segment"] < 1.0: + options["Target element density along longest segment"] = 1.0 + + dependentChanges = False + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options): + """ + Generate the base hermite-bilinear mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup, None + """ + parameterSetName = options['Base parameter set'] + isHuman = parameterSetName in ["Default", "Human 1", "Human 2 Coarse", "Human 2 Medium", "Human 2 Fine"] + + layoutRegion = region.createRegion() + networkLayout = options["Network layout"] + networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters + layoutAnnotationGroups = networkLayout.getAnnotationGroups() + networkMesh = networkLayout.getConstructionObject() + + annotationElementsCountsAround = [] + annotationElementsCountsAcross = [] + for layoutAnnotationGroup in layoutAnnotationGroups: + elementsCountAround = 0 + elementsCountAcrossMajor = 0 + name = layoutAnnotationGroup.getName() + if name in ["head", "torso", "arm", "leg"]: + elementsCountAround = options["Number of elements around " + name] + elementsCountAcrossMajor = options["Number of elements across " + name + " core major"] + annotationElementsCountsAround.append(elementsCountAround) + annotationElementsCountsAcross.append(elementsCountAcrossMajor) + + isCore = options["Use Core"] + + tubeNetworkMeshBuilder = TubeNetworkMeshBuilder( + networkMesh, + targetElementDensityAlongLongestSegment=options["Target element density along longest segment"], + defaultElementsCountAround=options["Number of elements around head"], + elementsCountThroughWall=options["Number of elements through wall"], + layoutAnnotationGroups=layoutAnnotationGroups, + annotationElementsCountsAround=annotationElementsCountsAround, + defaultElementsCountAcrossMajor=options['Number of elements across head core major'], + elementsCountTransition=options['Number of elements across core transition'], + annotationElementsCountsAcrossMajor=annotationElementsCountsAcross, + isCore=isCore) + + tubeNetworkMeshBuilder.build() + generateData = TubeNetworkMeshGenerateData( + region, 3, + isLinearThroughWall=False, + isShowTrimSurfaces=options["Show trim surfaces"]) + tubeNetworkMeshBuilder.generateMesh(generateData) + annotationGroups = generateData.getAnnotationGroups() + + return annotationGroups, None + + @classmethod + def defineFaceAnnotations(cls, region, options, annotationGroups): + """ + Add face annotation groups from the highest dimension mesh. + Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. + :param options: Dict containing options. See getDefaultOptions(). + :param annotationGroups: List of annotation groups for top-level elements. + New face annotation groups are appended to this list. + """ + + # create 2d surface mesh groups + fm = region.getFieldmodule() + mesh2d = fm.findMeshByDimension(2) + + skinGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_body_term("skin epidermis")) + + is_exterior = fm.createFieldIsExterior() + is_on_face_xi3_1 = fm.createFieldIsOnFace(Element.FACE_TYPE_XI3_1) + is_skin = fm.createFieldAnd(is_exterior, is_on_face_xi3_1) + + skinMeshGroup = skinGroup.getMeshGroup(mesh2d) + skinMeshGroup.addElementsConditional(is_skin) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index 519ec66e..68dfd0de 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -56,6 +56,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_uterus1 import MeshType_3d_uterus1 from scaffoldmaker.meshtypes.meshtype_3d_uterus2 import MeshType_3d_uterus2 from scaffoldmaker.meshtypes.meshtype_3d_wholebody1 import MeshType_3d_wholebody1 +from scaffoldmaker.meshtypes.meshtype_3d_wholebody2 import MeshType_3d_wholebody2 from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -114,7 +115,8 @@ def __init__(self): MeshType_3d_tubeseptum1, MeshType_3d_uterus1, MeshType_3d_uterus2, - MeshType_3d_wholebody1 + MeshType_3d_wholebody1, + MeshType_3d_wholebody2 ] def findScaffoldTypeByName(self, name): diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index c3ec1f19..713913ad 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -629,12 +629,35 @@ def __init__(self): self._nodeLayout6Way12_d3Defined = HermiteNodeLayout( [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]) + self._nodeLayout6WayBifurcation = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -1.0, 0.0], [0.0, -1.0, 0.0], + [0.0, 0.0, -1.0], [0.0, 0.0, 1.0], [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0]]) + self._nodeLayout6WayBifurcationTransition = HermiteNodeLayout( + [[0.0, 0.0, -1.0], [0.0, 0.0, 1.0], [1.0, 1.0, 0.0], [-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0]]) + self._nodeLayout6WayTriplePointTop1 = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, -1.0, 0.0], + [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]]) + self._nodeLayout6WayTriplePointTop2 = HermiteNodeLayout( + [[-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]]) + self._nodeLayout6WayTriplePointBottom1 = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, -1.0, 0.0], + [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 1.0, -1.0]]) + self._nodeLayout6WayTriplePointBottom2 = HermiteNodeLayout( + [[-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 1.0, -1.0]]) self._nodeLayout8Way12 = HermiteNodeLayout( [[1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [-1.0, 1.0], [-1.0, 0.0], [-1.0, -1.0], [0.0, -1.0], [1.0, -1.0]]) self._nodeLayout8Way12_d3Defined = HermiteNodeLayout( [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -1.0, 0.0], [0.0, -1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]) + self._nodeLayoutTriplePointTopLeft = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0]]) + self._nodeLayoutTriplePointTopRight = HermiteNodeLayout( + [[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 0.0, 1.0]]) + self._nodeLayoutTriplePointBottomLeft = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [-1.0, 0.0, -1.0]]) + self._nodeLayoutTriplePointBottomRight = HermiteNodeLayout( + [[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, -1.0]]) def getNodeLayoutRegularPermuted(self, d3Defined, limitDirections=None): """ @@ -652,13 +675,20 @@ def getNodeLayoutRegularPermuted(self, d3Defined, limitDirections=None): nodeLayout = HermiteNodeLayout(None, nodeLayout, limitDirections) return nodeLayout - def getNodeLayout6Way12(self, d3Defined): + def getNodeLayout6Way12(self, d3Defined, limitDirections=None): """ Get node layout for 6-way junction in 1-2 plane, including d1 + d2, -d1 - d2. :param d3Defined: Set to True to use tricubic variant with d3 defined, otherwise bicubic is used. + :param limitDirections: Optional list over element directions of lists of allowable weights for that + direction, or None to not filter. Default None for whole list does not filter any directions. + For example, with d3, [None, [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], [[0.0, 0.0, 1.0]]] places no + limits on the first derivative, but derivative 2 must be [0, +/-1, 0] and d3 must be [0, 0, 1]. :return: HermiteNodeLayout. """ - return self._nodeLayout6Way12_d3Defined if d3Defined else self._nodeLayout6Way12 + nodeLayout = self._nodeLayout6Way12_d3Defined if d3Defined else self._nodeLayout6Way12 + if limitDirections: + nodeLayout = HermiteNodeLayout(None, nodeLayout, limitDirections) + return nodeLayout def getNodeLayout8Way12(self, d3Defined): """ @@ -668,6 +698,42 @@ def getNodeLayout8Way12(self, d3Defined): """ return self._nodeLayout8Way12_d3Defined if d3Defined else self._nodeLayout8Way12 + def getNodeLayoutTriplePoint(self): + """ + Get node layout for triple-point corners of core box elements. There are four corners (Top Left, Top Right, + Bottom Left, and Bottom Right) each with its specific node layout. + :return: HermiteNodeLayout. + """ + nodeLayouts = [self._nodeLayoutTriplePointTopLeft, self._nodeLayoutTriplePointTopRight, + self._nodeLayoutTriplePointBottomLeft, self._nodeLayoutTriplePointBottomRight] + return nodeLayouts + + def getNodeLayout6WayBifurcation(self): + """ + Get node layout for a special case of 6-way bifurcation, used in conjunction with a node layout for + 6-way bifurcation transition and a node layout for 6-way triple point. + :return: HermiteNodeLayout. + """ + return self._nodeLayout6WayBifurcation + + def getNodeLayout6WayBifurcationTransition(self): + """ + Get node layout for a special case of 6-way bifurcation transition, used in conjunction with a node layout for + 6-way bifurcation and a node layout for 6-way triple point. + :return: HermiteNodeLayout. + """ + return self._nodeLayout6WayBifurcationTransition + + def getNodeLayout6WayTriplePoint(self): + """ + Get node layout for a special case where a node is located at both 6-way junction and one of the triple-point + corners of core box elements. Used in conjunction with 6-way bifurcation and 6-way bifurcation transition node + layouts. There are two corners (Top, and Bottom) each with its specific pair of node layouts. + :return: HermiteNodeLayout. + """ + nodeLayouts = [self._nodeLayout6WayTriplePointTop1, self._nodeLayout6WayTriplePointTop2, + self._nodeLayout6WayTriplePointBottom1, self._nodeLayout6WayTriplePointBottom2] + return nodeLayouts def determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts): """ diff --git a/src/scaffoldmaker/utils/tubenetworkmesh.py b/src/scaffoldmaker/utils/tubenetworkmesh.py index 8f81388a..24ba4813 100644 --- a/src/scaffoldmaker/utils/tubenetworkmesh.py +++ b/src/scaffoldmaker/utils/tubenetworkmesh.py @@ -1,15 +1,16 @@ """ Specialisation of Network Mesh for building 2-D and 3-D tube mesh networks. """ -from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub, rejection +from cmlibs.maths.vectorops import add, cross, dot, magnitude, mult, normalize, set_magnitude, sub, rejection from cmlibs.zinc.element import Element, Elementbasis from cmlibs.zinc.node import Node from scaffoldmaker.utils.eft_utils import determineCubicHermiteSerendipityEft, HermiteNodeLayoutManager from scaffoldmaker.utils.interpolation import ( computeCubicHermiteDerivativeScaling, DerivativeScalingMode, evaluateCoordinatesOnCurve, getCubicHermiteTrimmedCurvesLengths, interpolateCubicHermite, interpolateCubicHermiteDerivative, - interpolateLagrangeHermiteDerivative, interpolateSampleCubicHermite, sampleCubicHermiteCurvesSmooth, - smoothCubicHermiteDerivativesLine, smoothCubicHermiteDerivativesLoop, smoothCurveSideCrossDerivatives) + interpolateLagrangeHermiteDerivative, interpolateSampleCubicHermite, sampleCubicHermiteCurves, + sampleCubicHermiteCurvesSmooth, smoothCubicHermiteDerivativesLine, smoothCubicHermiteDerivativesLoop, + smoothCurveSideCrossDerivatives) from scaffoldmaker.utils.networkmesh import NetworkMesh, NetworkMeshBuilder, NetworkMeshGenerateData, \ NetworkMeshJunction, NetworkMeshSegment, pathValueLabels from scaffoldmaker.utils.tracksurface import TrackSurface @@ -60,6 +61,14 @@ def __init__(self, region, meshDimension, isLinearThroughWall, isShowTrimSurface self._nodeLayoutFlipD2 = self._nodeLayoutManager.getNodeLayoutRegularPermuted( d3Defined, limitDirections=[None, [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], [[0.0, 0.0, 1.0]]] if d3Defined else [None, [[0.0, 1.0], [0.0, -1.0]]]) + self._nodeLayout6WayTriplePoint = self._nodeLayoutManager.getNodeLayout6WayTriplePoint() + self._nodeLayoutBifrucation = self._nodeLayoutManager.getNodeLayout6WayBifurcation() + self._nodeLayoutTrifurcation = None + self._nodeLayoutTransition = self._nodeLayoutManager.getNodeLayoutRegularPermuted( + d3Defined, limitDirections=[None, [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], None]) + self._nodeLayoutTransitionTriplePoint = None + self._nodeLayoutBifurcationTransition = self._nodeLayoutManager.getNodeLayout6WayBifurcationTransition() + def getStandardEft(self): return self._standardEft @@ -76,6 +85,86 @@ def getNodeLayout8Way(self): def getNodeLayoutFlipD2(self): return self._nodeLayoutFlipD2 + def getNodeLayout6WayTriplePoint(self, location): + """ + Special node layout for generating core transition elements, where a node is located at both 6-way junction + and one of the triple point corners of core box elements. There are two layouts specific to top and bottom + corner: Top (location = 1); and bottom right (location = 2). + :param location: Location identifier. + :return: Node layout. + """ + nodeLayouts = self._nodeLayoutManager.getNodeLayout6WayTriplePoint() + location = abs(location) + # assert location in [1, 2] + + if location == 1: # "Top1" + nodeLayout = nodeLayouts[0] + elif location == 3: # "Top1" + nodeLayout = nodeLayouts[1] + elif location == 2: # "Bottom1" + nodeLayout = nodeLayouts[2] + elif location == 4: # "Bottom2" + nodeLayout = nodeLayouts[3] + + return nodeLayout + + def getNodeLayoutBifurcation(self): + """ + Special node layout for generating core elements for bifurcation. + """ + return self._nodeLayoutBifrucation + + def getNodeLayoutTrifurcation(self, location): + """ + Special node layout for generating core elements for trifurcation. There are two layouts specific to + left-hand side and right-hand side of the solid core cross-section: LHS (location = 1); and RHS (location = 2). + :param location: Location identifier. + :return: Node layout. + """ + if location == 1: # Left-hand side + limitDirections = [None, [[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]], None] + elif location == 2: # Right-hand side + limitDirections = [None, [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], None] + + self._nodeLayoutTrifurcation = self._nodeLayoutManager.getNodeLayout6Way12(True, limitDirections) + return self._nodeLayoutTrifurcation + + def getNodeLayoutTransition(self): + """ + Node layout for generating core transition elements, excluding at triple points. + """ + return self._nodeLayoutTransition + + def getNodeLayoutBifurcationTransition(self): + """ + Special node layout for generating core transition elements for bifurcation. + """ + return self._nodeLayoutBifurcationTransition + + def getNodeLayoutTransitionTriplePoint(self, location): + """ + Special node layout for generating core transition elements at triple points. + There are four layouts specific to each corner of the core box: Top left (location = 1); + top right (location = -1); bottom left (location = 2); and bottom right (location = -2). + :param location: Location identifier identifying four corners of solid core box. + :return: Node layout. + """ + nodeLayouts = self._nodeLayoutManager.getNodeLayoutTriplePoint() + assert location in [1, -1, 2, -2, 0] + if location == 1: # "Top Left" + nodeLayout = nodeLayouts[0] + elif location == -1: # "Top Right" + nodeLayout = nodeLayouts[1] + elif location == 2: # "Bottom Left" + nodeLayout = nodeLayouts[2] + elif location == -2: # "Bottom Right" + nodeLayout = nodeLayouts[3] + else: + nodeLayout = self._nodeLayoutTransition + + self._nodeLayoutTransitionTriplePoint = nodeLayout + return self._nodeLayoutTransitionTriplePoint + def getNodetemplate(self): return self._nodetemplate @@ -87,21 +176,34 @@ def isShowTrimSurfaces(self): class TubeNetworkMeshSegment(NetworkMeshSegment): - def __init__(self, networkSegment, pathParametersList, elementsCountAround, elementsCountThroughWall): + def __init__(self, networkSegment, pathParametersList, elementsCountAround, elementsCountThroughWall, + isCore=False, elementsCountAcrossMajor: int = 4, elementsCountAcrossMinor: int = 4, + elementsCountTransition: int = 1): """ :param networkSegment: NetworkSegment this is built from. :param pathParametersList: [pathParameters] if 2-D or [outerPathParameters, innerPathParameters] if 3-D :param elementsCountAround: Number of elements around this segment. :param elementsCountThroughWall: Number of elements between inner and outer tube if 3-D, 1 if 2-D. + :param isCore: True for generating a solid core inside the tube, False for regular tube network. + :param elementsCountAcrossMajor: Number of elements across major axis of an ellipse. + :param elementsCountAcrossMinor: Number of elements across minor axis of an ellipse. + :param elementsCountTranstion: Number of elements across transition zone between core box elements and + rim elements. """ super(TubeNetworkMeshSegment, self).__init__(networkSegment, pathParametersList) + self._isCore = isCore self._elementsCountAround = elementsCountAround + self._elementsCountAcrossMajor = elementsCountAcrossMajor + self._elementsCountAcrossMinor = elementsCountAcrossMinor + self._elementsCountTransition = elementsCountTransition + if self._isCore and self._elementsCountTransition > 1: + self._elementsCountAround = (elementsCountAround - 8 * (self._elementsCountTransition - 1)) assert elementsCountThroughWall > 0 self._elementsCountThroughWall = elementsCountThroughWall self._rawTubeCoordinatesList = [] self._rawTrackSurfaceList = [] for pathParameters in pathParametersList: - px, pd1, pd2, pd12 = getPathRawTubeCoordinates(pathParameters, elementsCountAround) + px, pd1, pd2, pd12 = getPathRawTubeCoordinates(pathParameters, self._elementsCountAround) self._rawTubeCoordinatesList.append((px, pd1, pd2, pd12)) nx, nd1, nd2, nd12 = [], [], [], [] for i in range(len(px)): @@ -116,12 +218,33 @@ def __init__(self, networkSegment, pathParametersList, elementsCountAround, elem self._rimNodeIds = None self._rimElementIds = None # [e2][e3][e1] + self._boxCoordinates = None + self._transitionCoordinates = None + self._boxNodeIds = None # [nAlong][nAcrossMajor][nAcrossMinor] + self._boxBoundaryNodeIds = None + # boxNodeIds that form the boundary of the solid core, rearranged in circular format + self._boxBoundaryNodeToBoxId = None + # lookup table that translates box boundary node ids in a circular format to box node ids in + # [nAlong][nAcrossMajor][nAcrossMinor] format. + def getElementsCountAround(self): return self._elementsCountAround def getRawTubeCoordinates(self, pathIndex=0): return self._rawTubeCoordinatesList[pathIndex] + def getIsCore(self): + return self._isCore + + def getElementsCountAcrossMajor(self): + return self._elementsCountAcrossMajor + + def getElementsCountAcrossMinor(self): + return self._elementsCountAcrossMinor + + def getElementsCountAcrossTransition(self): + return self._elementsCountTransition + def getRawTrackSurface(self, pathIndex=0): return self._rawTrackSurfaceList[pathIndex] @@ -152,10 +275,12 @@ def sample(self, targetElementLength): ix, id1, id2 = self._sampledTubeCoordinates[1][0:3] rx, rd1, rd2, rd3 = [], [], [], [] for n2 in range(elementsCountAlong + 1): + coreCentre, arcCentre = self._determineCentrePoints(n2) for r in (rx, rd1, rd2, rd3): r.append([]) otx, otd1, otd2 = ox[n2], od1[n2], od2[n2] itx, itd1, itd2 = ix[n2], id1[n2], id2[n2] + # wx, wd3 = self._determineWallCoordinates(otx, otd1, otd2, itx, itd1, itd2, coreCentre, arcCentre) wd3 = [mult(sub(otx[n1], itx[n1]), wallFactor) for n1 in range(self._elementsCountAround)] for n3 in range(self._elementsCountThroughWall + 1): oFactor = n3 / self._elementsCountThroughWall @@ -169,17 +294,651 @@ def sample(self, targetElementLength): x, d1, d2 = otx[n1], otd1[n1], otd2[n1] else: x = add(mult(itx[n1], iFactor), mult(otx[n1], oFactor)) + # x = wx[n3][n1] d1 = add(mult(itd1[n1], iFactor), mult(otd1[n1], oFactor)) d2 = add(mult(itd2[n1], iFactor), mult(otd2[n1], oFactor)) d3 = wd3[n1] + # d3 = wd3[n3][n1] for r, value in zip((rx, rd1, rd2, rd3), (x, d1, d2, d3)): r[n2][n3].append(value) self._rimCoordinates = rx, rd1, rd2, rd3 self._rimNodeIds = [None] * (elementsCountAlong + 1) self._rimElementIds = [None] * elementsCountAlong + if self._isCore: + # sample coordinates for the solid core + self._sampleCoreCoordinates(elementsCountAlong) + + def _sampleCoreCoordinates(self, elementsCountAlong): + """ + Black box function for sampling coordinates for the solid core. + :param elementsCountAlong: A number of elements along a segment. + """ + boxx, boxd1, boxd3 = [], [], [] + transx, transd1, transd3 = [], [], [] + for n2 in range(elementsCountAlong + 1): + coreCentre, arcCentre = self._determineCentrePoints(n2) + cbx, cbd1, cbd3, ctx, ctd1, ctd3 = self._generateCoreCoordinates(n2, coreCentre) + for lst, value in zip((boxx, boxd1, boxd3, transx, transd1, transd3), + (cbx, cbd1, cbd3, ctx, ctd1, ctd3)): + lst.append(value) + boxd2, transd2 = self._determineCoreD2Derivatives(boxx, boxd1, boxd3, transx, transd1, transd3) + self._boxCoordinates = boxx, boxd1, boxd2, boxd3 + self._transitionCoordinates = transx, transd1, transd2, transd3 + self._boxNodeIds = [None] * (elementsCountAlong + 1) + + def _determineCentrePoints(self, n2): + """ + Calculates coordinates for the centre of the solid core based on outer and inner tube coordinates. + :param n2: Index for elements along the tube. + :return: Coordinates of the solid core. + """ + ox = self._sampledTubeCoordinates[0][0][n2] + ix = self._sampledTubeCoordinates[1][0][n2] + cp = [] + + for x in [ox, ix]: + P0 = x[self._elementsCountAround // 4] + P1 = x[self._elementsCountAround // 4 * 3] + midpoint = mult(add(P0, P1), 0.5) + cp.append(midpoint) + + coreCentres, arcCentres = [], [] + for i in range(self._elementsCountAround): + tol = 1e-10 # tolerance to avoid float division zero + OP = magnitude(sub(ox[i], cp[0])) + IP = magnitude(sub(ix[i], cp[1])) + distBetweenOuterAndInner = magnitude(sub(cp[1], cp[0])) + distBetweenOuterAndInner = tol if distBetweenOuterAndInner == 0 else distBetweenOuterAndInner + outerBase = (OP ** 2 - IP ** 2 - distBetweenOuterAndInner ** 2) / (2 * distBetweenOuterAndInner) + circularArcRadius = math.sqrt(outerBase ** 2 + OP ** 2) + distBetweenCoreAndInner = circularArcRadius - outerBase - distBetweenOuterAndInner + + directionVector = sub(cp[1], cp[0]) + if directionVector[0] == 0 and directionVector[1] == 0 and directionVector[2] == 0: + directionVector = [tol, tol, tol] + scaledDV = mult(normalize(directionVector), (distBetweenOuterAndInner + distBetweenCoreAndInner)) + c = add(scaledDV, cp[0]) + coreCentres.append(c) + + dvi = [-d for d in directionVector] + mag = magnitude(dvi) + tol = 1e-10 + scaleDVI = mult(normalize(dvi), outerBase) if mag > tol else [0, 0, 0] + ac = add(scaleDVI, cp[0]) + arcCentres.append(ac) + + coreCentre = [sum(e) / len(e) for e in zip(*coreCentres)] + arcCentre = [sum(e) / len(e) for e in zip(*arcCentres)] + + return coreCentre, arcCentre + + def _createMirrorCurve(self, n2, centre): + """ + Generate coordinates and derivatives for the mirror curve. + :param n2: Index for elements along the tube. + :param centre: Coordinates of the solid core. + :return: Coordinates and derivatives for the mirror curve. + """ + ix = self._rimCoordinates[0][n2][0] + id3 = self._rimCoordinates[3][n2][0] + + # Create mirror curves + n2a = 0 + n2z = self._elementsCountAround // 2 + + rscx, rscd1 = [], [] + for n in range(2): + startIndex = [n2a, n2z][n] + tmdxStart = ix[startIndex] if n == 0 else centre + tmdxEnd = centre if n == 0 else ix[startIndex] + tmdd3 = id3[startIndex] + rcx = [tmdxStart, tmdxEnd] + mag = -1 if n == 0 else 1 + rcd3 = [set_magnitude(tmdd3, mag), set_magnitude(tmdd3, mag)] + elementsCountAcross = self._elementsCountAcrossMajor // 2 + + tx, td1 = sampleCubicHermiteCurves(rcx, rcd3, elementsCountAcross, arcLengthDerivatives=True)[0:2] + + rscx += tx[n::] + rscd1 += td1[n::] + + rscd1 = smoothCubicHermiteDerivativesLine(rscx, rscd1, fixStartDirection=True, fixEndDirection=True) + + # Determine d3 derivatives + rscd3 = [] + for n in range(len(rscx)): + d3 = normalize( + [ix[self._elementsCountAround // 4][c] - ix[(self._elementsCountAround // 4) * 3][c] for c in range(3)]) + rscd3.append(d3) + + return rscx, rscd1, rscd3 + + def _sampleCoreNodesAlongMinorAxis(self, n2, rscx, rscd1, rscd3): + """ + Samples nodes along minor axis (z-axis) between inner tube coordinate, centre and the opposing inner tube + coordinate by the number of elements across minor axis. + :param n2: Index for elements along the tube. + :param rscx, rscd1, rscd3: Lists of coordinates, d1 and d3 derivatives for mirror curve. + :return: Coordinates and derivatives for the solid core and the transition nodes around the core. + """ + ix = self._rimCoordinates[0][n2][0] + id1 = self._rimCoordinates[1][n2][0] + id3 = self._rimCoordinates[3][n2][0] + + # Create an empty list for the core with dimensions of (major - 1) x (minor - 1) + m = self._elementsCountAcrossMajor - 1 - 2 * (self._elementsCountTransition - 1) + n = self._elementsCountAcrossMinor - 1 - 2 * (self._elementsCountTransition - 1) + cbx, cbd1, cd2, cbd3, ctx, ctd1, ctd2, ctd3 = [], [], [], [], [], [], [], [] + + for i in range(m): + for lst in (cbx, cbd1, cd2, cbd3): + lst.append([[0, 0, 0] for _ in range(n)]) + + for i in range(self._elementsCountTransition - 1): + for lst in (ctx, ctd1, ctd2, ctd3): + lst.append([[0, 0, 0] for _ in range(self._elementsCountAround)]) + + # Create Regular Row Curves + n2a = 0 + n2b = self._elementsCountTransition - 1 + n2d = n1d = self._elementsCountAcrossMinor // 2 - n2b + n2m = self._elementsCountAcrossMajor - (4 + 2 * n2b) + n2d + n2z = self._elementsCountAcrossMinor + + ixTopHalf = ix[0:(self._elementsCountAround // 2) + 1] + id1TopHalf = id1[0:(self._elementsCountAround // 2) + 1] + id3TopHalf = id3[0:(self._elementsCountAround // 2) + 1] + ixBtmHalf = (ix[(self._elementsCountAround // 2)::])[::-1] + id1BtmHalf = (id1[(self._elementsCountAround // 2)::])[::-1] + id3BtmHalf = (id3[(self._elementsCountAround // 2)::])[::-1] + + tx, td1, td3 = [], [], [] + trx, trd1, trd3 = [], [], [] + for n in range(n2m + 1 - n2d): + n2 = list(range(n2d, n2m + 1))[n] + c = self._elementsCountTransition + 1 + c2 = list(range(c, c + (n2m + 1 - n2d)))[n] + nd1 = [set_magnitude(id3BtmHalf[n2 - 1], -1.0), set_magnitude(rscd3[c2], 1.0), + set_magnitude(id3TopHalf[n2], 1.0)] + txm, td3m, pe, pxi, psf = sampleCubicHermiteCurves([ixBtmHalf[n2 - 1], rscx[c2], ixTopHalf[n2]], nd1, + self._elementsCountAcrossMinor, + arcLengthDerivatives=True) + td1m = interpolateSampleCubicHermite([[-id1BtmHalf[n2 - 1][c] for c in range(3)], rscd1[c2], + id1TopHalf[n2]], [[0.0, 0.0, 0.0]] * 3, pe, pxi, psf)[0] + + td3m = smoothCubicHermiteDerivativesLine(txm, td3m, fixStartDirection=True, fixEndDirection=True) + + [lst1.append(lst2[n2a + 1: n2z]) for lst1, lst2 in zip((tx, td1, td3), (txm, td1m, td3m))] + + if self._elementsCountTransition > 1: + for i in range(len(tx)): + [lst.append([]) for lst in (trx, trd1, trd3)] + for j in range(self._elementsCountTransition - 1): + trx[-1].append([tx[i].pop(0), tx[i].pop(-1)]) + trd1[-1].append([td1[i].pop(0), td1[i].pop(-1)]) + trd3[-1].append([td3[i].pop(0), td3[i].pop(-1)]) + + # Store coordinate and derivative values into appropriate lists + for n2c in range(1, self._elementsCountAcrossMajor - 2 - 2 * (self._elementsCountTransition - 1)): + for n1 in range(len(tx[n2c - 1])): + for lst, value in zip([cbx[n2c], cbd1[n2c], cbd3[n2c]], + [tx[n2c - 1][n1], td1[n2c - 1][n1], td3[n2c - 1][n1]]): + lst[n1] = value + + if self._elementsCountTransition > 1: + for n2 in range(self._elementsCountTransition - 1): + for n in range(len(trx)): + n1 = -(n1d + n) + for lst, value in zip([ctx[n2], ctd1[n2], ctd3[n2]], + [trx[n][n2][0], trd1[n][n2][0], trd3[n][n2][0]]): + lst[n1] = value + for n in range(len(trx)): + n1 = n1d + n + for lst, value in zip([ctx[n2], ctd1[n2], ctd3[n2]], + [trx[n][n2][-1], trd1[n][n2][-1], trd3[n][n2][-1]]): + lst[n1] = value + + return cbx, cbd1, cbd3, ctx, ctd1, ctd3 + + def _sampleCoreNodeAlongMajorAxis(self, n2, cbx, cbd1, cbd3, ctx, ctd1, ctd3): + """ + Samples nodes along major axis (y-axis) between inner tube coordinate, centre and the opposing inner tube + coordinate by the number of elements across major axis. + :param n2: Index for elements along the tube. + :param cbx, cbd1, cbd3: Lists of coordinates, d1 and d3 derivatives for solid core. + :param ctx, ctd1, ctd3: Lists of coordinates, d1 and d3 derivatives for rims around the core. + :return: Coordinates and derivatives for the solid core and the rim(s) around the core. + """ + ix, id1, id3 = self._rimCoordinates[0][n2][0], self._rimCoordinates[1][n2][0], self._rimCoordinates[3][n2][0] + + # Create regular column curves + n1a, n1m, n1z = 0, self._elementsCountAround // 2, self._elementsCountAcrossMajor + ec = (self._elementsCountAcrossMinor - 4) // 2 - (self._elementsCountTransition - 1) + startIndexes = list(range(n1a - ec, n1a + ec + 1)) + endIndexes = list(range(n1m - ec, n1m + ec + 1))[::-1] + tx, td1, td3, trx, trd1, trd3 = [], [], [], [], [], [] + + elementsCountAcross = self._elementsCountAcrossMajor // 2 + nloop = (self._elementsCountAcrossMinor - 3 - 2 * (self._elementsCountTransition - 1)) \ + if self._elementsCountAcrossMinor > 4 else self._elementsCountAcrossMinor - 3 # This is affected by ellipse parameters issue + + for n in range(nloop): + n1, n2 = startIndexes[n], endIndexes[n] + m1 = n + 1 + m2 = (self._elementsCountAcrossMajor - 1 - 2 * (self._elementsCountTransition - 1)) // 2 + txm, td1m, td3m = [], [], [] + for n in range(2): + if n == 0: + startx, startd1 = ix[n1], id1[n1] + endx, endd1, endd3 = cbx[m2][m1], cbd1[m2][m1], cbd3[m2][m1] + startd3 = [-id3[n1][c] for c in range(3)] + nd1 = [set_magnitude(startd3, 1), set_magnitude(endd1, 1)] + nd3 = [startd1, endd3] + else: + startx, startd1, startd3 = cbx[m2][m1], cbd1[m2][m1], cbd3[m2][m1] + endx, endd3 = ix[n2], id3[n2] + endd1 = [-id1[n2][c] for c in range(3)] + nd1 = [set_magnitude(startd1, 1), set_magnitude(endd3, 1)] + nd3 = [startd3, endd1] + + tempx, tempd1, pe, pxi, psf = sampleCubicHermiteCurves([startx, endx], nd1, + elementsCountAcross, arcLengthDerivatives=True) + + tempd3 = interpolateSampleCubicHermite(nd3, [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + + for lst, value in zip([txm, td1m, td3m], [tempx, tempd1, tempd3]): + lst += value[n::] + + td1m = smoothCubicHermiteDerivativesLine(txm, td1m, fixStartDirection=True, fixEndDirection=True) + + for lst, value in zip([tx, td1, td3], [txm, td1m, td3m]): + lst.append(value[n1a + 1: n1z]) + + if self._elementsCountTransition > 1: + for i in range(len(tx)): + [lst.append([]) for lst in (trx, trd1, trd3)] + for j in range(self._elementsCountTransition - 1): + for lst, value in zip([trx[-1], trd1[-1], trd3[-1]], [tx[i], td1[i], td3[i]]): + lst.append([value.pop(0), value.pop(-1)]) + + for n2c in range(len(tx[0])): + for n1 in range(1, self._elementsCountAcrossMinor - 2 - 2 * (self._elementsCountTransition - 1)): + for lst, value in zip([cbx[n2c], cbd1[n2c], cbd3[n2c]], + [tx[n1 - 1][n2c], td1[n1 - 1][n2c], td3[n1 - 1][n2c]]): + lst[n1] = value + + if self._elementsCountTransition > 1: + for n2 in range(self._elementsCountTransition - 1): + for n in range(len(trx)): + n1 = startIndexes[n] + for lst, value in zip([ctx[n2], ctd1[n2], ctd3[n2]], + [trx[n][n2][0], trd1[n][n2][0], trd3[n][n2][0]]): + lst[n1] = value + for n in range(len(trx)): + n1 = endIndexes[n] + for lst, value in zip([ctx[n2], ctd1[n2], ctd3[n2]], + [trx[n][n2][-1], trd1[n][n2][-1], trd3[n][n2][-1]]): + lst[n1] = value + + # Fix derivatives for core rim + nSkip = self._elementsCountAcrossMinor // 2 - (1 + self._elementsCountTransition) + if ctx: + for n2 in range(self._elementsCountTransition - 1): + for n1 in range(-nSkip, nSkip + 1): + tempd1, tempd3 = copy.copy(ctd1[n2][n1]), copy.copy(ctd3[n2][n1]) + tempd1 = [-tempd1[c] for c in range(3)] + ctd1[n2][n1], ctd3[n2][n1] = tempd3, tempd1 + for n1 in range(self._elementsCountAround // 2 - nSkip, self._elementsCountAround // 2 + nSkip + 1): + tempd1, tempd3 = copy.copy(ctd1[n2][n1]), copy.copy(ctd3[n2][n1]) + tempd3 = [-tempd3[c] for c in range(3)] + ctd1[n2][n1], ctd3[n2][n1] = tempd3, tempd1 + for n1 in range(self._elementsCountAround // 4 * 3 - nSkip, + self._elementsCountAround // 4 * 3 + nSkip + 1): + tempd1, tempd3 = copy.copy(ctd1[n2][n1]), copy.copy(ctd3[n2][n1]) + tempd1, tempd3 = [-tempd1[c] for c in range(3)], [-tempd3[c] for c in range(3)] + ctd3[n2][n1], ctd1[n2][n1] = tempd3, tempd1 + + # Re-order the order of sublists, i.e. from inner to outer layer. + ctx = ctx[::-1] + + return cbx, cbd1, cbd3, ctx, ctd1, ctd3 + + def _determineCoreTriplePoints(self, n2, cbx, cbd1, cbd3, ctx, ctd1, ctd3): + """ + Compute coordinates and derivatives of points where 3 square elements merge. + :param n2: Index for elements along the tube. + :param cbx, cbd1, cbd3: Coordinates, d1 and d3 derivatives of core box. + :param ctx, ctd1, ctd3: Coordinates, d1 and d3 derivatives of core rim(s). + :return: Coordinates and derivatives of box and rim components of the core. + """ + ix, id1, id3 = self._rimCoordinates[0][n2][0], self._rimCoordinates[1][n2][0], self._rimCoordinates[3][n2][0] + + n1m = self._elementsCountAround // 2 + n2a = 0 + cix, cid1, cid3 = (ctx + [ix], ctd1 + [id1], ctd3 + [id3]) + minorSkip = self._elementsCountAcrossMinor // 2 - (1 + self._elementsCountTransition) + majorSkip = self._elementsCountAcrossMajor // 2 - (1 + self._elementsCountTransition) + + # Generate triple points for the core box + # Left top and bottom + for n in range(2): + ltx = [] + scalefactor = 1 if n == 0 else -1 + n1 = [-minorSkip, n1m + minorSkip][n] + n1a = 1 + n1z = self._elementsCountAround // 4 * 3 + (majorSkip * scalefactor) + n2 = [n2a, -1][n] + n2b = [1, -2][n] + tx, td1 = sampleCubicHermiteCurves([cix[0][n1], cbx[n2b][0]], + [[(-cid1[0][n1][c] * scalefactor - cid3[0][n1][c]) for c in range(3)], + set_magnitude(cbd1[n2b][0], scalefactor)], 2, arcLengthDerivatives=True)[ + 0:2] + ltx.append(tx[1]) + tx, td1 = sampleCubicHermiteCurves([cix[0][n1z], cbx[n2][n1a]], + [[(cid1[0][n1z][c] * scalefactor - cid3[0][n1z][c]) for c in range(3)], + cbd3[n2][n1a]], + 2, arcLengthDerivatives=True)[0:2] + ltx.append(tx[1]) + x = [(ltx[0][c] + ltx[1][c]) / 2.0 for c in range(3)] + cbx[n2][0] = x + cbd3[n2][0] = [(cbx[n2][1][c] - cbx[n2][0][c]) for c in range(3)] + cbd1[n2][0] = [(cbx[n2 + scalefactor][0][c] - cbx[n2][0][c]) * scalefactor for c in range(3)] + # Right + for n in range(2): + rtx = [] + scalefactor = 1 if n == 0 else -1 + n1 = [minorSkip, n1m - minorSkip][n] + n1b = -2 + n2 = [n2a, -1][n] + n2b = [1, -2][n] + n1z = self._elementsCountAround // 4 - (majorSkip * scalefactor) + tx, td1 = sampleCubicHermiteCurves([cix[0][n1], cbx[n2b][-1]], + [[(cid1[0][n1][c] * scalefactor - cid3[0][n1][c]) for c in range(3)], + [d * scalefactor for d in cbd1[n2b][-1]]], 2, + arcLengthDerivatives=True)[0:2] + rtx.append(tx[1]) + tx, td1 = sampleCubicHermiteCurves([cix[0][n1z], cbx[n2][n1b]], + [[(-cid1[0][n1z][c] * scalefactor - cid3[0][n1z][c]) for c in range(3)], + [-d for d in cbd3[n2][n1b]]], 2, arcLengthDerivatives=True)[0:2] + rtx.append(tx[1]) + x = [(rtx[0][c] + rtx[1][c]) / 2.0 for c in range(3)] + n2c = [0, - 1][n] + cbx[n2][-1] = x + cbd3[n2][-1] = [(cbx[n2c][-1][c] - cbx[n2c][-2][c]) for c in range(3)] + cbd1[n2][-1] = [(cbx[n2b][-1][c] - cbx[n2c][-1][c]) * scalefactor for c in range(3)] + + # Sample nodes along triple points + if self._elementsCountTransition > 1: + # Left + for n in range(2): + scalefactor = 1 if n == 0 else -1 + n1a = [-minorSkip, n1m + minorSkip][n] + n1c = [n1a - 1, n1a + 1][n] + n2 = [0, -1][n] + txm, td3m, pe, pxi, psf = sampleCubicHermiteCurves([ix[n1c], cbx[n2][0]], + [[-id3[n1c][c] for c in range(3)], + [cbd1[n2][0][c] * scalefactor + cbd3[n2][0][c] for c + in range(3)]], + self._elementsCountTransition, + arcLengthDerivatives=True) + td1m = interpolateSampleCubicHermite( + [id1[n1c], ([-cbd1[n2][0][c] + cbd3[n2][0][c] * scalefactor for c in range(3)])], + [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + + txm, td1m, td3m = txm[::-1], td1m[::-1], td3m[::-1] + for n3 in range(self._elementsCountTransition - 1): + ctx[n3][n1c] = txm[n3 + 1] + ctd1[n3][n1c] = td1m[n3 + 1] + ctd3[n3][n1c] = [-d for d in td3m[n3 + 1]] + # Right + for n in range(2): + scalefactor = 1 if n == 0 else -1 + n1a = [minorSkip, n1m - minorSkip][n] + n1c = [n1a + 1, n1a - 1][n] + n2 = [0, -1][n] + txm, td3m, pe, pxi, psf = sampleCubicHermiteCurves([ix[n1c], cbx[n2][-1]], + [[-id3[n1c][c] for c in range(3)], + [cbd1[n2][-1][c] * scalefactor - cbd3[n2][-1][c] for + c in range(3)]], + self._elementsCountTransition, + arcLengthDerivatives=True) + td1m = interpolateSampleCubicHermite([id1[n1c], + ([cbd1[n2][-1][c] + cbd3[n2][-1][c] * scalefactor for c in + range(3)])], + [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + + txm, td1m, td3m = txm[::-1], td1m[::-1], td3m[::-1] + for n3 in range(self._elementsCountTransition - 1): + ctx[n3][n1c] = txm[n3 + 1] + ctd1[n3][n1c] = td1m[n3 + 1] + ctd3[n3][n1c] = [-d for d in td3m[n3 + 1]] + + return cbx, cbd1, cbd3, ctx, ctd1, ctd3 + + def _smoothCoreDerivatives(self, cbx, cbd1, cbd3, ctx, ctd1, ctd3): + """ + Smooth d1 and d3 derivatives of the solid core. + :param cbx, cbd1, cbd3: Coordinates, d1 and d3 derivatives of the core box. + :param ctx, ctd1, ctd3: Coordinates, d1 and d3 derivatives of the core rim. + :return: Smoothed d1 and d3 derivatives of the solid core. + """ + # Smooth core box rows + for n2c in [0, -1]: + cbd3[n2c][1:-1] = smoothCubicHermiteDerivativesLine(cbx[n2c], cbd3[n2c])[1:-1] + + # Smooth core box columns + for n1 in [0, -1]: + tx = [] + td1 = [] + for n2 in range(len(cbx)): + tx.append(cbx[n2][n1]) + td1.append(cbd1[n2][n1]) + td1 = smoothCubicHermiteDerivativesLine(tx, td1) + for n2 in range(1, len(td1)): + cbd1[n2][n1] = td1[n2] + + # Smooth core transition derivatives + if self._elementsCountTransition > 1: + for m in range(self._elementsCountTransition - 1): + ctxTopHalf = ctx[m][0:(self._elementsCountAround // 2) + 1] + ctd1TopHalf = ctd1[m][0:(self._elementsCountAround // 2) + 1] + ctxBtmHalf = (ctx[m][(self._elementsCountAround // 2)::]) + ctd1BtmHalf = (ctd1[m][(self._elementsCountAround // 2)::]) + ctxHalf = [ctxTopHalf, ctxBtmHalf] + ctd1Half = [ctd1TopHalf, ctd1BtmHalf] + + for n in range(2): + ctd1m = smoothCubicHermiteDerivativesLine(ctxHalf[n], ctd1Half[n], + fixStartDirection=True, fixEndDirection=True) + if n == 0: + ctd1[m][0:(self._elementsCountAround // 2) + 1] = ctd1m + else: + ctd1[m][(self._elementsCountAround // 2)::] = ctd1m + + return cbd1, cbd3, ctd1, ctd3 + + def _generateCoreCoordinates(self, n2, coreCentre): + """ + Creates coordinates and derivatives for the solid core (and the rim(s) around the core, if exists) based on + outer and inner tube coordinates. + :param n2: Index for elements along the tube. + :param coreCentre: Coordinates at the centre point of solid core + :return: Lists of coordinates and derivatives for the solid core and rims around the core. + """ + # Determine mirror curves + rscx, rscd1, rscd3 = self._createMirrorCurve(n2, coreCentre) + # Create regular row curves + cbx, cbd1, cbd3, ctx, ctd1, ctd3 = self._sampleCoreNodesAlongMinorAxis(n2, rscx, rscd1, rscd3) + # Create regular column curves + cbx, cbd1, cbd3, ctx, ctd1, ctd3 = self._sampleCoreNodeAlongMajorAxis(n2, cbx, cbd1, cbd3, ctx, ctd1, ctd3) + # Get triple points + cbx, cbd1, cbd3, ctx, ctd1, ctd3 = self._determineCoreTriplePoints(n2, cbx, cbd1, cbd3, ctx, ctd1, ctd3) + # Smooth derivatives + cbd1, cbd3, ctd1, ctd3 = self._smoothCoreDerivatives(cbx, cbd1, cbd3, ctx, ctd1, ctd3) + + return cbx, cbd1, cbd3, ctx, ctd1, ctd3 + + def _determineCoreD2Derivatives(self, boxx, boxd1, boxd3, transx, transd1, transd3): + """ + Compute d2 derivatives for the solid core. + :param boxx, boxd1, boxd3: Coordinates and derivatives (d1 & d3) of the core box nodes. + :param transx, transd1, transd3: Coordinates and derivatives (d1 & d3) of the core transition nodes. + :return: D2 derivatives of box and rim components of the core. + """ + elementsCountAlong = len(boxx) + nodesCountAcrossMajor = len(boxx[0]) + nodesCountAcrossMinor = len(boxx[0][0]) + + boxd2 = [[[None for _ in range(nodesCountAcrossMinor)] for _ in range(nodesCountAcrossMajor)] + for _ in range(elementsCountAlong)] + transd2 = [[[None for _ in range(self._elementsCountAround)] for _ in range(self._elementsCountTransition - 1)] + for _ in range(elementsCountAlong)] + + for m in range(nodesCountAcrossMajor): + for n in range(nodesCountAcrossMinor): + tx, td2 = [], [] + for n2 in range(elementsCountAlong): + x = boxx[n2][m][n] + d2 = cross(boxd3[n2][m][n], boxd1[n2][m][n]) + tx.append(x) + td2.append(d2) + td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixStartDirection=False, fixEndDirection=False) + for n2 in range(elementsCountAlong): + boxd2[n2][m][n] = td2[n2] + + if self._elementsCountTransition > 1: + for n3 in range(self._elementsCountTransition - 1): + for n1 in range(self._elementsCountAround): + tx, td2 = [], [] + for n2 in range(elementsCountAlong): + x = transx[n2][n3][n1] + d2 = cross(transd3[n2][n3][n1], transd1[n2][n3][n1]) + tx.append(x) + td2.append(d2) + td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixStartDirection=False, fixEndDirection=True) + for n2 in range(elementsCountAlong): + transd2[n2][n3][n1] = td2[n2] + + return boxd2, transd2 + + def _determineWallCoordinates(self, ox, od1, od2, ix, id1, id2, coreCentre, arcCentre): + """ + Calculates rim coordinates and d3 derivatives based on the centre point of the solid core. + :param ox, od1, od2: Coordinates and (d1 and d2) derivatives for outermost rim. + :param ix, id1, id2: Coordinates and (d1 and d2) derivatives for innermost rim. + :param coreCentre: Centre point of the solid core. + :param arcCetnre: Centre point of the arc that passes through the core centre, inner rim and outer rim. + :return: Coordinates and d3 derivatives for rim nodes. + """ + wx, wd3 = [], [] + + # check if the cross-section of cylinder is regular shaped or irregular. + dist1a = sub(ix[0], coreCentre) + dist1b = sub(ix[self._elementsCountAround // 2], coreCentre) + dist2a = sub(ix[self._elementsCountAround // 4], coreCentre) + dist2b = sub(ix[self._elementsCountAround // 4 * 3], coreCentre) + tol = 1e-3 + if abs(magnitude(dist1a) - magnitude(dist1b)) > tol or \ + abs(magnitude(dist2a) - magnitude(dist2b)) > tol: + isRegular = False + else: + isRegular = True + + # Calculate d3 derivatives + tx, td3 = [], [] + for n1 in range(self._elementsCountAround): + if isRegular: + tol = 1e-10 + dist = sub(arcCentre, coreCentre) + if magnitude(dist) > tol: + if dist > [tol, tol, tol]: + oc = sub(ox[n1], arcCentre) + ic = sub(ix[n1], arcCentre) + else: + oc = add(mult(oc[n1], -1), arcCentre) + ic = add(mult(ic[n1], -1), arcCentre) + ot = cross(oc, od1[n1]) + it = cross(ic, id1[n1]) + else: + ot, it = cross(od1[n1], od2[n1]), cross(id1[n1], id2[n1]) + scalefactor = magnitude(sub(ox[n1], ix[n1])) / self._elementsCountThroughWall + od3 = mult(normalize(ot), scalefactor) + id3 = mult(normalize(it), scalefactor) + else: + wallFactor = 1.0 / self._elementsCountThroughWall + od3 = id3 = mult(sub(ox[n1], ix[n1]), wallFactor) + + txm, td3m, pe, pxi, psf = sampleCubicHermiteCurves( + [ix[n1], ox[n1]], [id3, od3], self._elementsCountThroughWall, arcLengthDerivatives=True) + + td3m = smoothCubicHermiteDerivativesLine(txm, td3m, fixStartDirection=True, fixEndDirection=True) + + tx.append(txm) + td3.append(td3m) + + for n3 in range(self._elementsCountThroughWall + 1): + wx.append([]) + wd3.append([]) + for n1 in range(self._elementsCountAround): + wx[n3].append(tx[n1][n3]) + wd3[n3].append(td3[n1][n3]) + + return wx, wd3 + + def _createBoxBoundaryNodeIdsList(self, startSkipCount=None, endSkipCount=None): + """ + Creates a list (in a circular format similar to other rim node id lists) of core box node ids that are + located at the boundary of the core. + This list is used to easily stitch inner rim nodes with box nodes. + :param startSkipCount: Row in from start that node ids are for. + :param endSkipCount: Row in from end that node ids are for. + :return: A list of box node ids stored in a circular format, and a lookup list that translates indexes used in + boxBoundaryNodeIds list to indexes that can be used in boxCoordinates list. + """ + boxBoundaryNodeIds = [] + boxBoundaryNodeToBoxId = [] + elementsCountAlong = len(self._rimCoordinates[0]) - 1 + + boxElementsCountRow = (self._elementsCountAcrossMajor - 2 * self._elementsCountTransition) + 1 + boxElementsCountColumn = (self._elementsCountAcrossMinor - 2 * self._elementsCountTransition) + 1 + for n2 in range(elementsCountAlong + 1): + if (n2 < startSkipCount) or (n2 > elementsCountAlong - endSkipCount) or self._boxNodeIds[n2] is None: + boxBoundaryNodeIds.append(None) + boxBoundaryNodeToBoxId.append(None) + continue + else: + boxBoundaryNodeIds.append([]) + boxBoundaryNodeToBoxId.append([]) + for n3 in range(boxElementsCountRow): + if n3 == 0 or n3 == boxElementsCountRow - 1: + ids = self._boxNodeIds[n2][n3] if n3 == 0 else self._boxNodeIds[n2][n3][::-1] + n1List = list(range(boxElementsCountColumn)) if n3 == 0 else ( + list(range(boxElementsCountColumn - 1, -1, -1))) + boxBoundaryNodeIds[n2] += [ids[c] for c in range(boxElementsCountColumn)] + for n1 in n1List: + boxBoundaryNodeToBoxId[n2].append([n3, n1]) + else: + for n1 in [-1, 0]: + boxBoundaryNodeIds[n2].append(self._boxNodeIds[n2][n3][n1]) + boxBoundaryNodeToBoxId[n2].append([n3, n1]) + + start = self._elementsCountAcrossMajor - 4 - 2 * (self._elementsCountTransition - 1) + idx = self._elementsCountAcrossMinor - 2 * (self._elementsCountTransition - 1) + for n in range(int(start), -1, -1): + boxBoundaryNodeIds[n2].append(boxBoundaryNodeIds[n2].pop(idx + 2 * n)) + boxBoundaryNodeToBoxId[n2].append(boxBoundaryNodeToBoxId[n2].pop(idx + 2 * n)) + + nloop = self._elementsCountAcrossMinor // 2 - self._elementsCountTransition + for _ in range(nloop): + boxBoundaryNodeIds[n2].insert(len(boxBoundaryNodeIds[n2]), boxBoundaryNodeIds[n2].pop(0)) + boxBoundaryNodeToBoxId[n2].insert(len(boxBoundaryNodeToBoxId[n2]), + boxBoundaryNodeToBoxId[n2].pop(0)) + + return boxBoundaryNodeIds, boxBoundaryNodeToBoxId + @classmethod - def blendSampledCoordinates(cls, segment1 , nodeIndexAlong1, segment2, nodeIndexAlong2): + def blendSampledCoordinates(cls, segment1, nodeIndexAlong1, segment2, nodeIndexAlong2): nodesCountAround = segment1._elementsCountAround nodesCountRim = len(segment1._rimCoordinates[0][0]) if ((nodesCountAround != segment2._elementsCountAround) or @@ -233,6 +992,135 @@ def getRimCoordinatesListAlong(self, n1, n2List, n3): paramsList.append(params) return paramsList + def getBoxCoordinatesListAlong(self, n1, n2List, n3): + """ + Get a list of parameters for solid core box for n2 indexes along segment at given n1, n3. + :param n1: Node index around segment. + :param n2List: List of node indexes along segment. + :param n3: Node index from inner to outer rim. + :return: [x[], d1[], d2[], d3[]]. + """ + paramsList = [] + for i in range(4): + params = [] + for n2 in n2List: + params.append(self._boxCoordinates[i][n2][n3][n1] if self._boxCoordinates[i] else None) + paramsList.append(params) + + return paramsList + + def getTransitionCoordinatesListAlong(self, n1, n2List, n3): + """ + Get a list of parameters for core transition nodes for n2 indexes along segment at given n1, n3. + :param n1: Node index around segment. + :param n2List: List of node indexes along segment. + :param n3: Node index from inner to outer rim. + :return: [x[], d1[], d2[], d3[]]. + """ + paramsList = [] + for i in range(4): + params = [] + for n2 in n2List: + params.append(self._transitionCoordinates[i][n2][n3][n1] if self._transitionCoordinates[i] + else None) + paramsList.append(params) + + return paramsList + + def getCoreNodesCountAcrossMajor(self): + return len(self._boxCoordinates[0][0]) + + def getCoreNodesCountAcrossMinor(self): + return len(self._boxCoordinates[0][0][0]) + + def getElementsCountTransition(self): + return self._elementsCountTransition + + def getBoxCoordinates(self, n1, n2, n3): + return (self._boxCoordinates[0][n2][n3][n1], self._boxCoordinates[1][n2][n3][n1], + self._boxCoordinates[2][n2][n3][n1], self._boxCoordinates[3][n2][n3][n1]) + + def getBoxNodeIds(self, n1, n2, n3): + """ + Get a box node ID for n2 index along segment at given n1, n3. + :param n1: Node index across major axis (y-axis). + :param n2: Node index along segment. + :param n3: Node index across minor axis (z-axis). + :return: Node identifier. + """ + return self._boxNodeIds[n2][n3][n1] + + def getBoxNodeIdsSlice(self, n2): + """ + Get slice of box node IDs at n2 index along segment. + :param n2: Node index along segment, including negative indexes from end. + :return: Node IDs arrays, or None if not set. + """ + return self._boxNodeIds[n2] + + def getBoxBoundaryNodeIds(self, n1, n2): + """ + Get a node ID around the core box for n2 index along segment at a given n1. + :param n1: Node index around the core box. + :param n2: Node index along segment. + :return: Node identifier. + """ + return self._boxBoundaryNodeIds[n2][n1] + + def getBoxBoundaryNodeToBoxId(self, n1, n2): + """ + Translate box boundary node indexes to core box node indexes in the form of major- and minor axes. + :param n1: Node index around the core box. + :param n2: Node index along segment. + :return: n3 (major axis) and n1 (minor axis) indexes used in boxCoordinates. + """ + return self._boxBoundaryNodeToBoxId[n2][n1] + + def getTriplePointIndexes(self): + """ + Get a node ID at triple points (special four corners) of the solid core. + :return: A list of circular (n1) indexes used to identify triple points. + """ + elementsCountAround = self._elementsCountAround + nodesCountAcrossMinorHalf = self.getCoreNodesCountAcrossMinor() // 2 + triplePointIndexesList = [] + + for n in range(0, elementsCountAround, elementsCountAround // 2): + triplePointIndexesList.append(n + nodesCountAcrossMinorHalf) + triplePointIndexesList.append((n - nodesCountAcrossMinorHalf) % elementsCountAround) + + return triplePointIndexesList + + def getTriplePointLocation(self, e1): + """ + Determines the location of a specific triple point relative to the solid core box. + There are four locations: Top left (location = 1); top right (location = -1); bottom left (location = 2); + and bottom right (location = -2). Location is None if not located at any of the four specified locations. + :return: Location identifier. + """ + em = (self._elementsCountAcrossMinor - 2) // 2 - (self._elementsCountTransition - 1) + eM = (self._elementsCountAcrossMajor - 2) // 2 - (self._elementsCountTransition - 1) + ec = self._elementsCountAround // 4 + + lftColumnElements = list(range(0, ec - eM)) + list(range(3 * ec + eM, self._elementsCountAround)) + topRowElements = list(range(ec - eM, ec + eM)) + rhtColumnElements = list((range(2 * ec - em, 2 * ec + em))) + btmRowElements = list(range(3 * ec - eM, 3 * ec + eM)) + + idx = len(lftColumnElements) // 2 + if e1 == topRowElements[0] or e1 == lftColumnElements[idx - 1]: + location = 1 # "TopLeft" + elif e1 == topRowElements[-1] or e1 == rhtColumnElements[0]: + location = -1 # "TopRight" + elif e1 == btmRowElements[-1] or e1 == lftColumnElements[idx]: + location = 2 # "BottomLeft" + elif e1 == btmRowElements[0] or e1 == rhtColumnElements[-1]: + location = -2 # "BottomRight" + else: + location = 0 + + return location + def getRimCoordinates(self, n1, n2, n3): """ Get rim parameters at a point. @@ -293,6 +1181,7 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): """ elementsCountAlong = len(self._rimCoordinates[0]) - 1 elementsCountRim = self.getElementsCountRim() + elementsCountTransition = self.getElementsCountTransition() coordinates = generateData.getCoordinates() fieldcache = generateData.getFieldcache() startSkipCount = 1 if (self._junctions[0].getSegmentsCount() > 2) else 0 @@ -304,15 +1193,25 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): nodetemplate = generateData.getNodetemplate() for n2 in range(elementsCountAlong + 1) if (n2Only is None) else [n2Only]: if (n2 < startSkipCount) or (n2 > elementsCountAlong - endSkipCount): + if self._isCore: + self._boxNodeIds[n2] = None self._rimNodeIds[n2] = None continue - if self._rimNodeIds[n2]: - continue + if self._isCore: + if self._rimNodeIds[n2] and self._boxNodeIds[n2]: + continue + else: + if self._rimNodeIds[n2]: + continue # get shared nodes from single adjacent segment, including loop on itself # only handles one in, one out if n2 == 0: if self._junctions[0].getSegmentsCount() == 2: segments = self._junctions[0].getSegments() + if self._isCore: + boxNodeIds = segments[0].getBoxNodeIdsSlice(-1) + if boxNodeIds: + self._boxNodeIds[n2] = boxNodeIds rimNodeIds = segments[0].getRimNodeIdsSlice(-1) if rimNodeIds: self._rimNodeIds[n2] = rimNodeIds @@ -320,19 +1219,54 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): if n2 == elementsCountAlong: if self._junctions[1].getSegmentsCount() == 2: segments = self._junctions[1].getSegments() + if self._isCore: + boxNodeIds = segments[1].getBoxNodeIdsSlice(0) + if boxNodeIds: + self._boxNodeIds[n2] = boxNodeIds rimNodeIds = segments[1].getRimNodeIdsSlice(0) if rimNodeIds: self._rimNodeIds[n2] = rimNodeIds continue - # create rim nodes + # create core box nodes + if self._boxCoordinates: + self._boxNodeIds[n2] = [] if self._boxNodeIds[n2] is None else self._boxNodeIds[n2] + nodesCountAcrossMajor = self.getCoreNodesCountAcrossMajor() + nodesCountAcrossMinor = self.getCoreNodesCountAcrossMinor() + for n3 in range(nodesCountAcrossMajor): + self._boxNodeIds[n2].append([]) + rx = self._boxCoordinates[0][n2][n3] + rd1 = self._boxCoordinates[1][n2][n3] + rd2 = self._boxCoordinates[2][n2][n3] + rd3 = self._boxCoordinates[3][n2][n3] + for n1 in range(nodesCountAcrossMinor): + nodeIdentifier = generateData.nextNodeIdentifier() + node = nodes.createNode(nodeIdentifier, nodetemplate) + fieldcache.setNode(node) + for nodeValue, rValue in zip([Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3], + [rx[n1], rd1[n1], rd2[n1], rd3[n1]]): + coordinates.setNodeParameters(fieldcache, -1, nodeValue, 1, rValue) + self._boxNodeIds[n2][n3].append(nodeIdentifier) + + # create rim nodes and transition nodes (if there are more than 1 layer of transition) nodesCountRim = len(self._rimCoordinates[0][0]) - self._rimNodeIds[n2] = [] - for n3 in range(nodesCountRim): - rx = self._rimCoordinates[0][n2][n3] - rd1 = self._rimCoordinates[1][n2][n3] - rd2 = self._rimCoordinates[2][n2][n3] - rd3 = None if isLinearThroughWall else self._rimCoordinates[3][n2][n3] + self._rimNodeIds[n2] = [] if self._rimNodeIds[n2] is None else self._rimNodeIds[n2] + nloop = nodesCountRim + (elementsCountTransition - 1) if self._isCore else nodesCountRim + for n3 in range(nloop): + n3p = n3 - (elementsCountTransition - 1) if self._isCore else n3 + if self._isCore and elementsCountTransition > 1 and n3 < (elementsCountTransition - 1): + # transition coordinates + rx = self._transitionCoordinates[0][n2][n3] + rd1 = self._transitionCoordinates[1][n2][n3] + rd2 = self._transitionCoordinates[2][n2][n3] + rd3 = self._transitionCoordinates[3][n2][n3] + else: + # rim coordinates + rx = self._rimCoordinates[0][n2][n3p] + rd1 = self._rimCoordinates[1][n2][n3p] + rd2 = self._rimCoordinates[2][n2][n3p] + rd3 = None if isLinearThroughWall else self._rimCoordinates[3][n2][n3p] ringNodeIds = [] for n1 in range(self._elementsCountAround): nodeIdentifier = generateData.nextNodeIdentifier() @@ -346,6 +1280,11 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): ringNodeIds.append(nodeIdentifier) self._rimNodeIds[n2].append(ringNodeIds) + # create a new list containing box node ids are located at the boundary + if self._isCore: + self._boxBoundaryNodeIds, self._boxBoundaryNodeToBoxId = ( + self._createBoxBoundaryNodeIdsList(startSkipCount, endSkipCount)) + if n2Only is not None: return @@ -357,9 +1296,68 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): for e2 in range(startSkipCount, elementsCountAlong - endSkipCount): self._rimElementIds[e2] = [] e2p = e2 + 1 + if self._isCore: + # create box elements + elementsCountAcrossMinor = self.getCoreNodesCountAcrossMinor() - 1 + elementsCountAcrossMajor = self.getCoreNodesCountAcrossMajor() - 1 + for e3 in range(elementsCountAcrossMajor): + e3p = e3 + 1 + for e1 in range(elementsCountAcrossMinor): + nids = [] + for n1 in [e1, e1 + 1]: + nids += [self._boxNodeIds[e2][e3][n1], self._boxNodeIds[e2][e3p][n1], + self._boxNodeIds[e2p][e3][n1], self._boxNodeIds[e2p][e3p][n1]] + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplateStd) + element.setNodesByIdentifier(eftStd, nids) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + + # create core transition elements + triplePointIndexesList = self.getTriplePointIndexes() + eftList = [None] * self._elementsCountAround + scalefactorsList = [None] * self._elementsCountAround + ringElementIds = [] + for e1 in range(self._elementsCountAround): + nids, nodeParameters, nodeLayouts = [], [], [] + n1p = (e1 + 1) % self._elementsCountAround + location = self.getTriplePointLocation(e1) + nodeLayoutTransition = generateData.getNodeLayoutTransition() + nodeLayoutTransitionTriplePoint = generateData.getNodeLayoutTransitionTriplePoint(location) + for n2 in [e2, e2 + 1]: + for n1 in [e1, n1p]: + nids += [self._boxBoundaryNodeIds[n2][n1]] + n3c, n1c = self._boxBoundaryNodeToBoxId[n2][n1] + nodeParameters.append(self.getBoxCoordinates(n1c, n2, n3c)) + nodeLayouts.append(nodeLayoutTransitionTriplePoint if n1 in triplePointIndexesList else + nodeLayoutTransition) + for n2 in [e2, e2 + 1]: + for n1 in [e1, n1p]: + nids += [self._rimNodeIds[n2][0][n1]] + nodeParameters.append(self.getRimCoordinates(n1, n2, 0)) + nodeLayouts.append(None) + eft = eftList[e1] + scalefactors = scalefactorsList[e1] + if not eft: + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + eftList[e1] = eft + scalefactorsList[e1] = scalefactors + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + ringElementIds.append(elementIdentifier) + self._rimElementIds[e2].append(ringElementIds) # create rim elements - for e3 in range(elementsCountRim): + nloop = elementsCountRim + (self._elementsCountTransition - 1) if self._isCore else elementsCountRim + for e3 in range(nloop): ringElementIds = [] for e1 in range(self._elementsCountAround): e1p = (e1 + 1) % self._elementsCountAround @@ -398,6 +1396,18 @@ def __init__(self, inSegments: list, outSegments: list): self._rimCoordinates = None # if set, (rx[], rd1[], rd2[], rd3[]) each over [n3][rim index] self._rimNodeIds = None # if set, nodeIdentifier[n3][rim index] + # parameters used for solid core + self._isCore = self._segments[0].getIsCore() + self._boxCoordinates = None # [nAlong][nAcrossMajor][nAcrossMinor] + self._transitionCoordinates = None + self._boxNodeIds = None + self._biSequence = None # sequence for bifurcation - either [1, 2, 3] or [1, 3, 2] + self._triSequence = None # sequence for trifurcation + self._boxIndexToSegmentNodeList = [] + # list[box index] giving list[(segment number, node index across major axis, node index across minor axis)] + self._segmentNodeToBoxIndex = [] + # list[segment number][node index across major axis][node index across minor axis] to boxIndex + def _calculateTrimSurfaces(self): """ Calculate surfaces for trimming adjacent segments so they can smoothly transition at junction. @@ -434,7 +1444,7 @@ def _calculateTrimSurfaces(self): weights = [] sumWeights = 0.0 maxWeight = 0.0 - phaseAngle = None + phaseAngle = None for os in range(self._segmentsCount): if os == s: continue @@ -514,7 +1524,7 @@ def _calculateTrimSurfaces(self): if lowestMaxProportionFromEnd == 0.0: xCentre = pathParameters[0][endIndex] else: - proportion = (1.0 - lowestMaxProportionFromEnd) if self._segmentsIn[s]\ + proportion = (1.0 - lowestMaxProportionFromEnd) if self._segmentsIn[s] \ else lowestMaxProportionFromEnd e = int(proportion) curveLocation = (e, proportion - e) @@ -709,7 +1719,6 @@ def _sampleMidPoint(self, segmentsParameterLists): elif segmentsCount == 4: si = sequence[1:3] - # print("sequence", sequence, "si", si) else: print("TubeNetworkMeshJunction._sampleMidPoint not fully implemented for segmentsCount =", segmentsCount) si = (1, 2) @@ -739,62 +1748,33 @@ def _sampleMidPoint(self, segmentsParameterLists): cd1, cd2 = cd2, cd1 return cx, cd1, cd2, cd3 - def sample(self, targetElementLength): + def _determineJunctionSequence(self): """ - Blend sampled d2 derivatives across 2-segment junctions with the same version. - Sample junction coordinates between second-from-end segment coordinates. - :param targetElementLength: Ignored here as always 2 elements across junction. + Determines the sequence of a junction. Currently only works for bifurcation and trifurcation. + There are two possible sequences for bifurcation - [0,1,2] and [0,2,1]. Sequence [0,1,2] indicates the second + segment is located on the right-hand side of the first segment, and [0,2,1] indicates the second segment is + located on the left-hand side of the first segment. + 2 1 1 2 + \ / \ / + 0 [0,1,2] 0 [0,2,1] + For trifurcation, only + shape is currently supported, not tetrahedral. The function determines plus sequence + [0,1,2,3], [0,1,3,2] or [0,2,3,1]. """ - if self._segmentsCount == 1: - return - if self._segmentsCount == 2: - TubeNetworkMeshSegment.blendSampledCoordinates( - self._segments[0], -1 if self._segmentsIn[0] else 0, - self._segments[1], -1 if self._segmentsIn[1] else 0) - return - - aroundCounts = [segment.getElementsCountAround() for segment in self._segments] - rimIndexesCount = 0 - - if self._segmentsCount == 3: - # numbers of elements directly connecting pairs of segments - # sequence = [0, 1, 2] - connectionCounts = [(aroundCounts[s] + aroundCounts[s - 2] - aroundCounts[s - 1]) // 2 for s in range(3)] - for s in range(3): - if (connectionCounts[s] < 1) or (aroundCounts[s] != (connectionCounts[s - 1] + connectionCounts[s])): - print("Can't make tube bifurcation between elements counts around", aroundCounts) - return - self._rimIndexToSegmentNodeList = [] # list[rim index] giving list[(segment index, node index around)] - self._segmentNodeToRimIndex = [[None] * aroundCounts[s] for s in range(self._segmentsCount)] - for s1 in range(3): - s2 = (s1 + 1) % 3 - startNodeIndex1 = (aroundCounts[s1] - connectionCounts[s1]) // 2 - startNodeIndex2 = connectionCounts[s1] // -2 - for n in range(connectionCounts[s1] + 1): - nodeIndex1 = startNodeIndex1 + (n if self._segmentsIn[s1] else (connectionCounts[s1] - n)) - if self._segmentNodeToRimIndex[s1][nodeIndex1] is None: - nodeIndex2 = startNodeIndex2 + ((connectionCounts[s1] - n) if self._segmentsIn[s2] else n) - if self._segmentNodeToRimIndex[s2][nodeIndex2] is None: - rimIndex = rimIndexesCount - # keep list in order from lowest s - self._rimIndexToSegmentNodeList.append( - [[s1, nodeIndex1], [s2, nodeIndex2]] if (s1 < s2) else - [[s2, nodeIndex2], [s1, nodeIndex1]]) - self._segmentNodeToRimIndex[s2][nodeIndex2] = rimIndex - rimIndexesCount += 1 - else: - rimIndex = self._segmentNodeToRimIndex[s2][nodeIndex2] - # keep list in order from lowest s - segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] - index = 0 - for i in range(len(segmentNodeList)): - if s1 < segmentNodeList[i][0]: - break - index += 1 - segmentNodeList.insert(index, [s1, nodeIndex1]) - self._segmentNodeToRimIndex[s1][nodeIndex1] = rimIndex + assert self._segmentsCount == 3 or self._segmentsCount == 4 + if self._segmentsCount == 3 and self._isCore: + # Find sequence order for bifurcation - only applies when core option is enabled: + # counter-clockwise sequence is [1, 2, 3] and clockwise sequence is [1, 3, 2] + pCoord = [] + for s in range(self._segmentsCount): + segment = self._segments[s] + n3 = (-1 if self._segmentsIn[s] else 0) if s >= 1 else 0 + p = segment.getBoxCoordinates(0, -2 if self._segmentsIn[s] else 1, n3)[0] + pCoord.append(p) + dist12 = magnitude(sub(pCoord[1], pCoord[0])) + dist13 = magnitude(sub(pCoord[2], pCoord[0])) + self._biSequence = [0, 1, 2] if dist12 > dist13 else [0, 2, 1] elif self._segmentsCount == 4: # only support plus + shape for now, not tetrahedral # determine plus + sequence [0, 1, 2, 3], [0, 1, 3, 2], [0, 2, 1, 3] or [0, 2, 3, 1] @@ -802,7 +1782,7 @@ def sample(self, targetElementLength): for s in range(self._segmentsCount): d1 = self._segments[s].getPathParameters()[1][-1 if self._segmentsIn[s] else 0] outDirections.append(normalize([-d for d in d1] if self._segmentsIn[s] else d1)) - ms = None + ns = None for s in range(1, self._segmentsCount): if dot(outDirections[0], outDirections[s]) > -0.9: ns = cross(outDirections[0], outDirections[s]) @@ -826,121 +1806,324 @@ def sample(self, targetElementLength): nexts = t angle = nextAngle sequence.append(nexts) - pairCount02 = aroundCounts[sequence[0]] + aroundCounts[sequence[2]] - pairCount13 = aroundCounts[sequence[1]] + aroundCounts[sequence[3]] + self._triSequence = sequence + else: + return + + def _sampleBifurcation(self, aroundCounts, acrossMajorCounts): + """ + Blackbox function for sampling bifurcation coordinates. The rim coordinates are first sampled, then the box + coordinates are sampled, if Core option is enabled. + :param aroundCounts: Number of elements around the tube. + :param acrossMajorCounts: Number of elements across major axis of the solid core. + :return rimIndexesCount, boxIndexesCount: Total number of rimIndexes and boxIndexes, respectively. + """ + assert self._segmentsCount == 3 + + # sample rim junction + rimIndexesCount = 0 + # numbers of elements directly connecting pairs of segments + # sequence = [0, 1, 2] + connectionCounts = [(aroundCounts[s] + aroundCounts[s - 2] - aroundCounts[s - 1]) // 2 for s in range(3)] + for s in range(3): + if (connectionCounts[s] < 1) or (aroundCounts[s] != (connectionCounts[s - 1] + connectionCounts[s])): + print("Can't make tube bifurcation between elements counts around", aroundCounts) + return + + self._rimIndexToSegmentNodeList = [] # list[rim index] giving list[(segment index, node index around)] + self._segmentNodeToRimIndex = [[None] * aroundCounts[s] for s in range(self._segmentsCount)] + for s1 in range(3): + s2 = (s1 + 1) % 3 + startNodeIndex1 = (aroundCounts[s1] - connectionCounts[s1]) // 2 + startNodeIndex2 = connectionCounts[s1] // -2 + for n in range(connectionCounts[s1] + 1): + nodeIndex1 = startNodeIndex1 + (n if self._segmentsIn[s1] else (connectionCounts[s1] - n)) + if self._segmentNodeToRimIndex[s1][nodeIndex1] is None: + nodeIndex2 = startNodeIndex2 + ((connectionCounts[s1] - n) if self._segmentsIn[s2] else n) + if self._segmentNodeToRimIndex[s2][nodeIndex2] is None: + rimIndex = rimIndexesCount + # keep list in order from lowest s + self._rimIndexToSegmentNodeList.append( + [[s1, nodeIndex1], [s2, nodeIndex2]] if (s1 < s2) else + [[s2, nodeIndex2], [s1, nodeIndex1]]) + self._segmentNodeToRimIndex[s2][nodeIndex2] = rimIndex + rimIndexesCount += 1 + else: + rimIndex = self._segmentNodeToRimIndex[s2][nodeIndex2] + # keep list in order from lowest s + segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] + index = 0 + for i in range(len(segmentNodeList)): + if s1 < segmentNodeList[i][0]: + break + index += 1 + segmentNodeList.insert(index, [s1, nodeIndex1]) + self._segmentNodeToRimIndex[s1][nodeIndex1] = rimIndex + + # sample box junction + boxIndexesCount = 0 + if self._isCore: + sequence = self._biSequence + connectionCounts = [(acrossMajorCounts[s] + acrossMajorCounts[s - 2] - acrossMajorCounts[s - 1]) // 2 for s + in range(3)] + midIndexes = [connectionCounts[s] - 1 for s in range(3)] + for s in range(3): + if acrossMajorCounts[s] != (connectionCounts[s - 1] + connectionCounts[s]): + print("Can't make core bifurcation between elements counts across major axis", acrossMajorCounts) + return + + nodesCountAcrossMinor = self._segments[0].getCoreNodesCountAcrossMinor() + nodesCountAcrossMajorList = [self._segments[s].getCoreNodesCountAcrossMajor() for s in range(3)] + self._segmentNodeToBoxIndex = \ + [[[None for _ in range(nodesCountAcrossMinor)] for _ in range(nodesCountAcrossMajorList[s])] + for s in range(self._segmentsCount)] + for s1 in range(3): + s2 = (s1 + 1) % 3 if sequence == [0, 1, 2] else (s1 - 1) % 3 + midIndex1 = midIndexes[s1] + midIndex2 = midIndexes[s2] + for m in range(connectionCounts[s1]): + for n in range(nodesCountAcrossMinor): + m1 = (m + midIndex1) if self._segmentsIn[s1] else (midIndex1 - m) % connectionCounts[s1] + m2 = (midIndex2 - m) % connectionCounts[s2] if self._segmentsIn[s2] else (m + midIndex2) + if m1 == midIndex1: + indexGroup = [[0, midIndexes[0], n], [1, midIndexes[1], n], [2, midIndexes[2], n]] + else: + indexGroup = [[s1, m1, n], [s2, m2, n]] if s1 < s2 else [[s2, m2, n], [s1, m1, n]] + if indexGroup not in self._boxIndexToSegmentNodeList: + self._boxIndexToSegmentNodeList.append(indexGroup) + boxIndexesCount += 1 + + for boxIndex in range(boxIndexesCount): + segmentNodeList = self._boxIndexToSegmentNodeList[boxIndex] + for segmentNode in segmentNodeList: + s, m, n = segmentNode + self._segmentNodeToBoxIndex[s][m][n] = boxIndex + + return rimIndexesCount, boxIndexesCount + + def _sampleTrifurcation(self, aroundCounts, acrossMajorCounts): + """ + Blackbox function for sampling trifurcation coordinates. The rim coordinates are first sampled, then the box + coordinates are sampled, if Core option is enabled. + :param aroundCounts: Number of elements around the tube. + :param acrossMajorCounts: Number of elements across major axis of the solid core. + :return rimIndexesCount, boxIndexesCount: Total number of rimIndexes and boxIndexes, respectively. + """ + assert self._segmentsCount == 4 + + # sample rim junction + rimIndexesCount = 0 + sequence = self._triSequence + pairCount02 = aroundCounts[sequence[0]] + aroundCounts[sequence[2]] + pairCount13 = aroundCounts[sequence[1]] + aroundCounts[sequence[3]] + throughCount02 = ((pairCount02 - pairCount13) // 2) if (pairCount02 > pairCount13) else 0 + throughCount13 = ((pairCount13 - pairCount02) // 2) if (pairCount13 > pairCount02) else 0 + throughCounts = [throughCount02, throughCount13, throughCount02, throughCount13] + # numbers of elements directly connecting pairs of segments + freeAroundCounts = [aroundCounts[sequence[s]] - throughCounts[s] for s in range(self._segmentsCount)] + if freeAroundCounts[0] == freeAroundCounts[2]: + count03 = freeAroundCounts[3] // 2 + count12 = freeAroundCounts[1] // 2 + connectionCounts = [count03, count12, count12, count03] + elif freeAroundCounts[1] == freeAroundCounts[3]: + count03 = freeAroundCounts[0] // 2 + count12 = freeAroundCounts[2] // 2 + connectionCounts = [count03, count12, count12, count03] + else: + connectionCounts = [((freeAroundCounts[s] + freeAroundCounts[(s + 1) % self._segmentsCount] + - freeAroundCounts[s - 1] + (s % 2)) // 2) for s in range(self._segmentsCount)] + + for s in range(self._segmentsCount): + if (aroundCounts[sequence[s]] != (connectionCounts[s - 1] + throughCounts[s] + connectionCounts[s])): + print("Can't make tube junction between elements counts around", aroundCounts) + return + + self._rimIndexToSegmentNodeList = [] # list[rim index] giving list[(segment index, node index around)] + self._segmentNodeToRimIndex = [[None] * aroundCounts[s] for s in range(self._segmentsCount)] + for os1 in range(self._segmentsCount): + os2 = (os1 + 1) % self._segmentsCount + s1 = sequence[os1] + s2 = sequence[os2] + s3 = sequence[(os1 + 2) % self._segmentsCount] + halfThroughCount = throughCounts[os1] // 2 + os1ConnectionCount = connectionCounts[os1] + os2ConnectionCount = connectionCounts[os2] + if self._segmentsIn[s1]: + startNodeIndex1 = (aroundCounts[s1] - os1ConnectionCount) // 2 + else: + startNodeIndex1 = os1ConnectionCount // -2 + if self._segmentsIn[s2]: + startNodeIndex2 = os1ConnectionCount // -2 + else: + startNodeIndex2 = (aroundCounts[s2] - os1ConnectionCount) // 2 + if self._segmentsIn[s3]: + startNodeIndex3h = os2ConnectionCount // -2 + startNodeIndex3l = startNodeIndex3h - (os2ConnectionCount - os1ConnectionCount) + else: + startNodeIndex3l = (aroundCounts[s3] - os2ConnectionCount) // 2 + startNodeIndex3h = startNodeIndex3l + (os2ConnectionCount - os1ConnectionCount) + + for n in range(-halfThroughCount, os1ConnectionCount + 1 + halfThroughCount): + n1 = startNodeIndex1 + (n if self._segmentsIn[s1] else (os1ConnectionCount - n)) + segmentIndexes = [s1] + nodeIndexes = [n1 % aroundCounts[s1]] + if 0 <= n <= os1ConnectionCount: + n2 = startNodeIndex2 + ((os1ConnectionCount - n) if self._segmentsIn[s2] else n) + segmentIndexes.append(s2) + nodeIndexes.append(n2 % aroundCounts[s2]) + if halfThroughCount and ((n <= 0) or (n >= os1ConnectionCount)): + n3 = ((startNodeIndex3l if n <= 0 else startNodeIndex3h) + + ((os2ConnectionCount - n) if self._segmentsIn[s3] else n)) + segmentIndexes.append(s3) + nodeIndexes.append(n3 % aroundCounts[s3]) + + rimIndex = None + for i in range(len(segmentIndexes)): + ri = self._segmentNodeToRimIndex[segmentIndexes[i]][nodeIndexes[i]] + if ri is not None: + rimIndex = ri + break + if rimIndex is None: + # new rim index + rimIndex = rimIndexesCount + rimIndexesCount += 1 + segmentNodeList = [] + self._rimIndexToSegmentNodeList.append(segmentNodeList) + else: + segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] + # build maps: rim index <--> segment index, node index + for i in range(len(segmentIndexes)): + segmentIndex = segmentIndexes[i] + nodeIndex = nodeIndexes[i] + if self._segmentNodeToRimIndex[segmentIndex][nodeIndex] is None: + # keep segment node list in order from lowest segment index + index = 0 + for j in range(len(segmentNodeList)): + if segmentIndex < segmentNodeList[j][0]: + break + index += 1 + segmentNodeList.insert(index, [segmentIndex, nodeIndex]) + self._segmentNodeToRimIndex[segmentIndex][nodeIndex] = rimIndex + + # sample box junction + boxIndexesCount = 0 + if self._isCore: + pCoord = [] + for s in range(self._segmentsCount): + segment = self._segments[s] + n3 = (-1 if self._segmentsIn[s] else 0) if s >= 1 else 0 + p = segment.getBoxCoordinates(0, -2 if self._segmentsIn[s] else 1, n3)[0] + pCoord.append(p) + dist12 = magnitude(sub(pCoord[sequence[1]], pCoord[0])) + dist13 = magnitude(sub(pCoord[sequence[-1]], pCoord[0])) + if sequence == [0, 1, 3, 2]: + sequence = [0, 2, 3, 1] if dist12 < dist13 else sequence + elif sequence == [0, 1, 2, 3]: + sequence = [0, 2, 1, 3] if dist12 < dist13 else sequence + + pairCount02 = acrossMajorCounts[sequence[0]] + acrossMajorCounts[sequence[2]] + pairCount13 = acrossMajorCounts[sequence[1]] + acrossMajorCounts[sequence[3]] throughCount02 = ((pairCount02 - pairCount13) // 2) if (pairCount02 > pairCount13) else 0 throughCount13 = ((pairCount13 - pairCount02) // 2) if (pairCount13 > pairCount02) else 0 throughCounts = [throughCount02, throughCount13, throughCount02, throughCount13] - # numbers of elements directly connecting pairs of segments - # print("sample sequence", sequence) - # for s in range(self._segmentsCount): - # print(s, ":", sequence[s], "=", aroundCounts[sequence[s]], sequence[(s - 1) % self._segmentsCount], '=', - # aroundCounts[sequence[(s - 1) % self._segmentsCount]], sequence[s - 1], "=", aroundCounts[sequence[s - 1]], throughCounts[s], (s % 2)) - freeAroundCounts = [aroundCounts[sequence[s]] - throughCounts[s] for s in range(self._segmentsCount)] - if freeAroundCounts[0] == freeAroundCounts[2]: - count03 = freeAroundCounts[3] // 2 - count12 = freeAroundCounts[1] // 2 + freeAcrossCounts = [acrossMajorCounts[sequence[s]] - throughCounts[s] for s in range(self._segmentsCount)] + if freeAcrossCounts[0] == freeAcrossCounts[2]: + count03 = freeAcrossCounts[3] // 2 + count12 = freeAcrossCounts[1] // 2 connectionCounts = [count03, count12, count12, count03] - elif freeAroundCounts[1] == freeAroundCounts[3]: - count03 = freeAroundCounts[0] // 2 - count12 = freeAroundCounts[2] // 2 + elif freeAcrossCounts[1] == freeAcrossCounts[3]: + count03 = freeAcrossCounts[0] // 2 + count12 = freeAcrossCounts[2] // 2 connectionCounts = [count03, count12, count12, count03] else: - connectionCounts = [((freeAroundCounts[s] + freeAroundCounts[(s + 1) % self._segmentsCount] - - freeAroundCounts[s - 1] + (s % 2)) // 2) for s in range(self._segmentsCount)] - # print("aroundCounts", aroundCounts) - # print("sequence", sequence) - # print("throughCounts", throughCounts) - # print("freeAroundCounts", freeAroundCounts) - # print("connectionCounts", connectionCounts) + connectionCounts = [((freeAcrossCounts[s] + freeAcrossCounts[(s + 1) % self._segmentsCount] + - freeAcrossCounts[s - 1] + (s % 2)) // 2) for s in range(self._segmentsCount)] for s in range(self._segmentsCount): - # print("s", s, "around", aroundCounts[sequence[s]], "connL", connectionCounts[s - 1], "through", throughCounts[s], "connR", connectionCounts[s]) - if (aroundCounts[sequence[s]] != (connectionCounts[s - 1] + throughCounts[s] + connectionCounts[s])): - print("Can't make tube junction between elements counts around", aroundCounts) + if (acrossMajorCounts[sequence[s]] != ( + connectionCounts[s - 1] + throughCounts[s] + connectionCounts[s])): + print("Can't make tube junction between elements counts around", acrossMajorCounts) return - self._rimIndexToSegmentNodeList = [] # list[rim index] giving list[(segment index, node index around)] - self._segmentNodeToRimIndex = [[None] * aroundCounts[s] for s in range(self._segmentsCount)] + nodesCountAcrossMinor = self._segments[0].getCoreNodesCountAcrossMinor() + nodesCountAcrossMajorList = [self._segments[s].getCoreNodesCountAcrossMajor() for s in range(4)] + midIndexes = [nodesCountAcrossMajor // 2 for nodesCountAcrossMajor in nodesCountAcrossMajorList] + halfThroughCounts = [throughCounts[s] // 2 for s in range(4)] + + connectingIndexesList = [[] for s in range(4)] + for s in range(4): + s1 = (s - 1) % self._segmentsCount + s2 = (s + 1) % self._segmentsCount + shift = (midIndexes[sequence[s2]] - midIndexes[sequence[s1]]) // 2 + midIndex = midIndexes[sequence[s]] + halfThroughCount = halfThroughCounts[s] + connectingIndexesList[sequence[s]] = \ + ([midIndex - halfThroughCount, midIndex + halfThroughCount] if halfThroughCount \ + else [midIndex + shift]) + + self._boxIndexToSegmentNodeList = [] + self._segmentNodeToBoxIndex = \ + [[[None for _ in range(nodesCountAcrossMinor)] for _ in range(nodesCountAcrossMajorList[s])] + for s in range(self._segmentsCount)] + for os1 in range(self._segmentsCount): os2 = (os1 + 1) % self._segmentsCount s1 = sequence[os1] s2 = sequence[os2] s3 = sequence[(os1 + 2) % self._segmentsCount] - halfThroughCount = throughCounts[os1] // 2 - os1ConnectionCount = connectionCounts[os1] - os2ConnectionCount = connectionCounts[os2] - if self._segmentsIn[s1]: - startNodeIndex1 = (aroundCounts[s1] - os1ConnectionCount) // 2 - else: - startNodeIndex1 = os1ConnectionCount // -2 - if self._segmentsIn[s2]: - startNodeIndex2 = os1ConnectionCount // -2 - else: - startNodeIndex2 = (aroundCounts[s2] - os1ConnectionCount) // 2 - if self._segmentsIn[s3]: - startNodeIndex3h = os2ConnectionCount // -2 - startNodeIndex3l = startNodeIndex3h - (os2ConnectionCount - os1ConnectionCount) - else: - startNodeIndex3l = (aroundCounts[s3] - os2ConnectionCount) // 2 - startNodeIndex3h = startNodeIndex3l + (os2ConnectionCount - os1ConnectionCount) - - # print(os1, "half", halfThroughCount, "start", startNodeIndex1, - # "in1", self._segmentsIn[s1], "in2", self._segmentsIn[s2], "in3", self._segmentsIn[s3]) - for n in range(-halfThroughCount, os1ConnectionCount + 1 + halfThroughCount): - n1 = startNodeIndex1 + (n if self._segmentsIn[s1] else (os1ConnectionCount - n)) - segmentIndexes = [s1] - nodeIndexes = [n1 % aroundCounts[s1]] - if 0 <= n <= os1ConnectionCount: - n2 = startNodeIndex2 + ((os1ConnectionCount - n) if self._segmentsIn[s2] else n) - segmentIndexes.append(s2) - nodeIndexes.append(n2 % aroundCounts[s2]) - if halfThroughCount and ((n <= 0) or (n >= os1ConnectionCount)): - n3 = ((startNodeIndex3l if n <= 0 else startNodeIndex3h) + - ((os2ConnectionCount - n) if self._segmentsIn[s3] else n)) - segmentIndexes.append(s3) - nodeIndexes.append(n3 % aroundCounts[s3]) - # print("seg node", [[segmentIndexes[i], nodeIndexes[i]] for i in range(len(segmentIndexes))]) - rimIndex = None - for i in range(len(segmentIndexes)): - ri = self._segmentNodeToRimIndex[segmentIndexes[i]][nodeIndexes[i]] - if ri is not None: - rimIndex = ri - # for j in range(len(segmentIndexes)): - # sn = [segmentIndexes[j], nodeIndexes[j]] - # if sn not in self._rimIndexToSegmentNodeList[rimIndex]: - # print(" add", sn, "to rim index", ri, "segment nodes", self._rimIndexToSegmentNodeList[rimIndex]) - break - if rimIndex is None: - # new rim index - rimIndex = rimIndexesCount - rimIndexesCount += 1 - segmentNodeList = [] - self._rimIndexToSegmentNodeList.append(segmentNodeList) - else: - segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] - # build maps: rim index <--> segment index, node index - for i in range(len(segmentIndexes)): - segmentIndex = segmentIndexes[i] - nodeIndex = nodeIndexes[i] - if self._segmentNodeToRimIndex[segmentIndex][nodeIndex] is None: - # keep segment node list in order from lowest segment index - index = 0 - for j in range(len(segmentNodeList)): - if segmentIndex < segmentNodeList[j][0]: - break - index += 1 - segmentNodeList.insert(index, [segmentIndex, nodeIndex]) - self._segmentNodeToRimIndex[segmentIndex][nodeIndex] = rimIndex - else: - print("Tube network mesh not implemented for", self._segmentsCount, "segments at junction") - return - - if not rimIndexesCount: - return - - # for rimIndex in range(rimIndexesCount): - # print("rim index", rimIndex, self._rimIndexToSegmentNodeList[rimIndex]) - + s4 = sequence[(os1 + 3) % self._segmentsCount] + aStartNodeIndex = midIndexes[s1] + nodesCountAcrossMajor1 = nodesCountAcrossMajorList[s1] + nodesCountAcrossMajor2 = nodesCountAcrossMajorList[s2] + connectingIndexes1 = connectingIndexesList[s1] + connectingIndexes2 = connectingIndexesList[s2] + connectingIndexes3 = connectingIndexesList[s3] + throughIndexes = list(range(connectingIndexes1[0] + 1, connectingIndexes1[1])) \ + if throughCounts[os1] else [None] + + for m in range(acrossMajorCounts[s1] // 2): + m1 = (m + aStartNodeIndex) if self._segmentsIn[s1] else (aStartNodeIndex - m) % connectionCounts[s1] + for n in range(nodesCountAcrossMinor): + if m1 in throughIndexes: + m3 = midIndexes[s3] + indexGroup = [[s1, m1, n], [s3, m3, n]] if s1 < s3 else [[s3, m3, n], [s1, m1, n]] + elif m1 in connectingIndexes1: + if all(v == 0 for v in throughCounts): + indexGroup = [[s1, m1, n], [s2, m1, n], [s3, m1, n], [s4, m1, n]] + else: + if throughCounts[os1]: + i = 1 if self._segmentsIn[s1] else 0 + m2 = connectingIndexes2[0] + m3 = connectingIndexes3[i] + indexGroup = [[s1, m1, n], [s2, m2, n], [s3, m3, n]] + else: + if self._segmentsIn[s1] == self._segmentsIn[s2]: + if throughCounts[os1]: + m2 = nodesCountAcrossMajor2 - m1 - 1 + else: + m2 = nodesCountAcrossMajor2 - (m1 + 1) if not self._segmentsIn[s1] else ( + nodesCountAcrossMajor1 - (m1 + 1)) + else: + m2 = m1 - throughCounts[s1] if self._segmentsIn[s1] else m1 + indexGroup = [[s1, m1, n], [s2, m2, n]] if s1 < s2 else [[s2, m2, n], [s1, m1, n]] + indexGroup = sorted(indexGroup, key=lambda x: (x[0], x[0]), reverse=False) + if indexGroup not in self._boxIndexToSegmentNodeList: + self._boxIndexToSegmentNodeList.append(indexGroup) + boxIndexesCount += 1 + + for boxIndex in range(boxIndexesCount): + segmentNodeList = self._boxIndexToSegmentNodeList[boxIndex] + for segmentNode in segmentNodeList: + s, m, n = segmentNode + self._segmentNodeToBoxIndex[s][m][n] = boxIndex + + return rimIndexesCount, boxIndexesCount + + def _optimiseRimIndexes(self, aroundCounts, rimIndexesCount): + """ + Iterates through a number of permutations to find the most optimised lookup table for rim indexes. + :param aroundCounts: Number of elements around the tube. + :param rimIndexesCount: Total number of rim indexes. + """ # get node indexes giving the lowest sum of distances between adjoining points on outer sampled tubes permutationCount = 1 for count in aroundCounts: @@ -974,8 +2157,6 @@ def sample(self, targetElementLength): break indexes[s] = 0 - # print("minIndexes", minIndexes) - # offset node indexes by minIndexes for rimIndex in range(rimIndexesCount): segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] @@ -985,23 +2166,340 @@ def sample(self, targetElementLength): self._segmentNodeToRimIndex[s][nodeIndex] = rimIndex segmentNode[1] = nodeIndex + def sample(self, targetElementLength): + """ + Blend sampled d2 derivatives across 2-segment junctions with the same version. + Sample junction coordinates between second-from-end segment coordinates. + :param targetElementLength: Ignored here as always 2 elements across junction. + """ + if self._segmentsCount == 1: + return + if self._segmentsCount == 2: + TubeNetworkMeshSegment.blendSampledCoordinates( + self._segments[0], -1 if self._segmentsIn[0] else 0, + self._segments[1], -1 if self._segmentsIn[1] else 0) + return + + aroundCounts = [segment.getElementsCountAround() for segment in self._segments] + acrossMajorCounts = [segment.getElementsCountAcrossMajor() for segment in self._segments] + + # determine junction sequence + self._determineJunctionSequence() + + if self._segmentsCount == 3: + rimIndexesCount, boxIndexesCount = self._sampleBifurcation(aroundCounts, acrossMajorCounts) + + elif self._segmentsCount == 4: + rimIndexesCount, boxIndexesCount = self._sampleTrifurcation(aroundCounts, acrossMajorCounts) + + else: + print("Tube network mesh not implemented for", self._segmentsCount, "segments at junction") + return + + if not rimIndexesCount: + return + + # optimise rim indexes + self._optimiseRimIndexes(aroundCounts, rimIndexesCount) + # sample rim coordinates - nodesCountRim = self._segments[0].getNodesCountRim() + elementsCountTransition = self._segments[0].getElementsCountTransition() + nodesCountRim = self._segments[0].getNodesCountRim() + (elementsCountTransition - 1) if self._isCore else ( + self._segments[0].getNodesCountRim()) rx, rd1, rd2, rd3 = [ [[None] * rimIndexesCount for _ in range(nodesCountRim)] for i in range(4)] self._rimCoordinates = (rx, rd1, rd2, rd3) for n3 in range(nodesCountRim): + n3p = n3 - (elementsCountTransition - 1) if self._isCore else n3 for rimIndex in range(rimIndexesCount): segmentNodeList = self._rimIndexToSegmentNodeList[rimIndex] # segments have been ordered from lowest to highest s index segmentsParameterLists = [] for s, n1 in segmentNodeList: + if self._isCore and n3 < (elementsCountTransition - 1): + segmentsParameterLists.append( + self._segments[s].getTransitionCoordinatesListAlong( + n1, [-2, -1] if self._segmentsIn[s] else [1, 0], n3)) + else: + segmentsParameterLists.append( + self._segments[s].getRimCoordinatesListAlong( + n1, [-2, -1] if self._segmentsIn[s] else [1, 0], n3p)) + rx[n3][rimIndex], rd1[n3][rimIndex], rd2[n3][rimIndex], rd3[n3][rimIndex] = \ + self._sampleMidPoint(segmentsParameterLists) + + # sample box coordinates + if self._isCore: + bx, bd1, bd2, bd3 = [[None] * boxIndexesCount for _ in range(4)] + self._boxCoordinates = (bx, bd1, bd2, bd3) + for boxIndex in range(boxIndexesCount): + segmentNodeList = self._boxIndexToSegmentNodeList[boxIndex] + segmentsParameterLists = [] + for s, n3, n1 in segmentNodeList: segmentsParameterLists.append( - self._segments[s].getRimCoordinatesListAlong( + self._segments[s].getBoxCoordinatesListAlong( n1, [-2, -1] if self._segmentsIn[s] else [1, 0], n3)) - rx[n3][rimIndex], rd1[n3][rimIndex], rd2[n3][rimIndex], rd3[n3][rimIndex] = \ + bx[boxIndex], bd1[boxIndex], bd2[boxIndex], bd3[boxIndex] = \ self._sampleMidPoint(segmentsParameterLists) + def _createBoxBoundaryNodeIdsList(self, s): + """ + Creates a list (in a circular format similar to other rim node id lists) of core box node ids that are + located at the boundary of the core. This list is used to easily stitch inner rim nodes with box nodes. + Used specifically for solid core at the junction. + :param s: Index for identifying segments. + :return: A list of box node ids stored in a circular format, and a lookup list that translates indexes used in + boxBoundaryNodeIds list to indexes that can be used in boxCoordinates list. + """ + boxBoundaryNodeIds = [] + boxBoundaryNodeToBoxId = [] + elementsCountAcrossMajor = self._segments[s].getElementsCountAcrossMajor() + elementsCountAcrossMinor = self._segments[s].getElementsCountAcrossMinor() + elementsCountTransition = self._segments[s].getElementsCountAcrossTransition() + boxElementsCountRow = (elementsCountAcrossMajor - 2 * elementsCountTransition) + 1 + boxElementsCountColumn = (elementsCountAcrossMinor - 2 * elementsCountTransition) + 1 + + for n3 in range(boxElementsCountRow): + if n3 == 0 or n3 == boxElementsCountRow - 1: + ids = self._boxNodeIds[s][n3] if n3 == 0 else self._boxNodeIds[s][n3][::-1] + n1List = list(range(boxElementsCountColumn)) if n3 == 0 else ( + list(range(boxElementsCountColumn - 1, -1, -1))) + boxBoundaryNodeIds += [ids[c] for c in range(boxElementsCountColumn)] + for n1 in n1List: + boxBoundaryNodeToBoxId.append([n3, n1]) + else: + for n1 in [-1, 0]: + boxBoundaryNodeIds.append(self._boxNodeIds[s][n3][n1]) + boxBoundaryNodeToBoxId.append([n3, n1]) + + start = elementsCountAcrossMajor - 4 - 2 * (elementsCountTransition - 1) + idx = elementsCountAcrossMinor - 2 * (elementsCountTransition - 1) + for n in range(int(start), -1, -1): + boxBoundaryNodeIds.append(boxBoundaryNodeIds.pop(idx + 2 * n)) + boxBoundaryNodeToBoxId.append(boxBoundaryNodeToBoxId.pop(idx + 2 * n)) + + nloop = elementsCountAcrossMinor // 2 - elementsCountTransition + for _ in range(nloop): + boxBoundaryNodeIds.insert(len(boxBoundaryNodeIds), boxBoundaryNodeIds.pop(0)) + boxBoundaryNodeToBoxId.insert(len(boxBoundaryNodeToBoxId), + boxBoundaryNodeToBoxId.pop(0)) + + return boxBoundaryNodeIds, boxBoundaryNodeToBoxId + + def _getBoxCoordinates(self, n1): + return (self._boxCoordinates[0][n1], self._boxCoordinates[1][n1], + self._boxCoordinates[2][n1], self._boxCoordinates[3][n1]) + + def _getRimCoordinates(self, n1): + return (self._rimCoordinates[0][0][n1], self._rimCoordinates[1][0][n1], + self._rimCoordinates[2][0][n1], self._rimCoordinates[3][0][n1]) + + def _generateBoxElements(self, s, n2, mesh, elementtemplate, coordinates, segment, generateData): + """ + Blackbox function for generating core box elements at a junction. + """ + annotationMeshGroups = generateData.getAnnotationMeshGroups(segment.getAnnotationTerms()) + boxElementsCountAcrossMinor = self._segments[0].getCoreNodesCountAcrossMinor() - 1 + boxElementsCountAcrossMajor = [self._segments[s].getCoreNodesCountAcrossMajor() - 1 + for s in range(self._segmentsCount)] + acrossMajorCounts = [segment.getElementsCountAcrossMajor() for segment in self._segments] + is6WayTriplePoint = True if (((max(acrossMajorCounts) - 2) // 2) == (min(acrossMajorCounts) - 2) + and (self._segmentsCount == 3)) else False + + eftList = [[None] * boxElementsCountAcrossMinor for _ in range(boxElementsCountAcrossMajor[s])] + scalefactorsList = [[None] * boxElementsCountAcrossMinor for _ in range(boxElementsCountAcrossMajor[s])] + + nodeLayout6Way = generateData.getNodeLayout6Way() + nodeLayout8Way = generateData.getNodeLayout8Way() + nodeLayoutFlipD2 = generateData.getNodeLayoutFlipD2() + nodeLayoutBifurcation = generateData.getNodeLayoutBifurcation() + + for e3 in range(boxElementsCountAcrossMajor[s]): + for e1 in range(boxElementsCountAcrossMinor): + e3p = (e3 + 1) + nids, nodeParameters, nodeLayouts = [], [], [] + for n1 in [e1, e1 + 1]: + for n3 in [e3, e3p]: + nids.append(segment.getBoxNodeIds(n1, n2, n3)) + boxCoordinates = segment.getBoxCoordinates(n1, n2, n3) + nodeParameters.append(boxCoordinates) + nodeLayouts.append(None) + for n3 in [e3, e3p]: + boxIndex = self._segmentNodeToBoxIndex[s][n3][n1] + nids.append(self._boxNodeIds[s][n3][n1]) + nodeParameters.append(self._getBoxCoordinates(boxIndex)) + segmentNodesCount = len(self._boxIndexToSegmentNodeList[boxIndex]) + if is6WayTriplePoint and (segmentNodesCount == 3) and self._segmentsCount == 3: + nodeLayouts.append(nodeLayoutBifurcation) + elif self._segmentsIn[s] and (segmentNodesCount == 3) and self._segmentsCount == 4: + location = 1 if e3 < boxElementsCountAcrossMajor[s] // 2 else 2 + nodeLayoutTrifurcation = generateData.getNodeLayoutTrifurcation(location) + nodeLayouts.append(nodeLayout6Way if self._triSequence == [0, 1, 3, 2] else + nodeLayoutTrifurcation) + else: + nodeLayouts.append(nodeLayoutFlipD2 if (segmentNodesCount == 2) else + nodeLayout6Way if (segmentNodesCount == 3) else + nodeLayout8Way) + + if not self._segmentsIn[s]: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft = eftList[e3][e1] + scalefactors = scalefactorsList[e3][e1] + if not eft: + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + eftList[e3][e1] = eft + scalefactorsList[e3][e1] = scalefactors + elementtemplate.defineField(coordinates, -1, eft) + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + + def _generateTransitionElements(self, s, n2, mesh, elementtemplate, coordinates, segment, generateData, + elementsCountAround, boxBoundaryNodeIds, boxBoundaryNodeToBoxId): + """ + Blackbox function for generating core transition elements at a junction. + """ + annotationMeshGroups = generateData.getAnnotationMeshGroups(segment.getAnnotationTerms()) + nodesCountAcrossMinor = self._segments[0].getCoreNodesCountAcrossMinor() + nodesCountAcrossMajor = [self._segments[s].getCoreNodesCountAcrossMajor() for s in + range(self._segmentsCount)] + acrossMajorCounts = [segment.getElementsCountAcrossMajor() for segment in self._segments] + eftList = [None] * elementsCountAround + scalefactorsList = [None] * elementsCountAround + + triplePointIndexesList = segment.getTriplePointIndexes() + is6WayTriplePoint = True if (((max(acrossMajorCounts) - 2) // 2) == (min(acrossMajorCounts) - 2) + and (self._segmentsCount == 3)) else False + pSegment = acrossMajorCounts.index(max(acrossMajorCounts)) + topMidIndex = (nodesCountAcrossMajor[pSegment] // 2) + (nodesCountAcrossMinor // 2) + bottomMidIndex = elementsCountAround - topMidIndex + midIndexes = [topMidIndex, bottomMidIndex] + + nodeLayout6Way = generateData.getNodeLayout6Way() + nodeLayout8Way = generateData.getNodeLayout8Way() + nodeLayoutFlipD2 = generateData.getNodeLayoutFlipD2() + nodeLayoutTransition = generateData.getNodeLayoutTransition() + nodeLayoutBifurcationTransition = generateData.getNodeLayoutBifurcationTransition() + + for e1 in range(elementsCountAround): + nids, nodeParameters, nodeLayouts = [], [], [] + n1p = (e1 + 1) % elementsCountAround + oLocation = segment.getTriplePointLocation(e1) + for n3 in range(2): + if n3 == 0: # core box region + for n1 in [e1, n1p]: + nids.append(segment.getBoxBoundaryNodeIds(n1, n2)) + n3c, n1c = segment.getBoxBoundaryNodeToBoxId(n1, n2) + nodeParameters.append(segment.getBoxCoordinates(n1c, n2, n3c)) + nodeLayoutTransitionTriplePoint = ( + generateData.getNodeLayoutTransitionTriplePoint(oLocation)) + nodeLayouts.append(nodeLayoutTransitionTriplePoint if n1 in triplePointIndexesList + else nodeLayoutTransition) + for n1 in [e1, n1p]: + nids.append(boxBoundaryNodeIds[n1]) + n3c, n1c = boxBoundaryNodeToBoxId[n1] + boxIndex = self._segmentNodeToBoxIndex[s][n3c][n1c] + nodeParameters.append(self._getBoxCoordinates(boxIndex)) + segmentNodesCount = len(self._boxIndexToSegmentNodeList[boxIndex]) + if segmentNodesCount == 3: # 6-way node + if is6WayTriplePoint: # Special 6-way triple point case + if (n1 in triplePointIndexesList or (s == pSegment and n1 in midIndexes)): + # 6-way AND triple-point node - only applies to bifurcations + location = (midIndexes.index(n1) + 1) if oLocation == 0 else oLocation + if (s == 1 and n1 == n1p) or (s == 2 and n1 == e1): + location = 3 if abs(location) == 1 else location + elif (s == 1 and n1 != n1p) or (s == 2 and n1 != e1): + location = 4 if abs(location) == 2 else location + nodeLayout = generateData.getNodeLayout6WayTriplePoint(location) + else: + nodeLayout = nodeLayoutBifurcationTransition + elif self._segmentsCount == 4 and self._segmentsIn[s]: # Trifurcation case + location = \ + 1 if (e1 < elementsCountAround // 4) or (e1 >= 3 * elementsCountAround // 4) else 2 + nodeLayoutTrifurcation = generateData.getNodeLayoutTrifurcation(location) + nodeLayout = nodeLayout6Way if self._triSequence == [0, 1, 3, 2] else ( + nodeLayoutTrifurcation) + else: + nodeLayout = nodeLayout6Way + elif segmentNodesCount == 4: # 8-way node + nodeLayout = nodeLayout8Way + elif n1 in triplePointIndexesList: # triple-point node + location = oLocation + if self._segmentsCount == 3: # bifurcation + sequence = self._biSequence + condition1 = oLocation > 0 + condition2 = oLocation < 0 + if self._segmentsIn.count(True) == 0: + condition = condition1 if sequence == [0, 2, 1] else condition2 + location *= -1 if s == 2 or (s == 1 and condition) else 1 + elif self._segmentsIn.count(True) == 1: + condition = condition1 if sequence == [0, 2, 1] else condition2 + location *= -1 if s == 2 and condition else 1 + elif self._segmentsIn.count(True) == 2: + condition = condition2 if sequence == [0, 2, 1] else condition1 + location *= -1 if s == 1 and condition else 1 + elif self._segmentsIn.count(True) == 3: + condition = condition2 if sequence == [0, 2, 1] else condition1 + location *= -1 if (s == 1 and condition) or (s == 2) else 1 + elif self._segmentsCount == 4: # trifurcation + sequence = self._triSequence + s0 = (s - 1) % self._segmentsCount + s1 = (s + 1) % self._segmentsCount + if sequence == [0, 1, 3, 2]: + if self._segmentsIn == [True, False, False, False] and self._segmentsIn[s1]: + location = (oLocation) * -1 + elif self._segmentsIn == [True, True, False, False] and \ + self._segmentsIn[s] != self._segmentsIn[s1]: + location = abs(oLocation) + elif sequence == [0, 1, 2, 3] and \ + (self._segmentsIn[s1] or all(not self._segmentsIn[n] for n in [s, s0, s1])): + location = abs(oLocation) + nodeLayout = generateData.getNodeLayoutTransitionTriplePoint(location) + else: + nodeLayout = nodeLayoutTransition + nodeLayouts.append(nodeLayout) + + if not self._segmentsIn[s]: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + else: # rim region + for n1 in [e1, n1p]: + nids.append(segment.getRimNodeId(n1, n2, 0)) + nodeParameters.append(segment.getRimCoordinates(n1, n2, 0)) + nodeLayouts.append(None) + for n1 in [e1, n1p]: + rimIndex = self._segmentNodeToRimIndex[s][n1] + nids.append(self._rimNodeIds[0][rimIndex]) + nodeParameters.append(self._getRimCoordinates(rimIndex)) + segmentNodesCount = len(self._rimIndexToSegmentNodeList[rimIndex]) + nodeLayouts.append(nodeLayoutFlipD2 if (segmentNodesCount == 2) else + nodeLayout6Way if (segmentNodesCount == 3) else + nodeLayout8Way) + if not self._segmentsIn[s]: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft = eftList[e1] + scalefactors = scalefactorsList[e1] + if not eft: + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + eftList[e1] = eft + scalefactorsList[e1] = scalefactors + elementtemplate.defineField(coordinates, -1, eft) + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + def generateMesh(self, generateData: TubeNetworkMeshGenerateData): if generateData.isShowTrimSurfaces(): dimension = generateData.getMeshDimension() @@ -1020,11 +2518,19 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData): return rimIndexesCount = len(self._rimIndexToSegmentNodeList) - nodesCountRim = self._segments[0].getNodesCountRim() + elementsCountTransition = self._segments[0].getElementsCountTransition() + nodesCountRim = self._segments[0].getNodesCountRim() + (elementsCountTransition - 1) if self._isCore else ( + self._segments[0].getNodesCountRim()) elementsCountRim = max(1, nodesCountRim - 1) if self._rimCoordinates: self._rimNodeIds = [[None] * rimIndexesCount for _ in range(nodesCountRim)] + if self._boxCoordinates: + nodesCountAcrossMinor = self._segments[0].getCoreNodesCountAcrossMinor() + acrossMajorCounts = [segment.getElementsCountAcrossMajor() for segment in self._segments] + self._boxNodeIds = [[[None for _ in range(nodesCountAcrossMinor)] for _ in range(acrossMajorCounts[s])] + for s in range(self._segmentsCount)] + coordinates = generateData.getCoordinates() fieldcache = generateData.getFieldcache() nodes = generateData.getNodes() @@ -1051,8 +2557,37 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData): elementsCountAround = segment.getElementsCountAround() + # Create nodes + if self._boxCoordinates: + # create box nodes + bx, bd1, bd2, bd3 = (self._boxCoordinates[0], self._boxCoordinates[1], + self._boxCoordinates[2], self._boxCoordinates[3]) + for n3 in range(acrossMajorCounts[s] - 1): + for n1 in range(nodesCountAcrossMinor): + boxIndex = self._segmentNodeToBoxIndex[s][n3][n1] + segmentNodeList = self._boxIndexToSegmentNodeList[boxIndex] + nodeIdentifiersCheck = [] + for segmentNodes in segmentNodeList: + sp, n3p, n1p = segmentNodes + nodeIdentifiersCheck.append(self._boxNodeIds[sp][n3p][n1p]) + if nodeIdentifiersCheck.count(None) == len(nodeIdentifiersCheck): + nodeIdentifier = generateData.nextNodeIdentifier() + node = nodes.createNode(nodeIdentifier, nodetemplate) + fieldcache.setNode(node) + + for nodeValue, bValue in zip([Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3], + [bx, bd1, bd2, bd3]): + coordinates.setNodeParameters(fieldcache, -1, nodeValue, 1, bValue[boxIndex]) + else: + nodeIdentifier = next(id for id in nodeIdentifiersCheck if id is not None) + + self._boxNodeIds[s][n3][n1] = nodeIdentifier + + boxBoundaryNodeIds, boxBoundaryNodeToBoxId = self._createBoxBoundaryNodeIdsList(s) + if self._rimCoordinates: - # create rim nodes + # create rim nodes (including core transition nodes) for n3 in range(nodesCountRim): rx = self._rimCoordinates[0][n3] rd1 = self._rimCoordinates[1][n3] @@ -1074,6 +2609,14 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData): coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, 1, rd3[rimIndex]) layerNodeIds[rimIndex] = nodeIdentifier + # Create elements + if self._isCore: + # create box elements + self._generateBoxElements(s, n2, mesh, elementtemplate, coordinates, segment, generateData) + # create core transition elements + self._generateTransitionElements(s, n2, mesh, elementtemplate, coordinates, segment, generateData, + elementsCountAround, boxBoundaryNodeIds, boxBoundaryNodeToBoxId) + if self._rimCoordinates: # create rim elements annotationMeshGroups = generateData.getAnnotationMeshGroups(segment.getAnnotationTerms()) @@ -1091,7 +2634,8 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData): if e3 == 0: rimCoordinates = self._segments[s].getRimCoordinates(n1, n2, n3) nodeParameters.append(rimCoordinates if d3Defined else - (rimCoordinates[0], rimCoordinates[1], rimCoordinates[2], None)) + (rimCoordinates[0], rimCoordinates[1], rimCoordinates[2], + None)) nodeLayouts.append(None) for n1 in [e1, n1p]: rimIndex = self._segmentNodeToRimIndex[s][n1] @@ -1133,17 +2677,22 @@ class TubeNetworkMeshBuilder(NetworkMeshBuilder): def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSegment: float, defaultElementsCountAround: int, elementsCountThroughWall: int, layoutAnnotationGroups: list = [], annotationElementsCountsAround: list = [], - annotationElementsCountsAcrossMajor: list = []): + defaultElementsCountAcrossMajor: int = 4, elementsCountTransition: int = 1, + annotationElementsCountsAcrossMajor: list = [], isCore=False): super(TubeNetworkMeshBuilder, self).__init__( networkMesh, targetElementDensityAlongLongestSegment, layoutAnnotationGroups) self._defaultElementsCountAround = defaultElementsCountAround self._elementsCountThroughWall = elementsCountThroughWall + self._layoutAnnotationGroups = layoutAnnotationGroups self._annotationElementsCountsAround = annotationElementsCountsAround - self._annotationElementsCountsAcrossMajor = annotationElementsCountsAcrossMajor layoutFieldmodule = self._layoutRegion.getFieldmodule() self._layoutInnerCoordinates = layoutFieldmodule.findFieldByName("inner coordinates").castFiniteElement() if not self._layoutInnerCoordinates.isValid(): self._layoutInnerCoordinates = None + self._isCore = isCore + self._defaultElementsCountAcrossMajor = defaultElementsCountAcrossMajor + self._elementsCountTransition = elementsCountTransition + self._annotationElementsCountsAcrossMajor = annotationElementsCountsAcrossMajor def createSegment(self, networkSegment): pathParametersList = [get_nodeset_path_ordered_field_parameters( @@ -1154,6 +2703,9 @@ def createSegment(self, networkSegment): self._layoutNodes, self._layoutInnerCoordinates, pathValueLabels, networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions())) elementsCountAround = self._defaultElementsCountAround + elementsCountAcrossMajor = self._defaultElementsCountAcrossMajor + elementsCountAcrossMinor = (((elementsCountAround - 4) // 4 - elementsCountAcrossMajor // 2) * 2 + 6) + i = 0 for layoutAnnotationGroup in self._layoutAnnotationGroups: if i >= len(self._annotationElementsCountsAround): @@ -1163,8 +2715,27 @@ def createSegment(self, networkSegment): elementsCountAround = self._annotationElementsCountsAround[i] break i += 1 + if self._isCore: + annotationElementsCountAcrossMinor = [] + i = 0 + for layoutAnnotationGroup in self._layoutAnnotationGroups: + if i >= len(self._annotationElementsCountsAcrossMajor): + break + if self._annotationElementsCountsAcrossMajor[i] > 0: + if networkSegment.hasLayoutElementsInMeshGroup( + layoutAnnotationGroup.getMeshGroup(self._layoutMesh)): + elementsCountAcrossMajor = self._annotationElementsCountsAcrossMajor[i] + elementsCountAcrossMinor = ( + ((elementsCountAround - 4) // 4 - elementsCountAcrossMajor // 2) * 2 + 6) + annotationElementsCountAcrossMinor.append(elementsCountAcrossMinor) + break + i += 1 + # elements count across minor must be the same for all annotation groups + assert all(value == annotationElementsCountAcrossMinor[0] for value in annotationElementsCountAcrossMinor) + return TubeNetworkMeshSegment(networkSegment, pathParametersList, elementsCountAround, - self._elementsCountThroughWall) + self._elementsCountThroughWall, self._isCore, elementsCountAcrossMajor, + elementsCountAcrossMinor, self._elementsCountTransition) def createJunction(self, inSegments, outSegments): """ @@ -1322,7 +2893,7 @@ def resampleTubeCoordinates(rawTubeCoordinates, fixedElementsCountAlong=None, if not endCurveLocation: meanEndLocation += endPointLocation endCurveLocations.append(endCurveLocation) - startLength, length, endLength =\ + startLength, length, endLength = \ getCubicHermiteTrimmedCurvesLengths(cx, cd2, startCurveLocation, endCurveLocation)[0:3] sumLengths += length startLengths.append(startLength) diff --git a/tests/test_general.py b/tests/test_general.py index 3db263d3..89a4f8ee 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -1600,8 +1600,7 @@ def test_2d_tube_intersections_bifurcation(self): for networkSegment in networkSegments: pathParameters = get_nodeset_path_ordered_field_parameters( nodes, coordinates, valueLabels, networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions()) - tubeSegment = TubeNetworkMeshSegment(networkSegment, [pathParameters], - elementsCountAround=elementsCountAround, elementsCountThroughWall=1) + tubeSegment = TubeNetworkMeshSegment(networkSegment, [pathParameters], elementsCountAround, 1) tubeSegments.append(tubeSegment) trackSurfaces.append(tubeSegment.getRawTrackSurface()) diff --git a/tests/test_network.py b/tests/test_network.py index 075a5c43..660c04f4 100644 --- a/tests/test_network.py +++ b/tests/test_network.py @@ -248,13 +248,17 @@ def test_3d_tube_network_bifurcation(self): networkLayoutScaffoldPackage = settings["Network layout"] networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() self.assertTrue(networkLayoutSettings["Define inner coordinates"]) - self.assertEqual(7, len(settings)) + self.assertEqual(11, len(settings)) self.assertEqual(8, settings["Elements count around"]) self.assertEqual(1, settings["Elements count through wall"]) self.assertEqual([0], settings["Annotation elements counts around"]) self.assertEqual(4.0, settings["Target element density along longest segment"]) self.assertFalse(settings["Use linear through wall"]) self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(4, settings["Number of elements across core major"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation elements counts across major"]) context = Context("Test") region = context.getDefaultRegion() @@ -306,6 +310,69 @@ def test_3d_tube_network_bifurcation(self): self.assertAlmostEqual(outerSurfaceArea, 1.928821019338746, delta=X_TOL) self.assertAlmostEqual(innerSurfaceArea, 0.995733838660512, delta=X_TOL) + def test_3d_tube_network_bifurcation_core(self): + """ + Test bifurcation 3-D tube network with solid core is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Bifurcation") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(11, len(settings)) + self.assertEqual(8, settings["Elements count around"]) + self.assertEqual(1, settings["Elements count through wall"]) + self.assertEqual([0], settings["Annotation elements counts around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertFalse(settings["Use linear through wall"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(4, settings["Number of elements across core major"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation elements counts across major"]) + settings["Core"] = True + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + mesh3d = fieldmodule.findMeshByDimension(3) + + self.assertEqual((8 * 4 * 3) * 2 + (4 * 4 * 3), mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual((8 * 4 * 3 + 3 * 3 + 2) * 2 + (9 * 4 * 3 + 3 * 4), nodes.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + X_TOL = 1.0E-6 + + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [0.0, -0.5894427190999916, -0.10000000000000002], X_TOL) + assertAlmostEqualList(self, maximums, [2.044721359549996, 0.5894427190999916, 0.10000000000000002], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(4) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.09946683712947964, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 1.928821019338746, delta=X_TOL) + def test_3d_tube_network_sphere_cube(self): """ Test sphere cube 3-D tube network is generated correctly. @@ -315,7 +382,7 @@ def test_3d_tube_network_sphere_cube(self): networkLayoutScaffoldPackage = settings["Network layout"] networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() self.assertTrue(networkLayoutSettings["Define inner coordinates"]) - self.assertEqual(7, len(settings)) + self.assertEqual(11, len(settings)) self.assertEqual(8, settings["Elements count around"]) self.assertEqual(1, settings["Elements count through wall"]) self.assertEqual([0], settings["Annotation elements counts around"]) @@ -399,6 +466,96 @@ def test_3d_tube_network_sphere_cube(self): self.assertAlmostEqual(outerSurfaceArea, 4.045008760308934, delta=X_TOL) self.assertAlmostEqual(innerSurfaceArea, 3.3328595903228115, delta=X_TOL) + def test_3d_tube_network_sphere_cube_core(self): + """ + Test sphere cube 3-D tube network with solid core is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Sphere cube") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(11, len(settings)) + self.assertEqual(8, settings["Elements count around"]) + self.assertEqual(1, settings["Elements count through wall"]) + self.assertEqual([0], settings["Annotation elements counts around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertFalse(settings["Use linear through wall"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(4, settings["Number of elements across core major"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation elements counts across major"]) + settings["Elements count through wall"] = 2 + settings["Core"] = True + + context = Context("Test") + region = context.getDefaultRegion() + + # set custom inner coordinates + tmpRegion = region.createRegion() + networkLayoutScaffoldPackage.generate(tmpRegion) + networkMesh = networkLayoutScaffoldPackage.getConstructionObject() + functionOptions = { + "To field": {"coordinates": False, "inner coordinates": True}, + "From field": {"coordinates": True, "inner coordinates": False}, + "Mode": {"Scale": True, "Offset": False}, + "D2 value": 0.8, + "D3 value": 0.8} + editGroupName = "meshEdits" + MeshType_1d_network_layout1.assignCoordinates(tmpRegion, networkLayoutSettings, networkMesh, + functionOptions, editGroupName=editGroupName) + # put edited coordinates into scaffold package + sir = tmpRegion.createStreaminformationRegion() + srm = sir.createStreamresourceMemory() + sir.setResourceGroupName(srm, editGroupName) + sir.setResourceFieldNames(srm, ["coordinates", "inner coordinates"]) + tmpRegion.write(sir) + result, meshEditsString = srm.getBuffer() + self.assertEqual(RESULT_OK, result) + networkLayoutScaffoldPackage.setMeshEdits(meshEditsString) + + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual((8 * 3 + 4) * 4 * 12, mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual((3 * 3 + 8 * 3) * 3 * 12 + (11 * 3 + 3 * 4) * 8, nodes.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(4224, mesh2d.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + X_TOL = 1.0E-6 + + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.5664610069377635, -0.5965021612010833, -0.5985868755975445], X_TOL) + assertAlmostEqualList(self, maximums, [0.5664609474985409, 0.5965021612010833, 0.5985868966530402], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(4) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.21482044353689586, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 4.045008760308934, delta=X_TOL) + + def test_3d_tube_network_trifurcation_cross(self): """ Test trifurcation cross 3-D tube network is generated correctly with variable elements count around. @@ -408,7 +565,7 @@ def test_3d_tube_network_trifurcation_cross(self): networkLayoutScaffoldPackage = settings["Network layout"] networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() self.assertTrue(networkLayoutSettings["Define inner coordinates"]) - self.assertEqual(7, len(settings)) + self.assertEqual(11, len(settings)) self.assertEqual(8, settings["Elements count around"]) self.assertEqual(1, settings["Elements count through wall"]) self.assertEqual([0], settings["Annotation elements counts around"]) @@ -492,6 +649,97 @@ def test_3d_tube_network_trifurcation_cross(self): self.assertAlmostEqual(outerSurfaceArea, 2.59759659324524, delta=X_TOL) self.assertAlmostEqual(innerSurfaceArea, 1.3635516941376224, delta=X_TOL) + def test_3d_tube_network_trifurcation_cross_core(self): + """ + Test trifurcation cross 3-D tube network with solid coreis generated correctly with + variable elements count around. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Trifurcation cross") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(11, len(settings)) + self.assertEqual(8, settings["Elements count around"]) + self.assertEqual(1, settings["Elements count through wall"]) + self.assertEqual([0], settings["Annotation elements counts around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertFalse(settings["Use linear through wall"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(4, settings["Number of elements across core major"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation elements counts across major"]) + settings["Core"] = True + settings["Annotation elements counts around"] = [12] # requires annotation group below + settings["Annotation elements counts across major"] = [6] + + context = Context("Test") + region = context.getDefaultRegion() + + # add a user-defined annotation group to network layout to vary elements count around. Must generate first + tmpRegion = region.createRegion() + tmpFieldmodule = tmpRegion.getFieldmodule() + networkLayoutScaffoldPackage.generate(tmpRegion) + + annotationGroup1 = networkLayoutScaffoldPackage.createUserAnnotationGroup(("straight", "STRAIGHT:1")) + group = annotationGroup1.getGroup() + mesh1d = tmpFieldmodule.findMeshByDimension(1) + meshGroup = group.createMeshGroup(mesh1d) + mesh_group_add_identifier_ranges(meshGroup, [[1, 1], [4, 4]]) + self.assertEqual(2, meshGroup.getSize()) + self.assertEqual(1, annotationGroup1.getDimension()) + identifier_ranges_string = identifier_ranges_to_string(mesh_group_to_identifier_ranges(meshGroup)) + self.assertEqual("1,4", identifier_ranges_string) + networkLayoutScaffoldPackage.updateUserAnnotationGroups() + + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(1, len(annotationGroups)) + + fieldmodule = region.getFieldmodule() + + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(416, mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(569, nodes.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + # check annotation group transferred to 3D tube + annotationGroup = annotationGroups[0] + self.assertEqual("straight", annotationGroup.getName()) + self.assertEqual("STRAIGHT:1", annotationGroup.getId()) + self.assertEqual(256, annotationGroup.getMeshGroup(fieldmodule.findMeshByDimension(3)).getSize()) + + X_TOL = 1.0E-6 + + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.0447213595499958, -0.5894427190999916, -0.1], X_TOL) + assertAlmostEqualList(self, maximums, [2.044721359549996, 0.5894427190999916, 0.10000000000000002], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(4) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.1355916886131598, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 2.596646206538057, delta=X_TOL) + def test_3d_box_network_bifurcation(self): """ Test 3-D box network bifurcation is generated correctly. @@ -634,6 +882,56 @@ def test_3d_tube_network_loop(self): self.assertEqual(result, RESULT_OK) self.assertAlmostEqual(innerSurfaceArea, 0.9844505573970027, delta=1.0E-6) + def test_3d_tube_network_loop_core(self): + """ + Test loop 3-D tube network with solid core is generated correctly. + This has one segment which loops back on itself so nodes are common at start and end. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Loop") + settings = scaffoldPackage.getScaffoldSettings() + settings["Target element density along longest segment"] = 8.0 + settings["Core"] = True + + context = Context("Test") + region = context.getDefaultRegion() + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(160, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(512, mesh2d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(200, nodes.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.6, -0.6, -0.1], 1.0E-8) + assertAlmostEqualList(self, maximums, [0.6, 0.6, 0.1], 1.0E-8) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(4) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 0.09823844907582693, delta=1.0E-6) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 1.9689027258731782, delta=1.0E-6) + + def test_3d_tube_network_loop_two_segments(self): """ Test loop 3-D tube network is generated with 2 segments with fixed element boundary between them. diff --git a/tests/test_wholebody2.py b/tests/test_wholebody2.py new file mode 100644 index 00000000..7126ddbb --- /dev/null +++ b/tests/test_wholebody2.py @@ -0,0 +1,215 @@ +import math +import unittest + +from cmlibs.utils.zinc.finiteelement import evaluateFieldNodesetRange +from cmlibs.utils.zinc.general import ChangeManager + +from cmlibs.zinc.context import Context +from cmlibs.zinc.element import Element +from cmlibs.zinc.field import Field +from cmlibs.zinc.result import RESULT_OK + +from scaffoldmaker.annotation.annotationgroup import getAnnotationGroupForTerm, findAnnotationGroupByName +from scaffoldmaker.annotation.body_terms import get_body_term +from scaffoldmaker.meshtypes.meshtype_3d_wholebody2 import MeshType_3d_wholebody2 + +from testutils import assertAlmostEqualList + + +class WholeBody2ScaffoldTestCase(unittest.TestCase): + + def test_wholebody2_core(self): + """ + Test creation of Whole-body scaffold with solid core. + """ + scaffold = MeshType_3d_wholebody2 + parameterSetNames = scaffold.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default", "Human 1", "Human 2 Coarse", "Human 2 Medium", + "Human 2 Fine"]) + options = scaffold.getDefaultOptions("Default") + self.assertEqual(15, len(options)) + self.assertEqual(12, options["Number of elements around head"]) + self.assertEqual(12, options["Number of elements around torso"]) + self.assertEqual(8, options["Number of elements around arm"]) + self.assertEqual(8, options["Number of elements around leg"]) + self.assertEqual(1, options["Number of elements through wall"]) + self.assertEqual(5.0, options["Target element density along longest segment"]) + self.assertEqual(False, options["Show trim surfaces"]) + self.assertEqual(True, options["Use Core"]) + self.assertEqual(6, options["Number of elements across head core major"]) + self.assertEqual(6, options["Number of elements across torso core major"]) + self.assertEqual(4, options["Number of elements across arm core major"]) + self.assertEqual(4, options["Number of elements across leg core major"]) + self.assertEqual(1, options["Number of elements across core transition"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold.generateMesh(region, options)[0] + self.assertEqual(5, len(annotationGroups)) # Needs updating as we add more annotation groups + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(784, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(2562, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(2809, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(1032, nodes.getSize()) + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # Check coordinates range, volume + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + tol = 1.0E-6 + assertAlmostEqualList(self, minimums, [0.0, -1.976, -0.6113507244855241], tol) + assertAlmostEqualList(self, maximums, [8.748998777022969, 1.968, 0.7035640657448223], tol) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 6.419528010968665, delta=tol) + self.assertAlmostEqual(surfaceArea, 35.880031911102506, delta=tol) + + # check some annotationGroups: + expectedSizes2d = { + "skin epidermis": (308, 35.880031911102506) + } + for name in expectedSizes2d: + term = get_body_term(name) + annotationGroup = getAnnotationGroupForTerm(annotationGroups, term) + size = annotationGroup.getMeshGroup(mesh2d).getSize() + self.assertEqual(expectedSizes2d[name][0], size, name) + + surfaceMeshGroup = annotationGroup.getMeshGroup(mesh2d) + surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol) + + + def test_wholebody2_tube(self): + """ + Test creation of Whole-body scaffold without solid core. + """ + scaffold = MeshType_3d_wholebody2 + parameterSetNames = scaffold.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default", "Human 1", "Human 2 Coarse", "Human 2 Medium", + "Human 2 Fine"]) + options = scaffold.getDefaultOptions("Default") + self.assertEqual(15, len(options)) + self.assertEqual(12, options["Number of elements around head"]) + self.assertEqual(12, options["Number of elements around torso"]) + self.assertEqual(8, options["Number of elements around arm"]) + self.assertEqual(8, options["Number of elements around leg"]) + self.assertEqual(1, options["Number of elements through wall"]) + self.assertEqual(5.0, options["Target element density along longest segment"]) + self.assertEqual(False, options["Show trim surfaces"]) + self.assertEqual(True, options["Use Core"]) + self.assertEqual(6, options["Number of elements across head core major"]) + self.assertEqual(6, options["Number of elements across torso core major"]) + self.assertEqual(4, options["Number of elements across arm core major"]) + self.assertEqual(4, options["Number of elements across leg core major"]) + self.assertEqual(1, options["Number of elements across core transition"]) + + options["Use Core"] = False + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold.generateMesh(region, options)[0] + self.assertEqual(5, len(annotationGroups)) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(308, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(1254, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(1603, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(654, nodes.getSize()) + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # Check coordinates range, volume + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + tol = 1.0E-6 + assertAlmostEqualList(self, minimums, [0.0, -1.976, -0.6113507244855241], tol) + assertAlmostEqualList(self, maximums, [8.748998777022969, 1.968, 0.7035640657448223], tol) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + isExteriorXi3_0 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_0)) + isExteriorXi3_1 = fieldmodule.createFieldAnd( + isExterior, fieldmodule.createFieldIsOnFace(Element.FACE_TYPE_XI3_1)) + mesh2d = fieldmodule.findMeshByDimension(2) + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + outerSurfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) + outerSurfaceAreaField.setNumbersOfPoints(4) + result, outerSurfaceArea = outerSurfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + innerSurfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_0, coordinates, mesh2d) + innerSurfaceAreaField.setNumbersOfPoints(4) + result, innerSurfaceArea = innerSurfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 3.1692513375640607, delta=tol) + self.assertAlmostEqual(outerSurfaceArea, 35.880031911102506, delta=tol) + self.assertAlmostEqual(innerSurfaceArea, 25.851093527837623, delta=tol) + + # check some annotationGroups: + expectedSizes2d = { + "skin epidermis": (308, 35.880031911102506) + } + for name in expectedSizes2d: + term = get_body_term(name) + annotationGroup = getAnnotationGroupForTerm(annotationGroups, term) + size = annotationGroup.getMeshGroup(mesh2d).getSize() + self.assertEqual(expectedSizes2d[name][0], size, name) + + surfaceMeshGroup = annotationGroup.getMeshGroup(mesh2d) + surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol) + + +if __name__ == "__main__": + unittest.main()