diff --git a/src/scaffoldmaker/annotation/kidney_terms.py b/src/scaffoldmaker/annotation/kidney_terms.py new file mode 100644 index 00000000..f9e535a2 --- /dev/null +++ b/src/scaffoldmaker/annotation/kidney_terms.py @@ -0,0 +1,73 @@ +""" +Common resource for kidney annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +kidney_terms = [ + ("anterior surface of kidney", "UBERON:0035368", "ILX:0724840"), + ("anterior surface of left kidney", ""), + ("anterior surface of right kidney", ""), + ("cortex of kidney", "UBERON:0001225", "ILX:0726853"), + ("cortex of left kidney", "ILX:0791219"), + ("cortex of right kidney", "ILX:0791182"), + ("dorsal surface of kidney", ""), + ("dorsal surface of left kidney", ""), + ("dorsal surface of right kidney", ""), + ("hilum of kidney", "UBERON:0008716", "ILX:0731719"), + ("hilum of left kidney", ""), + ("hilum of right kidney", ""), + ("inferior pole of left kidney", "FMA:15609"), + ("inferior pole of right kidney", "FMA:15608"), + ("inner medulla of left kidney", "ILX:0784932"), + ("inner medulla of right kidney", "ILX:0791193"), + ("juxtamedullary cortex", "UBERON:0005271", "ILX:0730126"), + ("juxtamedullary cortex surface of kidney", ""), + ("juxtamedullary cortex surface of left kidney", ""), + ("juxtamedullary cortex surface of right kidney", ""), + ("kidney", "UBERON:0002113", "ILX:0735723"), + ("kidney capsule", "UBERON:0002015", "ILX:0733912"), + ("lateral edge of kidney", ""), + ("lateral edge of left kidney", ""), + ("lateral edge of right kidney", ""), + ("lateral surface of kidney", ""), + ("lateral surface of left kidney", ""), + ("lateral surface of right kidney", ""), + ("left kidney", "UBERON:0004538", "ILX:0725163"), + ("left kidney capsule", ""), + ("major calyx", "UBERON:0001226", "ILX:0730785"), + ("medial edge of kidney", ""), + ("medial edge of left kidney", ""), + ("medial edge of right kidney", ""), + ("medial surface of kidney", ""), + ("medial surface of left kidney", ""), + ("medial surface of right kidney", ""), + ("medulla of left kidney", ""), + ("medulla of right kidney", ""), + ("minor calyx", "UBERON:0001227", "ILX:0730473"), + ("outer medulla of left kidney", ""), + ("outer medulla of right kidney", ""), + ("renal medulla", "UBERON:0000362", "ILX:0729114"), + ("renal pelvis", "UBERON:0001224", "ILX:0723968"), + ("renal pyramid", "UBERON:0004200", "ILX:0727514"), + ("right kidney", "UBERON:0004539", "ILX:0735697"), + ("right kidney capsule", ""), + ("posterior surface of kidney", "UBERON:0035471", "ILX:0724479"), + ("posterior surface of left kidney", ""), + ("posterior surface of right kidney", ""), + ("superior pole of left kidney", "FMA:15607"), + ("superior pole of right kidney", "FMA:15606"), + ("ventral surface of kidney", ""), + ("ventral surface of left kidney", ""), + ("ventral surface of right kidney", "") + ] + +def get_kidney_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in kidney_terms: + if name in term: + return (term[0], term[1]) + raise NameError("Kidney annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_kidney1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_kidney1.py new file mode 100644 index 00000000..be491615 --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_kidney1.py @@ -0,0 +1,802 @@ +""" +Generates a 3D kidney using tube network mesh. +""" +import math + +from cmlibs.maths.vectorops import mult, set_magnitude, cross +from cmlibs.utils.zinc.field import find_or_create_field_coordinates, findOrCreateFieldCoordinates +from cmlibs.zinc.field import Field + +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, \ + getAnnotationGroupForTerm, findAnnotationGroupByName +from scaffoldmaker.annotation.kidney_terms import get_kidney_term +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.utils.interpolation import smoothCubicHermiteDerivativesLine, sampleCubicHermiteCurves, \ + smoothCurveSideCrossDerivatives +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from scaffoldmaker.utils.networkmesh import NetworkMesh +from scaffoldmaker.utils.tubenetworkmesh import KidneyTubeNetworkMeshBuilder, TubeNetworkMeshGenerateData +from cmlibs.zinc.node import Node + +from scaffoldmaker.utils.zinc_utils import translate_nodeset_coordinates + + +class MeshType_1d_kidney_network_layout1(MeshType_1d_network_layout1): + """ + Defines kidney network layout. + """ + + showKidneys = [False, False] + + @classmethod + def getName(cls): + return "1D Kidney Network Layout 1" + + @classmethod + def getParameterSetNames(cls): + return ["Default"] + + @classmethod + def getDefaultOptions(cls, parameterSetName="Default"): + options = {} + options["Base parameter set"] = "Human 1" if (parameterSetName == "Default") else parameterSetName + options["Define inner coordinates"] = True + options["Left kidney"] = True + options["Right kidney"] = True + options["Kidney length"] = 1.0 + options["Kidney width"] = 0.5 + options["Kidney thickness"] = 0.4 + options["Medulla to cortex proportion"] = 0.6 + return options + + @classmethod + def getOrderedOptionNames(cls): + return [ + "Left kidney", + "Right kidney", + "Kidney length", + "Kidney width", + "Kidney thickness", + "Medulla to cortex proportion" + ] + + @classmethod + def checkOptions(cls, options): + dependentChanges = False + for key in [ + "Kidney length", + "Kidney width", + "Kidney thickness" + ]: + if options[key] < 0.1: + options[key] = 0.1 + + if options["Medulla to cortex proportion"] < 0.1: + options["Medulla to cortex proportion"] = 0.1 + elif options["Medulla to cortex proportion"] > 0.9: + options["Medulla to cortex proportion"] = 0.9 + + if not options["Left kidney"] and not options["Right kidney"]: + dependentChanges = True + options["Left kidney"] = True + + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options): + """ + Generate the unrefined mesh. + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: [] empty list of AnnotationGroup, NetworkMesh + """ + # parameters + structure = options["Structure"] = cls.getLayoutStructure(options) + isLeftKidney = options["Left kidney"] + isRightKidney = options["Right kidney"] + kidneyLength = options["Kidney length"] + halfKidneyLength = 0.5 * kidneyLength + halfKidneyWidth = 0.5 * options["Kidney width"] + halfKidneyThickness = 0.5 * options["Kidney thickness"] + innerProportionDefault = options["Medulla to cortex proportion"] + cls.setShowKidneys(options) + + networkMesh = NetworkMesh(structure) + networkMesh.create1DLayoutMesh(region) + + fieldmodule = region.getFieldmodule() + mesh = fieldmodule.findMeshByDimension(1) + + # set up element annotations + kidneyGroup = AnnotationGroup(region, get_kidney_term("kidney")) + kidneyMeshGroup = kidneyGroup.getMeshGroup(mesh) + + leftKidneyGroup = AnnotationGroup(region, get_kidney_term("left kidney")) + leftKidneyMeshGroup = leftKidneyGroup.getMeshGroup(mesh) + + rightKidneyGroup = AnnotationGroup(region, get_kidney_term("right kidney")) + rightKidneyMeshGroup = rightKidneyGroup.getMeshGroup(mesh) + + annotationGroups = [kidneyGroup, leftKidneyGroup, rightKidneyGroup] + meshGroups = [kidneyMeshGroup, leftKidneyMeshGroup, rightKidneyMeshGroup] + + # set coordinates (outer) + fieldcache = fieldmodule.createFieldcache() + coordinates = find_or_create_field_coordinates(fieldmodule) + # need to ensure inner coordinates are at least defined: + cls.defineInnerCoordinates(region, coordinates, options, networkMesh, innerProportion=innerProportionDefault) + innerCoordinates = find_or_create_field_coordinates(fieldmodule, "inner coordinates") + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + + # Kidney + nodeIdentifier = 1 + elementIdentifier = 1 + kidneyElementsCount = 2 + capRadius = cls.getCapRadius(halfKidneyWidth, halfKidneyThickness) * ( + halfKidneyWidth * 0.45 + halfKidneyThickness * 0.55) + extensionLength = 0.5 * (halfKidneyWidth * 0.45 + halfKidneyThickness * 0.55) + halfLayoutLength = (halfKidneyLength - capRadius - extensionLength) + kidneyScale = 2 * halfLayoutLength / kidneyElementsCount + + leftKidney, rightKidney = 0, 1 + kidneys = [kidney for show, kidney in [(isLeftKidney, leftKidney), (isRightKidney, rightKidney)] if show] + for kidney in kidneys: + mx = [0.0, 0.0, 0.0] + d1 = [kidneyScale, 0.0, 0.0] + d3 = [0.0, 0.0, halfKidneyThickness] + id3 = mult(d3, innerProportionDefault) + + tx = halfLayoutLength + sx = [-tx, 0.0, 0.0] if kidney is leftKidney else [-tx, 0.0, 0.0] + ex = [tx, 0.0, 0.0] if kidney is leftKidney else [tx, 0.0, 0.0] + sd1 = mult([-1.0, 0.0, 0.0], kidneyScale) + ed1 = [-sd1[0], sd1[1], sd1[2]] + nx, nd1 = sampleCubicHermiteCurves([sx, mx, ex], [sd1, d1, ed1], kidneyElementsCount)[0:2] + nd1 = smoothCubicHermiteDerivativesLine(nx, nd1) + + sd2_list = [] + sd3_list = [] + sNodeIdentifiers = [] + for e in range(kidneyElementsCount + 1): + sNodeIdentifiers.append(nodeIdentifier) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + sd2 = set_magnitude(cross(d3, nd1[e]), halfKidneyWidth) + sid2 = mult(sd2, innerProportionDefault) + sd2_list.append(sd2) + sd3_list.append(d3) + for field, derivatives in ((coordinates, (nd1[e], sd2, d3)), (innerCoordinates, (nd1[e], sid2, id3))): + setNodeFieldParameters(field, fieldcache, nx[e], *derivatives) + nodeIdentifier += 1 + + sd12 = smoothCurveSideCrossDerivatives(nx, nd1, [sd2_list])[0] + sd13 = smoothCurveSideCrossDerivatives(nx, nd1, [sd3_list])[0] + for e in range(kidneyElementsCount + 1): + node = nodes.findNodeByIdentifier(sNodeIdentifiers[e]) + fieldcache.setNode(node) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, sd12[e]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, sd13[e]) + sid12 = mult(sd12[e], innerProportionDefault) + sid13 = mult(sd13[e], innerProportionDefault) + innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, sid12) + innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, sid13) + + # add annotations + for e in range(kidneyElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + meshGroups[0].addElement(element) + if kidney is leftKidney and isLeftKidney: + meshGroups[1].addElement(element) + if kidney is rightKidney and isRightKidney: + meshGroups[2].addElement(element) + elementIdentifier += 1 + + return annotationGroups, networkMesh + + @classmethod + def getLayoutStructure(cls, options): + """ + Generate 1D layout structure based on the number of elements count along. + :param options: Dict containing options. See getDefaultOptions(). + :return: string version of the 1D layout structure + """ + nodes_count = 3 + assert nodes_count > 1 + + left = f"({'-'.join(map(str, range(1, nodes_count + 1)))})" + if options["Left kidney"] and options["Right kidney"]: + right = f"({'-'.join(map(str, range(nodes_count + 1, 2 * nodes_count + 1)))})" + return f"{left},{right}" + return left + + @classmethod + def getCapRadius(cls, majorRadius, minorRadius): + """ + Calculate the radius of the cap mesh based on the major radius and the minor radius of tube cross-section. + :param majorRadius: The radius of a tube in the major-axis. + :param minorRadius: The radius of a tube in the minor-axis. + :return: Cap radius + """ + if majorRadius > minorRadius: + return math.pow((majorRadius / minorRadius), 1 / 3) + elif majorRadius < minorRadius: + return math.pow((minorRadius / majorRadius), 1 / 3) + else: + return majorRadius + + @classmethod + def getShowKidneys(cls): + return cls.showKidneys + + @classmethod + def setShowKidneys(cls, options): + cls.showKidneys[0] = True if options["Left kidney"] else False + cls.showKidneys[1] = True if options["Right kidney"] else False + + +class MeshType_3d_kidney1(Scaffold_base): + """ + Generates a 3-D Kidney. + """ + + @classmethod + def getName(cls): + return "3D Kidney 1" + + @classmethod + def getParameterSetNames(cls): + return [ + "Default", + "Human 1" + ] + + @classmethod + def getDefaultOptions(cls, parameterSetName='Default'): + options = {} + useParameterSetName = "Human 1" if (parameterSetName == "Default") else parameterSetName + options["Base parameter set"] = useParameterSetName + options["Kidney network layout"] = ScaffoldPackage(MeshType_1d_kidney_network_layout1) + options["Number of elements around"] = 8 + options["Number of elements through cortex"] = 1 + options["Annotation numbers of elements around"] = [0] + options["Target element density along longest segment"] = 2.0 + options["Number of elements across core box minor"] = 2 + options["Number of elements across core transition"] = 1 + options["Annotation numbers of elements across core box minor"] = [0] + options["Kidney spacing"] = 1.0 + options["Kidney curvature"] = 1.0 + options["Refine"] = False + options["Refine number of elements"] = 4 + return options + + @classmethod + def getOrderedOptionNames(cls): + optionNames = [ + "Kidney network layout", + "Number of elements around", + "Number of elements through cortex", + "Annotation numbers of elements around", + "Target element density along longest segment", + "Number of elements across core box minor", + "Number of elements across core transition", + "Annotation numbers of elements across core box minor", + "Kidney spacing", + "Kidney curvature", + "Refine", + "Refine number of elements" + ] + return optionNames + + @classmethod + def getOptionValidScaffoldTypes(cls, optionName): + if optionName == "Kidney network layout": + return [MeshType_1d_kidney_network_layout1] + return [] + + @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 == "Kidney network layout": + if not parameterSetName: + parameterSetName = "Default" + return ScaffoldPackage(MeshType_1d_kidney_network_layout1, defaultParameterSetName=parameterSetName) + assert False, cls.__name__ + ".getOptionScaffoldPackage: Option " + optionName + " is not a scaffold" + + @classmethod + def checkOptions(cls, options): + dependentChanges = False + if (options["Kidney network layout"].getScaffoldType() not in + cls.getOptionValidScaffoldTypes("Kidney network layout")): + options["Kidney network layout"] = ScaffoldPackage(MeshType_1d_kidney_network_layout1) + + if options["Number of elements around"] < 8: + options["Number of elements around"] = 8 + elif options["Number of elements around"] % 4: + options["Number of elements around"] += 4 - (options["Number of elements around"] % 4) + + if options["Number of elements through cortex"] < 1: + options["Number of elements through cortex"] = 1 + + if options["Number of elements across core transition"] < 1: + options["Number of elements across core transition"] = 1 + + minElementsCountAround = options["Number of elements around"] + maxElementsCountCoreBoxMinor = minElementsCountAround // 2 - 2 + if options["Number of elements across core box minor"] < 2: + options["Number of elements across core box minor"] = 2 + elif options["Number of elements across core box minor"] > maxElementsCountCoreBoxMinor: + options["Number of elements across core box minor"] = maxElementsCountCoreBoxMinor + dependentChanges = True + elif options["Number of elements across core box minor"] % 2: + options["Number of elements across core box minor"] += options["Number of elements across core box minor"] % 2 + + annotationElementsCountsAround = options["Annotation numbers of elements around"] + if len(annotationElementsCountsAround) == 0: + options["Annotation numbers of elements around"] = [0] + else: + for i in range(len(annotationElementsCountsAround)): + if annotationElementsCountsAround[i] <= 0: + annotationElementsCountsAround[i] = 0 + else: + if annotationElementsCountsAround[i] < 8: + annotationElementsCountsAround[i] = 8 + elif annotationElementsCountsAround[i] % 4: + annotationElementsCountsAround[i] += 4 - (annotationElementsCountsAround[i] % 4) + if annotationElementsCountsAround[i] < minElementsCountAround: + minElementsCountAround = annotationElementsCountsAround[i] + + annotationCoreBoxMinorCounts = options["Annotation numbers of elements across core box minor"] + if len(annotationCoreBoxMinorCounts) == 0: + annotationCoreBoxMinorCounts = options["Annotation numbers of elements across core box minor"] = [0] + if len(annotationCoreBoxMinorCounts) > len(annotationElementsCountsAround): + annotationCoreBoxMinorCounts = options["Annotation numbers of elements across core box minor"] = \ + annotationCoreBoxMinorCounts[:len(annotationElementsCountsAround)] + dependentChanges = True + for i in range(len(annotationCoreBoxMinorCounts)): + aroundCount = annotationElementsCountsAround[i] if annotationElementsCountsAround[i] \ + else options["Number of elements around"] + maxCoreBoxMinorCount = aroundCount // 2 - 2 + if annotationCoreBoxMinorCounts[i] <= 0: + annotationCoreBoxMinorCounts[i] = 0 + # this may reduce the default + if maxCoreBoxMinorCount < options["Number of elements across core box minor"]: + options["Number of elements across core box minor"] = maxCoreBoxMinorCount + dependentChanges = True + elif annotationCoreBoxMinorCounts[i] < 2: + annotationCoreBoxMinorCounts[i] = 2 + elif annotationCoreBoxMinorCounts[i] > maxCoreBoxMinorCount: + annotationCoreBoxMinorCounts[i] = maxCoreBoxMinorCount + dependentChanges = True + elif annotationCoreBoxMinorCounts[i] % 2: + annotationCoreBoxMinorCounts[i] += 1 + + if options["Target element density along longest segment"] < 2.0: + options["Target element density along longest segment"] = 2.0 + + if options['Refine number of elements'] < 1: + options['Refine number of elements'] = 1 + 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 + """ + networkLayout = options["Kidney network layout"] + layoutRegion = region.createRegion() + networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters + layoutAnnotationGroups = networkLayout.getAnnotationGroups() + networkMesh = networkLayout.getConstructionObject() + showKidneys = getShowKidneysSettings() + isLeftKidney = showKidneys[0] + isRightKidney = showKidneys[1] + + kidneyTubeNetworkMeshBuilder = KidneyTubeNetworkMeshBuilder( + networkMesh, + targetElementDensityAlongLongestSegment=options["Target element density along longest segment"], + defaultElementsCountAround=options["Number of elements around"], + elementsCountThroughShell=options["Number of elements through cortex"], + layoutAnnotationGroups=layoutAnnotationGroups, + isCore=True, + elementsCountTransition=options["Number of elements across core transition"], + defaultElementsCountCoreBoxMinor=options["Number of elements across core box minor"], + annotationElementsCountsCoreBoxMinor=options["Annotation numbers of elements across core box minor"], + showKidneys=showKidneys + ) + + kidneyTubeNetworkMeshBuilder.build() + generateData = TubeNetworkMeshGenerateData( + region, 3, + isLinearThroughShell=False) + kidneyTubeNetworkMeshBuilder.generateMesh(generateData) + annotationGroups = generateData.getAnnotationGroups() + + # add kidney-specific annotation groups + fm = region.getFieldmodule() + coordinates = findOrCreateFieldCoordinates(fm) + mesh = generateData.getMesh() + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + + coreGroup = getAnnotationGroupForTerm(annotationGroups, ("core", "")).getGroup() + shellGroup = getAnnotationGroupForTerm(annotationGroups, ("shell", "")).getGroup() + openingGroup = getAnnotationGroupForTerm(annotationGroups, ("opening", "")).getGroup() + + kidneyGroup = AnnotationGroup(region, get_kidney_term("kidney")) + kidneyNodesetGroup = kidneyGroup.getNodesetGroup(nodes) + + leftKidneyGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("left kidney")) + leftKidneyNodesetGroup = leftKidneyGroup.getNodesetGroup(nodes) + + rightKidneyGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("right kidney")) + rightKidneyNodesetGroup = rightKidneyGroup.getNodesetGroup(nodes) + + hilumGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("hilum of kidney")) + hilumGroup.getMeshGroup(mesh).addElementsConditional(openingGroup) + + cortexGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("cortex of kidney")) + tempGroup = fm.createFieldSubtract(shellGroup, openingGroup) + cortexGroup.getMeshGroup(mesh).addElementsConditional(tempGroup) + + medullaGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("renal medulla")) + medullaGroup.getMeshGroup(mesh).addElementsConditional(coreGroup) + + for term in ["core", "shell", "opening"]: + annotationGroups.remove(findAnnotationGroupByName(annotationGroups, term)) + + # marker points + leftSuperiorPoleGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("superior pole of left kidney")) + leftInferiorPoleGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("inferior pole of left kidney")) + + rightSuperiorPoleGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("superior pole of right kidney")) + rightInferiorPoleGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("inferior pole of right kidney")) + + markerList = [] + elementsCountAround = options["Number of elements around"] + elementsCountAlong = int(options["Target element density along longest segment"] + 2) # extra sections from the cap mesh + elementsCountThroughShell = options["Number of elements through cortex"] + elementsCountCoreBoxMinor = options["Number of elements across core box minor"] + elementsCountCoreBoxMajor = (elementsCountAround // 2) - elementsCountCoreBoxMinor + elementsCountTransition = options["Number of elements across core transition"] + + box_count = elementsCountCoreBoxMinor * elementsCountCoreBoxMajor + depth = elementsCountThroughShell + elementsCountTransition + offset = elementsCountCoreBoxMajor // 2 * elementsCountCoreBoxMinor + elementsCountCoreBoxMinor // 2 + 1 + cap_count = box_count * (depth + 1) + elementsCountAround * depth + tube_section_count = box_count + elementsCountAround * depth + tube_count = tube_section_count * elementsCountAlong + kidney_elements_count = cap_count * 2 + tube_count if showKidneys[0] else 0 + + if isLeftKidney: + idx = box_count * depth + offset + markerList.append({"group": leftSuperiorPoleGroup, "elementId": idx, "xi": [0.0, 0.0, 1.0]}) + + idx = cap_count + tube_count + (box_count * depth + offset) + markerList.append({"group": leftInferiorPoleGroup, "elementId": idx, "xi": [0.0, 1.0, 1.0]}) + + if isRightKidney: + idx = kidney_elements_count + box_count * depth + offset + markerList.append({"group": rightSuperiorPoleGroup, "elementId": idx, "xi": [0.0, 0.0, 1.0]}) + + idx = kidney_elements_count + cap_count + tube_count + (box_count * depth + offset) + markerList.append({"group": rightInferiorPoleGroup, "elementId": idx, "xi": [0.0, 1.0, 1.0]}) + + nodeIdentifier = generateData.nextNodeIdentifier() + for marker in markerList: + annotationGroup = marker["group"] + markerNode = annotationGroup.createMarkerNode( + nodeIdentifier, element=mesh.findElementByIdentifier(marker["elementId"]), xi=marker["xi"]) + annotationGroup.setMarkerMaterialCoordinates(coordinates) + kidneyNodesetGroup.addNode(markerNode) + nodeIdentifier += 1 + + # transformation + leftKidney, rightKidney = 0, 1 + kidneys = [kidney for show, kidney in [(isLeftKidney, leftKidney), (isRightKidney, rightKidney)] if show] + for kidney in kidneys: + isLeft = True if kidney == leftKidney else False + isRight = True if kidney == rightKidney else False + spacing = options["Kidney spacing"] / 2 if isLeft else -options["Kidney spacing"] / 2 + curvature = -options["Kidney curvature"] if isLeft else options["Kidney curvature"] + + kidneyNodeset = leftKidneyNodesetGroup if isLeft else rightKidneyNodesetGroup + + if curvature != 0.0: + if isLeft: + bendKidneyMeshAroundZAxis(curvature, fm, coordinates, kidneyNodeset, stationaryPointXY=[0.05, 0.0]) + if isRight: + bendKidneyMeshAroundZAxis(curvature, fm, coordinates, kidneyNodeset, stationaryPointXY=[-0.05, 0.0]) + + translate_nodeset_coordinates(kidneyNodeset, coordinates, [0.0, spacing, 0.0]) + + return annotationGroups, None + + + @classmethod + def refineMesh(cls, meshRefinement, options): + """ + Refine source mesh into separate region, with change of basis. + :param meshRefinement: MeshRefinement, which knows source and target region. + :param options: Dict containing options. See getDefaultOptions(). + """ + assert isinstance(meshRefinement, MeshRefinement) + refineElementsCount = options['Refine number of elements'] + meshRefinement.refineAllElementsCubeStandard3d(refineElementsCount, refineElementsCount, refineElementsCount) + + + @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. + """ + show_kidneys = getShowKidneysSettings() + + # Initialize field module and meshes + fm = region.getFieldmodule() + mesh1d = fm.findMeshByDimension(1) + mesh2d = fm.findMeshByDimension(2) + mesh3d = fm.findMeshByDimension(3) + is_exterior = fm.createFieldIsExterior() + + # Get base kidney group + kidney_group = getAnnotationGroupForTerm(annotationGroups, get_kidney_term("kidney")).getGroup() + + # Create side groups and tracking dictionaries + side_groups = {} + side_kidney_groups = {"left": {}, "right": {}} + + for side in ["left", "right"]: + group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(f"{side} kidney")) + side_groups[side] = group.getGroup() + side_kidney_groups[side][f"{side} kidney"] = group + + # Create kidney part groups (cortex, hilum, medulla) + create_kidney_part_groups(fm, mesh3d, annotationGroups, region, side_groups, side_kidney_groups) + + # Create capsule groups + create_capsule_groups(fm, mesh2d, annotationGroups, region, kidney_group, side_groups, side_kidney_groups, is_exterior) + + # Create surface and edge annotation groups + create_surface_and_edge_groups(fm, mesh2d, mesh1d, annotationGroups, region, side_groups, side_kidney_groups, is_exterior) + + # Remove groups based on kidney visibility settings + remove_hidden_kidney_groups(show_kidneys, side_kidney_groups, annotationGroups) + + +def getShowKidneysSettings(): + return MeshType_1d_kidney_network_layout1.getShowKidneys() + + +def create_kidney_part_groups(fm, mesh3d, annotationGroups, region, side_groups, side_kidney_groups): + """ + Create cortex, hilum, and medulla groups for each kidney side. + """ + kidney_parts = ["cortex", "hilum", "medulla"] + + for part in kidney_parts: + # Get the anatomical term for the part + arb_term = f"renal {part}" if part == "medulla" else f"{part} of kidney" + arb_group = getAnnotationGroupForTerm(annotationGroups, get_kidney_term(arb_term)).getGroup() + + # Create side-specific part groups + for side in ["left", "right"]: + part_term = f"{part} of {side} kidney" + part_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(part_term)) + part_group.getMeshGroup(mesh3d).addElementsConditional(fm.createFieldAnd(arb_group, side_groups[side])) + side_kidney_groups[side][part_term] = part_group + + +def create_capsule_groups(fm, mesh2d, annotationGroups, region, kidney_group, side_groups, side_kidney_groups, is_exterior): + """ + Create kidney capsule surface groups. + """ + # General kidney capsule group + kidney_capsule_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term("kidney capsule")) + kidney_exterior = fm.createFieldAnd(kidney_group, is_exterior) + kidney_capsule_group.getMeshGroup(mesh2d).addElementsConditional(kidney_exterior) + + # Side-specific capsule groups + for side in ["left", "right"]: + capsule_term = f"{side} kidney capsule" + capsule_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(capsule_term)) + capsule_exterior = fm.createFieldAnd(side_groups[side], is_exterior) + capsule_group.getMeshGroup(mesh2d).addElementsConditional(capsule_exterior) + side_kidney_groups[side][capsule_term] = capsule_group + + +def create_surface_and_edge_groups(fm, mesh2d, mesh1d, annotationGroups, region, side_groups, side_kidney_groups, is_exterior): + """ + Create surface and edge annotation groups. + """ + surface_types = ["anterior", "posterior", "lateral", "medial", "dorsal", "ventral", "juxtamedullary cortex"] + surface_fields = {} + + # Create surface groups + for surface_type in surface_types: + if surface_type == "juxtamedullary cortex": + cortex_group = getAnnotationGroupForTerm(annotationGroups, get_kidney_term("cortex of kidney")).getGroup() + medulla_group = getAnnotationGroupForTerm(annotationGroups, get_kidney_term("renal medulla")).getGroup() + surface_exterior = fm.createFieldAnd(medulla_group, cortex_group) + else: + base_group = getAnnotationGroupForTerm(annotationGroups, (surface_type, "")) + surface_field = base_group.getGroup() + surface_exterior = fm.createFieldAnd(surface_field, is_exterior) + surface_fields[surface_type] = surface_field + + # General kidney surface group + general_term = f"{surface_type} surface of kidney" + general_surface_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(general_term)) + general_surface_group.getMeshGroup(mesh2d).addElementsConditional(surface_exterior) + + # Side-specific surface groups + for side in ["left", "right"]: + side_term = f"{surface_type} surface of {side} kidney" + side_surface_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(side_term)) + side_surface_field = fm.createFieldAnd(surface_exterior, side_groups[side]) + side_surface_group.getMeshGroup(mesh2d).addElementsConditional(side_surface_field) + side_kidney_groups[side][side_term] = side_surface_group + + # Create edge groups at dorsal-ventral intersection + dorsal_ventral_border = fm.createFieldAnd( + fm.createFieldAnd(surface_fields["dorsal"], surface_fields["ventral"]), is_exterior) + + edge_types = ["lateral", "medial"] + + for edge_type in edge_types: + # General edge group + edge_term = f"{edge_type} edge of kidney" + edge_field = fm.createFieldAnd(surface_fields[edge_type], dorsal_ventral_border) + general_edge_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(edge_term)) + general_edge_group.getMeshGroup(mesh1d).addElementsConditional(edge_field) + + # Side-specific edge groups + for side in ["left", "right"]: + side_edge_term = f"{edge_type} edge of {side} kidney" + side_edge_group = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_kidney_term(side_edge_term)) + side_edge_field = fm.createFieldAnd(edge_field, side_groups[side]) + side_edge_group.getMeshGroup(mesh1d).addElementsConditional(side_edge_field) + side_kidney_groups[side][side_edge_term] = side_edge_group + + +def remove_hidden_kidney_groups(show_kidneys, side_kidney_groups, annotationGroups): + """ + Remove annotation groups for kidneys that should not be shown. + """ + if not show_kidneys[0]: # Left kidney + for group in side_kidney_groups["left"].values(): + annotationGroups.remove(group) + + if not show_kidneys[1]: # Right kidney + for group in side_kidney_groups["right"].values(): + annotationGroups.remove(group) + + +def bendKidneyMeshAroundZAxis(curvature, fm, coordinates, kidneyNodeset, stationaryPointXY): + """ + Transform coordinates by bending with curvature about a centre point the radius in + x direction from stationaryPointXY. + :param curvature: 1/radius. Must be non-zero. + :param fm: Field module being worked with. + :param coordinates: The coordinate field, initially circular in y-z plane. + :param kidneyNodeset: Zinc NodesetGroup containing nodes to transform. + :param stationaryPointXY: Coordinates x, y which are not displaced by bending. + """ + rotateKidneyMeshAboutZAxis(90, fm, coordinates, kidneyNodeset) + + radius = 1.0 / curvature + scale = fm.createFieldConstant([-1.0, -curvature, -1.0]) + centreOffset = [stationaryPointXY[0] - radius, stationaryPointXY[1], 0.0] + centreOfCurvature = fm.createFieldConstant(centreOffset) + polarCoordinates = (centreOfCurvature - coordinates) * scale + polarCoordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_CYLINDRICAL_POLAR) + rcCoordinates = fm.createFieldCoordinateTransformation(polarCoordinates) + rcCoordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN) + newCoordinates = rcCoordinates + centreOfCurvature + + fieldassignment = coordinates.createFieldassignment(newCoordinates) + fieldassignment.setNodeset(kidneyNodeset) + fieldassignment.assign() + + rotateKidneyMeshAboutZAxis(-90, fm, coordinates, kidneyNodeset) + + +def rotateKidneyMeshAboutZAxis(rotateAngle, fm, coordinates, kidneyNodeset): + """ + Rotates the mesh coordinates about a specified axis using the right-hand rule. + :param rotateAngle: Angle of rotation in degrees. + :param fm: Field module being worked with. + :param coordinates: The coordinate field, initially circular in y-z plane. + :param kidneyNodeset: Zinc NodesetGroup containing nodes to transform. + :return: None + """ + rotateAngle = -math.radians(rotateAngle) # negative value due to right handed rule + rotateMatrix = fm.createFieldConstant([math.cos(rotateAngle), math.sin(rotateAngle), 0.0, + -math.sin(rotateAngle), math.cos(rotateAngle), 0.0, + 0.0, 0.0, 1.0]) + + rotated_coordinates = fm.createFieldMatrixMultiply(3, rotateMatrix, coordinates) + + fieldassignment = coordinates.createFieldassignment(rotated_coordinates) + fieldassignment.setNodeset(kidneyNodeset) + fieldassignment.assign() + + +def setNodeFieldParameters(field, fieldcache, x, d1, d2, d3, d12=None, d13=None): + """ + Assign node field parameters x, d1, d2, d3 of field. + :param field: Field parameters to assign. + :param fieldcache: Fieldcache with node set. + :param x: Parameters to set for Node.VALUE_LABEL_VALUE. + :param d1: Parameters to set for Node.VALUE_LABEL_D_DS1. + :param d2: Parameters to set for Node.VALUE_LABEL_D_DS2. + :param d3: Parameters to set for Node.VALUE_LABEL_D_DS3. + :param d12: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS2. + :param d13: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS3. + :return: + """ + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, x) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, 1, d3) + if d12: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d12) + if d13: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d13) + + +def setNodeFieldVersionDerivatives(field, fieldcache, version, d1, d2, d3, d12=None, d13=None): + """ + Assign node field parameters d1, d2, d3 of field. + :param field: Field to assign parameters of. + :param fieldcache: Fieldcache with node set. + :param version: Version of d1, d2, d3 >= 1. + :param d1: Parameters to set for Node.VALUE_LABEL_D_DS1. + :param d2: Parameters to set for Node.VALUE_LABEL_D_DS2. + :param d3: Parameters to set for Node.VALUE_LABEL_D_DS3. + :param d12: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS2. + :param d13: Optional parameters to set for Node.VALUE_LABEL_D2_DS1DS3. + :return: + """ + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, version, d1) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, version, d2) + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, version, d3) + if d12: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, d12) + if d13: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, d13) + + +def setHilumGroupThreshold(fm, coordinates, halfLength): + """ + Creates a field to identify hilum elements based on x-coordinate threshold. + :param fm: Field module used for creating and managing fields. + :param coordinates: The coordinate field. + :param halfLength: Half-length of tube. + :return is_within_threshold: True for elements between the positive and negative x-threshold. + """ + x_component = fm.createFieldComponent(coordinates, [1]) + x_threshold = 0.25 * halfLength + minus_one = fm.createFieldConstant(-1) + + x_threshold_field = fm.createFieldConstant(x_threshold) + is_less_than_threshold = fm.createFieldLessThan(x_component, x_threshold_field) + is_greater_than_threshold = fm.createFieldGreaterThan(x_component, x_threshold_field * minus_one) + is_within_threshold = fm.createFieldAnd(is_less_than_threshold, is_greater_than_threshold) + + return is_within_threshold diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index f6bd8a5b..2489398f 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -44,6 +44,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_nerve1 import MeshType_3d_nerve1 from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1 from scaffoldmaker.meshtypes.meshtype_3d_ostium2 import MeshType_3d_ostium2 +from scaffoldmaker.meshtypes.meshtype_3d_kidney1 import MeshType_1d_kidney_network_layout1, MeshType_3d_kidney1 from scaffoldmaker.meshtypes.meshtype_3d_smallintestine1 import MeshType_3d_smallintestine1 from scaffoldmaker.meshtypes.meshtype_3d_solidcylinder1 import MeshType_3d_solidcylinder1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 @@ -103,6 +104,7 @@ class Scaffolds(object): MeshType_3d_heartventricles3, MeshType_3d_heartventriclesbase1, MeshType_3d_heartventriclesbase2, + MeshType_3d_kidney1, MeshType_3d_lens1, MeshType_3d_lung1, MeshType_3d_lung2, @@ -133,6 +135,7 @@ class Scaffolds(object): MeshType_1d_human_body_network_layout1, MeshType_1d_human_spinal_nerve_network_layout1, MeshType_1d_human_trigeminal_nerve_network_layout1, + MeshType_1d_kidney_network_layout1, MeshType_1d_uterus_network_layout1 ] diff --git a/src/scaffoldmaker/utils/capmesh.py b/src/scaffoldmaker/utils/capmesh.py new file mode 100644 index 00000000..7c19501f --- /dev/null +++ b/src/scaffoldmaker/utils/capmesh.py @@ -0,0 +1,1854 @@ +""" +Specialisation of Tube Network Mesh for building 3-D cap mesh. +""" + +import math + +from cmlibs.maths.vectorops import magnitude, sub, add, set_magnitude, normalize, rotate_vector_around_vector, cross, \ + angle, mult, div +from cmlibs.zinc.element import Element +from cmlibs.zinc.node import Node +from scaffoldmaker.utils.eft_utils import determineCubicHermiteSerendipityEft +from scaffoldmaker.utils.interpolation import smoothCubicHermiteDerivativesLine, sampleCubicHermiteCurves, \ + interpolateSampleCubicHermite, DerivativeScalingMode +from scaffoldmaker.utils.spheremesh import calculate_arc_length, local_to_global_coordinates, spherical_to_cartesian, \ + calculate_azimuth + + +class CapMesh: + """ + Cap mesh generator. + """ + + def __init__(self, elementsCountAround, elementsCountCoreBoxMajor, elementsCountCoreBoxMinor, + elementsCountThroughShell, elementsCountTransition, networkPathParameters, isCap, isCore): + """ + :param elementsCountAround: Number of elements around this segment. + :param elementsCountCoreBoxMajor: Number of elements across core box major axis. + :param elementsCountCoreBoxMinor: Number of elements across core box minor axis. + :param elementsCountThroughShell: Number of elements between inner and outer tube. + :param elementsCountTransition: Number of elements across transition zone between core box elements and + rim elements. + :param networkPathParameters: List containing path parameters of a tube network. + :param isCap: List [startCap, endCap] with boolean values. True if the tube segment requires a cap at the + start of a segment, or at the end of a segment, respectively. [True, True] if the segment requires cap at both + ends. + :param isCore: True for tube network with a solid core, False for regular tube network. + """ + self._isCap = isCap + self._isCore = isCore + + self._elementsCountAround = elementsCountAround + self._elementsCountCoreBoxMajor = elementsCountCoreBoxMajor + self._elementsCountCoreBoxMinor = elementsCountCoreBoxMinor + self._elementsCountThroughShell = elementsCountThroughShell + self._elementsCountTransition = elementsCountTransition + + self._networkPathParameters = networkPathParameters + self._tubeBoxCoordinates = None # tube box coordinates + self._tubeTransitionCoordinates = None # tube transition coordinates + self._tubeShellCoordinates = None # tube rim coordinates + + self._isStartCap = None + self._generateData = None + + self._boxExtCoordinates = None + # coordinates and derivatives for box nodes extended from the tube segment + # list[startCap, endCap][x, d1, d2, d3][nAcrossMajor][nAcrossMinor] + self._transitionExtCoordinates = None + self._shellExtCoordinates = None + # coordinates and derivatives for shell nodes extended from the tube segment + # list[startCap, endCap][x, d1, d2, d3][nThroughWall][nAround] + self._boxExtNodeIds = None + self._rimExtNodeIds = None + + self._boxCoordinates = None + # list[startCap, endCap][[x, d1, d2, d3][nAcrossMajor][nAcrossMinor] + self._shellCoordinates = None + # list[startCap, endCap][[x, d1, d2, d3][nThroughWall][nAcrossMajor][nAcrossMinor]. + self._startCapNodeIds = None + # capNodeIds that form the cap at the start of a tube segment. + # list[nThroughWall][nAcrossMajor][nAcrossMinor]. + self._endCapNodeIds = None + # capNodeIds that form the cap at the end of a tube segment. + # list structure is identical to startCapNodeIds. + self._startCapElementIds = None + # elementIds that form the cap at the start of a tube segment. + # list[box, rim]: [box] and [rim] sublists have different structures. + # [box][core, transition, shield][nAcrossMajor][nAcrossMinor] + # [rim][base, shield][nAround] + self._endCapElementIds = None + # elementIds that form the cap at the end of a tube segment. + self._startExtElementIds = None + self._endExtElementIds = None + + # annotation groups created if core: + self._coreGroup = None + self._shellGroup = None + + def _extendTubeEnds(self): + """ + Add additional tube sections with smaller element size along the tube at either ends of the tube with smaller + D2 derivatives. This function is to minimise the effect of large difference in D2 derivatives between the cap + mesh and the tube mesh. + """ + isStartCap = self._isStartCap + idx = 0 if isStartCap else -1 + + layoutD1 = self._networkPathParameters[0][1][idx] + ext = self._getExtensionLength() + unitVector = normalize(layoutD1) + signValue = -1 if isStartCap else 1 + + coreBoxMajorNodesCount = self._elementsCountCoreBoxMajor + 1 + coreBoxMinorNodesCount = self._elementsCountCoreBoxMinor + 1 + + boxCoordinates, transitionCoordinates, shellCoordinates = [], [], [] + for nx in range(4): + boxCoordinates.append([]) + transitionCoordinates.append([]) + shellCoordinates.append([]) + for m in range(coreBoxMajorNodesCount): + boxCoordinates[nx].append([]) + for n in range(coreBoxMinorNodesCount): + boxCoordinates[nx][m].append([]) + if self._elementsCountTransition > 1: + for n3 in range(self._elementsCountTransition - 1): + transitionCoordinates[nx].append([]) + for n3 in range(self._elementsCountThroughShell + 1): + shellCoordinates[nx].append([]) + + for m in range(coreBoxMajorNodesCount): + xList, d2List = [], [] + x = self._tubeBoxCoordinates[0][idx][m] + for n in range(coreBoxMinorNodesCount): + tx = add(x[n], set_magnitude(unitVector, ext * signValue)) + td2 = mult(sub(tx, x[n]), signValue) + xList.append(tx) + d2List.append(td2) + boxCoordinates[0][m] = xList + boxCoordinates[1][m] = self._tubeBoxCoordinates[1][idx][m] + boxCoordinates[2][m] = d2List + boxCoordinates[3][m] = self._tubeBoxCoordinates[3][idx][m] + + if self._elementsCountTransition > 1: + for n3 in range(self._elementsCountTransition - 1): + xList, d2List = [], [] + x = self._tubeTransitionCoordinates[0][idx][n3] + for nx in range(self._elementsCountAround): + tx = add(x[nx], set_magnitude(unitVector, ext * signValue)) + td2 = mult(sub(tx, x[nx]), signValue) + xList.append(tx) + d2List.append(td2) + transitionCoordinates[0][n3] = xList + transitionCoordinates[1][n3] = self._tubeTransitionCoordinates[1][idx][n3] + transitionCoordinates[2][n3] = d2List + transitionCoordinates[3][n3] = self._tubeTransitionCoordinates[3][idx][n3] + + for n3 in range(self._elementsCountThroughShell + 1): + xList, d2List = [], [] + x = self._tubeShellCoordinates[0][idx][n3] + for nx in range(self._elementsCountAround): + tx = add(x[nx], set_magnitude(unitVector, ext * signValue)) + td2 = mult(sub(tx, x[nx]), signValue) + xList.append(tx) + d2List.append(td2) + shellCoordinates[0][n3] = xList + shellCoordinates[1][n3] = self._tubeShellCoordinates[1][idx][n3] + shellCoordinates[2][n3] = d2List + shellCoordinates[3][n3] = self._tubeShellCoordinates[3][idx][n3] + + self._boxExtCoordinates = [None, None] if self._boxExtCoordinates is None else self._boxExtCoordinates + self._boxExtCoordinates[idx] = boxCoordinates + if self._elementsCountTransition > 1: + self._transitionExtCoordinates = [None, None] if self._transitionExtCoordinates is None \ + else self._transitionExtCoordinates + self._transitionExtCoordinates[idx] = transitionCoordinates + + self._shellExtCoordinates = [None, None] if self._shellExtCoordinates is None else self._shellExtCoordinates + self._shellExtCoordinates[idx] = shellCoordinates + + def _remapCapCoordinates(self): + """ + Remap box and rim coordinates of the cap nodes based on the scale of tube extension. + """ + isStartCap = self._isStartCap + idx = 0 if isStartCap else -1 + + layoutD1 = self._networkPathParameters[0][1][idx] + ext = self._getExtensionLength() + unitVector = normalize(layoutD1) + signValue = -1 if isStartCap else 1 + + coreBoxMajorNodesCount = self._getNodesCountCoreBoxMajor() + coreBoxMinorNodesCount = self._getNodesCountCoreBoxMinor() + nodesCountRim = self._getNodesCountRim() + + for m in range(coreBoxMajorNodesCount): + xList = [] + x = self._boxCoordinates[idx][0][m] + for n in range(coreBoxMinorNodesCount): + tx = add(x[n], set_magnitude(unitVector, ext * signValue)) + xList.append(tx) + self._boxCoordinates[idx][0][m] = xList + for n3 in range(nodesCountRim): + for m in range(coreBoxMajorNodesCount): + xList = [] + x = self._shellCoordinates[idx][0][n3][m] + for n in range(coreBoxMinorNodesCount): + tx = add(x[n], set_magnitude(unitVector, ext * signValue)) + xList.append(tx) + self._shellCoordinates[idx][0][n3][m] = xList + + def _sampleCapCoordinates(self, s): + """ + Blackbox function for calculating coordinates and derivatives for the cap elements. + It first calculates the coordinates for shell nodes, then calculates for box nodes. + nodes, and then calculates the coordinates for rim nodes on the shell surface. + :param s: Index for isCap list. 0 indicates start cap and 1 indicates end cap. + """ + + def _getRatioBetweenTwoRadii(rList, oList): + """ + Calculates the ratio between the original radius of a sphere and the new radius of an ellipsoid. + :param rList: List of new radius in each direction. + :param oList: List of original radius in each direction. + :return: List of ratio between two radii in x, y, and z-direction. + """ + return [rList[c] / oList[c] for c in range(3)] + + self._isStartCap = isStartCap = True if self._isCap[0] and s == 0 else False + idx = 0 if isStartCap else -1 + centre = self._networkPathParameters[0][0][idx] + + self._extendTubeEnds() # extend tube end + + # shell nodes + nodesCountRim = self._getNodesCountRim() + for n3 in range(nodesCountRim): + ox = self._getRimExtCoordinatesAround(n3)[0] + radius = self._getRadius(ox) + radii = self._getTubeRadii(centre, n3, idx) # radii for spheroid + oRadii = [1.0, radius, radius] # original radii used to create the sphere + ratio = _getRatioBetweenTwoRadii(radii, oRadii) + # ratio between original radii for the sphere and the new radii for spheroid + self._calculateMajorAndMinorNodesCoordinates(n3, centre, ratio) + self._calculateShellQuadruplePoints(n3, centre, radius) + self._calculateShellRegularNodeCoordinates(n3, centre) + self._sphereToSpheroid(n3, ratio, centre) + self._determineShellDerivatives() + + # box nodes + self._calculateBoxQuadruplePoints(centre) + self._calculateBoxMajorAndMinorNodes() + self._determineBoxDerivatives() + + self._remapCapCoordinates() + self._smoothDerivatives(ratio) + + def _createShellCoordinatesList(self): + """ + Creates an empty list for storing rim coordinates. Only applies when the solid core is active. + """ + self._shellCoordinates = [] if self._shellCoordinates is None else self._shellCoordinates + elementsCountRim = self._getElementsCountRim() + for s in range(2): + self._shellCoordinates.append([] if self._isCap[s] else None) + if self._shellCoordinates[s] is not None: + for nx in range(4): + self._shellCoordinates[s].append([]) + for n3 in range(elementsCountRim): + self._shellCoordinates[s][nx].append([]) + self._shellCoordinates[s][nx][n3] = [[] for _ in range(self._elementsCountCoreBoxMajor + 1)] + for m in range(self._elementsCountCoreBoxMajor + 1): + self._shellCoordinates[s][nx][n3][m] = \ + [None for _ in range(self._elementsCountCoreBoxMinor + 1)] + + def _getOuterShellRadius(self): + """ + Calculates the radius of an outer shell. It takes the average of a half-distance between two opposing nodes on + the outer shell of a tube segment. + :return: Radius of the cap shell. + """ + ox = self._tubeShellCoordinates[0][0][-1] if self._isStartCap else self._tubeShellCoordinates[0][-1][-1] + radii = [magnitude(sub(ox[i], ox[i + self._elementsCountAround // 2])) / 2 for i in + range(self._elementsCountAround // 2)] + return sum(radii) / len(radii) + + def _getRadius(self, ox): + """ + Calculates the radius of a shell. It takes the average of a half-distance between two opposing nodes around a + tube segment. + :param ox: Coordinates of shell nodes around a tube segment. + :return: Radius of the cap shell. + """ + radii = [magnitude(sub(ox[i], ox[i + self._elementsCountAround // 2])) / 2 for i in + range(self._elementsCountAround // 2)] + return sum(radii) / len(radii) + + def _getExtensionLength(self): + """ + Calculates the length of extended tube segment. Currently set to half of the outer tube radius. + :return: Length of extended tube segment. + """ + outerRadius = self._getOuterShellRadius() + return outerRadius / 2 + + def _calculateMajorAndMinorNodesCoordinates(self, n3, centre, ratio): + """ + Calculates coordinates and derivatives for major and minor axis nodes on the surface of a cap shell by rotating + the major and minor axis nodes on the rim of a tube segment. + :param n3: Node index from inner to outer rim. + :param centre: Centre coordinates of a tube segment at either ends. + :param ratio: List of ratios between original circular radii and new radii if the tube is non-circular. + [x-axis, major axis, minor axis]. The values should equal 1.0 if the tube cross-section is circular. + """ + idx = 0 if self._isStartCap else -1 + layoutD2 = self._networkPathParameters[0][2][idx] + layoutD3 = self._networkPathParameters[0][4][idx] + + elementsCountAcrossMajor = self._elementsCountCoreBoxMajor + 2 + elementsCountAcrossMinor = self._elementsCountCoreBoxMinor + 2 + + refAxis = normalize(layoutD2) + rotateAngle = (math.pi / elementsCountAcrossMinor) if self._isStartCap else \ + -(math.pi / elementsCountAcrossMinor) + minorAxisNodesCoordinates = [[], []] # [startCap, endCap] + n1 = self._elementsCountAround * 3 // 4 + ix = self._getTubeRimCoordinates(n1, idx, n3) + for n in range(1, elementsCountAcrossMinor): + for nx in [0, 1]: # x coord and d1 derivatives + vi = sub(ix[nx], centre) if nx == 0 else mult(ix[nx], -1) + vi = div(vi, ratio[2]) + vr = rotate_vector_around_vector(vi, refAxis, n * rotateAngle) + vr = add(vr, centre) if nx == 0 else vr + minorAxisNodesCoordinates[nx].append(vr) + + refAxis = normalize(layoutD3) + rotateAngle = (math.pi / elementsCountAcrossMajor) if self._isStartCap else \ + -(math.pi / elementsCountAcrossMajor) + majorAxisNodesCoordinates = [[], []] # [startCap, endCap] + ix = self._getTubeRimCoordinates(0, idx, n3) + for m in range(1, elementsCountAcrossMajor): + for nx in [0, 1]: # [x, d1] + vi = sub(ix[nx], centre) if nx == 0 else ix[nx] + vi = div(vi, ratio[1]) + vr = rotate_vector_around_vector(vi, refAxis, m * rotateAngle) + vr = add(vr, centre) if nx == 0 else vr + majorAxisNodesCoordinates[nx].append(vr) + + midMajorIndex = elementsCountAcrossMajor // 2 - 1 + midMinorIndex = elementsCountAcrossMinor // 2 - 1 + for n in range(elementsCountAcrossMinor - 1): + for i in [0, 1]: + nx = [0, 1][i] + if self._shellCoordinates[idx][nx][n3][midMajorIndex][n] is None: + self._shellCoordinates[idx][nx][n3][midMajorIndex][n] = minorAxisNodesCoordinates[i][n] + for m in range(elementsCountAcrossMajor - 1): + for i in [0, 1]: + nx = [0, 2][i] + if self._shellCoordinates[idx][nx][n3][m][midMinorIndex] is None: + self._shellCoordinates[idx][nx][n3][m][midMinorIndex] = majorAxisNodesCoordinates[i][m] + + # derivatives + for n in range(elementsCountAcrossMinor - 1): + if self._shellCoordinates[idx][2][n3][midMajorIndex][n] is None: + self._shellCoordinates[idx][2][n3][midMajorIndex][n] = majorAxisNodesCoordinates[1][midMajorIndex] + tx = self._shellCoordinates[idx][0][n3][midMajorIndex] + td2 = self._shellCoordinates[idx][2][n3][midMajorIndex] + self._shellCoordinates[idx][2][n3][midMajorIndex] = smoothCubicHermiteDerivativesLine(tx, td2) + + for m in range(elementsCountAcrossMajor - 1): + if self._shellCoordinates[idx][1][n3][m][midMinorIndex] is None: + self._shellCoordinates[idx][1][n3][m][midMinorIndex] = minorAxisNodesCoordinates[1][midMinorIndex] + tx = [self._shellCoordinates[idx][0][n3][m][midMinorIndex] for m in range(elementsCountAcrossMajor - 1)] + td1 = [self._shellCoordinates[idx][1][n3][m][midMinorIndex] for m in range(elementsCountAcrossMajor - 1)] + sd1 = smoothCubicHermiteDerivativesLine(tx, td1) + for m in range(elementsCountAcrossMajor - 1): + self._shellCoordinates[idx][1][n3][m][midMinorIndex] = sd1[m] + + def _calculateShellQuadruplePoints(self, n3, centre, radius): + """ + Calculate coordinates and derivatives of the quadruple point on the surface, where 3 hex elements merge. + :param n3: Node index from inner to outer rim. + :param centre: Centre coordinates of a tube segment at either ends. + :param radius: Shell radius. + """ + idx = 0 if self._isStartCap else -1 + + layoutD1 = self._networkPathParameters[0][1][idx] + layoutD2 = self._networkPathParameters[0][2][idx] + layoutD3 = self._networkPathParameters[0][4][idx] + + axesList = [[mult(layoutD1, -1), layoutD2, mult(layoutD3, -1)], + [layoutD2, mult(layoutD1, -1), layoutD3], + [mult(layoutD2, -1), mult(layoutD1, -1), mult(layoutD3, -1)], + [mult(layoutD1, -1), mult(layoutD2, -1), layoutD3]] + + elementsCountUp = 2 + elementsCountAcrossMajor = self._elementsCountCoreBoxMajor + 2 + elementsCountAcrossMinor = self._elementsCountCoreBoxMinor + 2 + signValue = 1 if self._isStartCap else -1 + for counter, (m, n) in enumerate([(0, 0), (0, -1), (-1, 0), (-1, -1)]): + elementsCount = ( + [elementsCountUp, elementsCountAcrossMajor // 2, elementsCountAcrossMinor // 2] if m == n else + [elementsCountAcrossMajor // 2, elementsCountUp, elementsCountAcrossMinor // 2] + ) + + radiansPerElementAroundEllipse12 = math.pi / (2 * (elementsCount[0] + elementsCount[1] - 2)) + radiansPerElementAroundEllipse13 = math.pi / (2 * (elementsCount[0] + elementsCount[2] - 2)) + + theta_2 = (elementsCount[2] - 1) * radiansPerElementAroundEllipse13 + theta_3 = (elementsCount[1] - 1) * radiansPerElementAroundEllipse12 + phi_3 = calculate_azimuth(theta_3, theta_2) + + local_x = spherical_to_cartesian(radius, theta_3, phi_3) + c = counter if self._isStartCap else -(counter + 1) + axes = [mult(axis, signValue) for axis in axesList[c]] + x = local_to_global_coordinates(local_x, axes, centre) + + self._shellCoordinates[idx][0][n3][m][n] = x + + def _calculateShellRegularNodeCoordinates(self, n3, centre): + """ + Calculate coordinates and derivatives of all other shell nodes on the cap surface. + :param n3: Node index from inner to outer rim. + :param centre: Centre coordinates of a tube segment at either ends. + """ + idx = 0 if self._isStartCap else -1 + elementsCountAcrossMajor = self._elementsCountCoreBoxMajor + 2 + elementsCountAcrossMinor = self._elementsCountCoreBoxMinor + 2 + midMajorIndex = elementsCountAcrossMajor // 2 - 1 + midMinorIndex = elementsCountAcrossMinor // 2 - 1 + + elementsOut = self._elementsCountCoreBoxMinor // 2 + for m in range(self._elementsCountCoreBoxMajor + 1): + for n in [0, -1]: + x1 = self._shellCoordinates[idx][0][n3][m][n] + x2 = self._shellCoordinates[idx][0][n3][m][midMinorIndex] + if x1 is None: + continue + nx, nd1 = (sampleCurvesOnSphere(x1, x2, centre, elementsOut) if n == 0 else + sampleCurvesOnSphere(x2, x1, centre, elementsOut)) + start, end = ( + (n + 1, midMinorIndex) if n == 0 else (midMinorIndex + 1, self._elementsCountCoreBoxMinor)) + for c in range(start, end): + idx_c = c % (len(nx) - 1) + self._shellCoordinates[idx][0][n3][m][c] = nx[idx_c] + self._shellCoordinates[idx][1][n3][m][c] = nd1[idx_c] + self._shellCoordinates[idx][2][n3][m][c] = [0, 0, 0] + + elementsOut = self._elementsCountCoreBoxMajor // 2 + for n in range(self._elementsCountCoreBoxMinor + 1): + for m in [0, -1]: + x1 = self._shellCoordinates[idx][0][n3][m][n] + x2 = self._shellCoordinates[idx][0][n3][midMajorIndex][n] + if x1 is None: + continue + nx, nd2 = (sampleCurvesOnSphere(x1, x2, centre, elementsOut) if m == 0 else + sampleCurvesOnSphere(x2, x1, centre, elementsOut)) + start, end = ( + (m + 1, midMajorIndex) if m == 0 else (midMajorIndex + 1, self._elementsCountCoreBoxMajor)) + for c in range(start, end): + idx_c = c % (len(nx) - 1) + self._shellCoordinates[idx][0][n3][c][n] = nx[idx_c] + self._shellCoordinates[idx][1][n3][c][n] = [0, 0, 0] + self._shellCoordinates[idx][2][n3][c][n] = nd2[idx_c] + + def _sphereToSpheroid(self, n3, ratio, centre): + """ + Transform the sphere to ellipsoid using the radius in each direction. + :param n3: Node index from inner to outer rim. + :param ratio: List of ratios between original circular radii and new radii if the tube is non-circular. + [x-axis, major axis, minor axis]. The values should equal 1.0 if the tube cross-section is circular. + :param centre: Centre coordinates of a tube segment at either ends. + """ + idx = 0 if self._isStartCap else -1 + + # rotation angles to use in case when the cap is tilted from xyz axes. + layoutD2 = normalize(self._networkPathParameters[0][2][idx]) + layoutD3 = normalize(self._networkPathParameters[0][4][idx]) + thetaD2 = angle(layoutD2, [0.0, 1.0, 0.0]) + thetaD3 = angle(layoutD3, [0.0, 0.0, 1.0]) + + if layoutD2[0] < 0.0: + thetaD2 *= -1.0 + + mCount = self._elementsCountCoreBoxMajor + 1 + nCount = self._elementsCountCoreBoxMinor + 1 + # process shell coordinates + for m in range(mCount): + for n in range(nCount): + btx = self._shellCoordinates[idx][0][n3][m][n] + btx = sub(btx, centre) + # apply forward transformations + for vec, theta in [(layoutD3, thetaD2), (layoutD2, thetaD3)]: + btx = rotate_vector_around_vector(btx, vec, theta) + # scale by ratios + btx = [ratio[c] * btx[c] for c in range(3)] + # apply inverse transformations + for vec, theta in [(layoutD2, -thetaD3), (layoutD3, -thetaD2)]: + btx = rotate_vector_around_vector(btx, vec, theta) + self._shellCoordinates[idx][0][n3][m][n] = add(btx, centre) + + def _getTubeRadii(self, centre, n3, idx): + """ + Calculates the radius of a tube segment in major axis and minor axis. + :param centre: + :param n3: Node index from inner to outer rim. + :param idx: 0 if calculating for the start segment, -1 if calculating for the end segment. + :return: List of radii in major and minor axes. The radius in x-direction is set to 1.0 by default because the + radius in this direction is constant. + """ + n1m, n1n = 0, self._elementsCountAround // 4 + ixm, ixn = (self._getTubeRimCoordinates(n, idx, n3)[0] for n in [n1m, n1n]) + majorRadius, minorRadius = (magnitude(sub(coord, centre)) for coord in [ixm, ixn]) + xRadius = 1.0 + if majorRadius > minorRadius: + xRadius = math.pow((majorRadius / minorRadius), 1 / 3) + elif majorRadius < minorRadius: + xRadius = math.pow((minorRadius / majorRadius), 1 / 3) + return [xRadius, majorRadius, minorRadius] + + def _determineShellDerivatives(self): + """ + Compute d1, d2, and d3 derivatives for the shell nodes. + """ + idx = 0 if self._isStartCap else -1 + nodesCountRim = self._getNodesCountRim() + for n3 in range(nodesCountRim): + for m in [0, -1]: + for n in [0, -1]: + # initialise derivatives for quadruple points + mp = m + 1 if m == 0 else m - 1 + np = n + 1 if n == 0 else n - 1 + d1 = sub(self._shellCoordinates[idx][0][n3][mp][n], self._shellCoordinates[idx][0][n3][m][n]) + self._shellCoordinates[idx][1][n3][m][n] = mult(d1, -1) if m == -1 else d1 + d2 = sub(self._shellCoordinates[idx][0][n3][m][np], self._shellCoordinates[idx][0][n3][m][n]) + self._shellCoordinates[idx][2][n3][m][n] = mult(d2, -1) if n == -1 else d2 + + for n3 in range(nodesCountRim): + for m in range(self._elementsCountCoreBoxMajor + 1): + signValue = 1 if self._isStartCap else -1 + tx = self._shellCoordinates[idx][0][n3][m] + td2 = self._shellCoordinates[idx][2][n3][m] + sd2 = smoothCubicHermiteDerivativesLine(tx, td2) + for n in range(self._elementsCountCoreBoxMinor + 1): + self._shellCoordinates[idx][2][n3][m][n] = mult(sd2[n], signValue) + for n in range(self._elementsCountCoreBoxMinor + 1): + tx = [self._shellCoordinates[idx][0][n3][m][n] for m in range(self._elementsCountCoreBoxMajor + 1)] + td1 = [self._shellCoordinates[idx][1][n3][m][n] for m in range(self._elementsCountCoreBoxMajor + 1)] + sd1 = smoothCubicHermiteDerivativesLine(tx, td1) + for m in range(self._elementsCountCoreBoxMajor + 1): + self._shellCoordinates[idx][1][n3][m][n] = sd1[m] + + elementsCountRim = self._elementsCountThroughShell + self._elementsCountTransition - 1 + for m in range(self._elementsCountCoreBoxMajor + 1): + for n in range(self._elementsCountCoreBoxMinor + 1): + otx = self._shellCoordinates[idx][0][-1][m][n] + itx = self._shellCoordinates[idx][0][0][m][n] + shellFactor = 1.0 / elementsCountRim + sd3 = mult(sub(otx, itx), shellFactor) + for n3 in range(nodesCountRim): + self._shellCoordinates[idx][3][n3][m][n] = sd3 + + def _calculateBoxQuadruplePoints(self, centre): + """ + Calculate coordinates and derivatives of the quadruple point for the box elements, where 3 hex elements merge. + :param centre: Centre coordinates of a tube segment at either ends. + """ + idx = 0 if self._isStartCap else -1 + capBoxCoordinates = [] + for nx in range(4): + capBoxCoordinates.append([]) + capBoxCoordinates[nx] = [[] for _ in range(self._elementsCountCoreBoxMajor + 1)] + for m in range(self._elementsCountCoreBoxMajor + 1): + capBoxCoordinates[nx][m] = [None for _ in range(self._elementsCountCoreBoxMinor + 1)] + + boxCoordinates = self._tubeBoxCoordinates[0][idx] + rimCoordinates = self._tubeTransitionCoordinates[0][idx][0] if self._elementsCountTransition > 1 \ + else self._tubeShellCoordinates[0][idx][0] + capCoordinates = self._shellCoordinates[idx][0][0] + + elementsCountAround = self._elementsCountAround + nodesCountAcrossMinorHalf = len(boxCoordinates[0]) // 2 + triplePointIndexesList = [] + for n in range(0, elementsCountAround, elementsCountAround // 2): + triplePointIndexesList.append((n - nodesCountAcrossMinorHalf) % elementsCountAround) + triplePointIndexesList.append(n + nodesCountAcrossMinorHalf) + triplePointIndexesList[-2], triplePointIndexesList[-1] = triplePointIndexesList[-1], triplePointIndexesList[-2] + + for counter, (m, n) in enumerate([(0, 0), (0, -1), (-1, 0), (-1, -1)]): + tpIndex = triplePointIndexesList[counter] + x1, x2, x3 = boxCoordinates[m][n], rimCoordinates[tpIndex], capCoordinates[m][n] + ts = magnitude(sub(x1, x2)) + ra = sub(x3, centre) + radius = magnitude(ra) + local_x = mult(ra, (1 - ts / radius)) + capBoxCoordinates[0][m][n] = add(local_x, centre) + capBoxCoordinates[1][m][n] = self._tubeBoxCoordinates[1][idx][m][n] + capBoxCoordinates[3][m][n] = self._tubeBoxCoordinates[3][idx][m][n] + + if self._boxCoordinates is None: + self._boxCoordinates = [None] * 2 + self._boxCoordinates[idx] = capBoxCoordinates + + def _calculateBoxMajorAndMinorNodes(self): + """ + Calculate coordinates and derivatives for box nodes along the central major and minor axes. + """ + idx = 0 if self._isStartCap else -1 + midMajorIndex = self._elementsCountCoreBoxMajor // 2 + midMinorIndex = self._elementsCountCoreBoxMinor // 2 + + # box side nodes + for m in [0, -1]: + nx = [self._boxCoordinates[idx][0][m][n] for n in [0, -1]] + nd1 = [self._boxCoordinates[idx][1][m][n] for n in [0, -1]] + nd3 = [self._boxCoordinates[idx][3][m][n] for n in [0, -1]] + tx, td3, pe, pxi, psf = sampleCubicHermiteCurves(nx, nd3, self._elementsCountCoreBoxMinor, + arcLengthDerivatives=True) + td1 = interpolateSampleCubicHermite(nd1, [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + for n in range(1, self._elementsCountCoreBoxMinor): + self._boxCoordinates[idx][0][m][n] = tx[n] + self._boxCoordinates[idx][1][m][n] = td1[n] + self._boxCoordinates[idx][3][m][n] = td3[n] + for n in [0, -1]: + nx = [self._boxCoordinates[idx][0][m][n] for m in [0, -1]] + nd1 = [self._boxCoordinates[idx][1][m][n] for m in [0, -1]] + nd3 = [self._boxCoordinates[idx][3][m][n] for m in [0, -1]] + tx, td1, pe, pxi, psf = sampleCubicHermiteCurves(nx, nd1, self._elementsCountCoreBoxMajor, + arcLengthDerivatives=True) + td3 = interpolateSampleCubicHermite(nd3, [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + for m in range(1, self._elementsCountCoreBoxMajor): + self._boxCoordinates[idx][0][m][n] = tx[m] + self._boxCoordinates[idx][1][m][n] = td1[m] + self._boxCoordinates[idx][3][m][n] = td3[m] + + # box major and minor nodes + nx = [self._boxCoordinates[idx][0][midMajorIndex][n] for n in [0, -1]] + nd1 = [self._boxCoordinates[idx][1][midMajorIndex][n] for n in [0, -1]] + nd3 = [self._boxCoordinates[idx][3][midMajorIndex][n] for n in [0, -1]] + tx, td3, pe, pxi, psf = sampleCubicHermiteCurves(nx, nd3, self._elementsCountCoreBoxMinor, + arcLengthDerivatives=True) + td1 = interpolateSampleCubicHermite(nd1, [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + for n in range(1, self._elementsCountCoreBoxMinor): + self._boxCoordinates[idx][0][midMajorIndex][n] = tx[n] + self._boxCoordinates[idx][1][midMajorIndex][n] = td1[n] + self._boxCoordinates[idx][3][midMajorIndex][n] = td3[n] + + nx = [self._boxCoordinates[idx][0][m][midMinorIndex] for m in [0, -1]] + nd1 = [self._boxCoordinates[idx][1][m][midMinorIndex] for m in [0, -1]] + tx, td1, pe, pxi, psf = sampleCubicHermiteCurves(nx, nd1, self._elementsCountCoreBoxMajor, + arcLengthDerivatives=True) + td3 = interpolateSampleCubicHermite(nd3, [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + for m in range(1, self._elementsCountCoreBoxMajor): + self._boxCoordinates[idx][0][m][midMinorIndex] = tx[m] + self._boxCoordinates[idx][1][m][midMinorIndex] = td1[m] + self._boxCoordinates[idx][3][m][midMinorIndex] = td3[m] + + # remaining nodes + for m in range(self._elementsCountCoreBoxMajor): + for n in range(self._elementsCountCoreBoxMinor): + if self._boxCoordinates[idx][0][m][n] is None: + nx = [self._boxCoordinates[idx][0][i][n] for i in [0, midMajorIndex, -1]] + nd1 = [self._boxCoordinates[idx][1][i][n] for i in [0, midMajorIndex, -1]] + nd3 = [self._boxCoordinates[idx][3][i][n] for i in [0, midMajorIndex, -1]] + tx, td1, pe, pxi, psf = sampleCubicHermiteCurves(nx, nd1, self._elementsCountCoreBoxMajor, + arcLengthDerivatives=True) + td3 = interpolateSampleCubicHermite(nd3, [[0.0, 0.0, 0.0]] * 3, pe, pxi, psf)[0] + for mi in range(1, self._elementsCountCoreBoxMajor): + self._boxCoordinates[idx][0][mi][n] = tx[mi] + self._boxCoordinates[idx][1][mi][n] = td1[mi] + self._boxCoordinates[idx][3][mi][n] = td3[mi] + + # smooth derivatives + for m in range(self._elementsCountCoreBoxMajor + 1): + nx = self._boxCoordinates[idx][0][m] + nd3 = self._boxCoordinates[idx][3][m] + sd3 = smoothCubicHermiteDerivativesLine(nx, nd3) + for n in range(self._elementsCountCoreBoxMinor + 1): + self._boxCoordinates[idx][3][m][n] = sd3[n] + for n in range(self._elementsCountCoreBoxMinor + 1): + nx = [self._boxCoordinates[idx][0][m][n] for m in range(self._elementsCountCoreBoxMajor + 1)] + nd1 = [self._boxCoordinates[idx][1][m][n] for m in range(self._elementsCountCoreBoxMajor + 1)] + sd1 = smoothCubicHermiteDerivativesLine(nx, nd1) + for m in range(self._elementsCountCoreBoxMajor + 1): + self._boxCoordinates[idx][1][m][n] = sd1[m] + + def _determineBoxDerivatives(self): + """ + Calculate d2 derivatives of box nodes. + """ + idx = 0 if self._isStartCap else -1 + signValue = 1 if self._isStartCap else -1 + for m in range(self._elementsCountCoreBoxMajor + 1): + for n in range(self._elementsCountCoreBoxMinor + 1): + otx = self._tubeBoxCoordinates[0][idx][m][n] + itx = self._boxCoordinates[idx][0][m][n] + d2 = mult(sub(otx, itx), signValue) + self._boxCoordinates[idx][2][m][n] = d2 + + def _smoothDerivatives(self, ratio): + """ + Smooths derivatives to eliminate zero Jacobian contours at the cap-tube joint. + """ + nodesCountCoreBoxMajor = self._getNodesCountCoreBoxMajor() + nodesCountCoreBoxMinor = self._getNodesCountCoreBoxMinor() + nodesCountRim = self._getNodesCountRim() + nloop = self._getNodesCountRim() + + idx = 0 if self._isStartCap else -1 + if self._isCap[idx]: + boxCoordinates = self._boxCoordinates[idx] + boxExtCoordinates = self._boxExtCoordinates[idx] + shellCoordinates = self._shellCoordinates[idx] + tubeBoxCoordinates = self._tubeBoxCoordinates + boxBoundaryNodeToBoxId = self._createBoxBoundaryNodeIdsList() + midMajorIndex = self._elementsCountCoreBoxMajor // 2 + # smooth derivatives along d2 for box + for m in range(nodesCountCoreBoxMajor): + for n in range(nodesCountCoreBoxMinor): + nx = [boxCoordinates[0][m][n], boxExtCoordinates[0][m][n], tubeBoxCoordinates[0][idx][m][n]] + nd2 = [boxCoordinates[2][m][n], boxExtCoordinates[2][m][n], tubeBoxCoordinates[2][idx][m][n]] + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + boxCoordinates[2][m][n] = sd2[0] + boxExtCoordinates[2][m][n] = sd2[1] + + # smooth derivatives along d2 for shell + if ratio != [1.0, 1.0, 1.0]: + for n3 in range(nloop): + for m in range(nodesCountCoreBoxMajor): + cx = [shellCoordinates[0][n3][m][n] for n in range(nodesCountCoreBoxMinor)] + cd2 = [shellCoordinates[2][n3][m][n] for n in range(nodesCountCoreBoxMinor)] + + start = -self._elementsCountCoreBoxMinor // 2 - m + end = self._elementsCountCoreBoxMinor // 2 + m + rimIndex = [start, end] if idx == 0 else [end, start] + rimExt = self._getRimExtCoordinatesAround(n3) + ex = [rimExt[0][n2] for n2 in rimIndex] + ed2 = [rimExt[2][n2] for n2 in rimIndex] + + if idx == -1: + cx.reverse() + cd2.reverse() + + nx = [ex[0]] + cx + [ex[1]] + if idx == 0: + nd2 = [mult(ed2[0], -1.0)] + cd2 + [ed2[1]] + else: + nd2 = [ed2[0]] + cd2 + [mult(ed2[1], -1.0)] + + if m == midMajorIndex: + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)[ + 1:-1] + else: + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN)[ + 1:-1] + + if idx == -1: + sd2.reverse() + + for n in range(nodesCountCoreBoxMinor): + shellCoordinates[2][n3][m][n] = sd2[n] + + # smooth derivatives along d3 for shell + if self._isCore: + for m in range(nodesCountCoreBoxMajor): + for n in range(nodesCountCoreBoxMinor): + bx = boxCoordinates[0][m][n] + bd3 = sub(shellCoordinates[0][0][m][n], bx) + nx = [bx] + [shellCoordinates[0][n3][m][n] for n3 in range(nloop)] + nd3 = [bd3] + [shellCoordinates[3][n3][m][n] for n3 in range(nloop)] + sd3 = smoothCubicHermiteDerivativesLine(nx, nd3, fixAllDirections=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + for n3 in range(nloop): + shellCoordinates[3][n3][m][n] = sd3[n3 + 1] + + # smooth transition/shell ext coordinates along d2 + for n2 in range(self._elementsCountAround): + m, n = boxBoundaryNodeToBoxId[n2] + for n3 in range(nodesCountRim - 1): + tx = self._getRimExtCoordinatesAround(n3)[0][n2] + bx = shellCoordinates[0][n3][m][n] + bd2 = sub(bx, tx) + nx = [bx, tx] + nd2 = [bd2, self._getRimExtCoordinatesAround(n3)[2][n2]] + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2, fixAllDirections=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + if self._elementsCountTransition - 1 > n3: + self._transitionExtCoordinates[idx][2][n3][n2] = sd2[1] + else: + n3p = n3 - 1 - self._elementsCountTransition if self._elementsCountTransition > 1 else n3 + self._shellExtCoordinates[idx][2][n3p][n2] = sd2[1] + + # smooth shell ext coordinates along d3 + for n2 in range(self._elementsCountAround): + m, n = boxBoundaryNodeToBoxId[n2] + bx = boxCoordinates[0][m][n] + tx = [self._getRimExtCoordinatesAround(n3)[0][n2] for n3 in range(nodesCountRim)] + bd3 = sub(bx, tx[0]) + td3 = [self._getRimExtCoordinatesAround(n3)[3][n2] for n3 in range(nodesCountRim)] + nx = [bx] + tx + nd3 = [bd3] + td3 + sd3 = smoothCubicHermiteDerivativesLine(nx, nd3, fixAllDirections=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + self._shellExtCoordinates[idx][3][0][n2] = sd3[1] + + def _createBoxBoundaryNodeIdsList(self): + """ + 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. + :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. + """ + boxBoundaryNodeToBoxId = [] + elementsCountCoreBoxMajor = self._elementsCountCoreBoxMajor + elementsCountCoreBoxMinor = self._elementsCountCoreBoxMinor + coreBoxMajorNodesCount = elementsCountCoreBoxMajor + 1 + coreBoxMinorNodesCount = elementsCountCoreBoxMinor + 1 + + for n3 in range(coreBoxMajorNodesCount): + if n3 == 0 or n3 == coreBoxMajorNodesCount - 1: + n1List = list(range(coreBoxMinorNodesCount)) if n3 == 0 else ( + list(range(coreBoxMinorNodesCount - 1, -1, -1))) + for n1 in n1List: + boxBoundaryNodeToBoxId.append([n3, n1]) + else: + for n1 in [-1, 0]: + boxBoundaryNodeToBoxId.append([n3, n1]) + + start = elementsCountCoreBoxMajor - 2 + idx = elementsCountCoreBoxMinor + 2 + for n in range(int(start), -1, -1): + boxBoundaryNodeToBoxId.append(boxBoundaryNodeToBoxId.pop(idx + 2 * n)) + + nloop = elementsCountCoreBoxMinor // 2 + for _ in range(nloop): + boxBoundaryNodeToBoxId.insert(len(boxBoundaryNodeToBoxId), boxBoundaryNodeToBoxId.pop(0)) + + return boxBoundaryNodeToBoxId + + def _createBoundaryNodeIdsList(self, nodeIds): + """ + Creates a list (in a circular format similar to other rim node id lists) of box node ids that are + located at the boundary of the box component. + This list is used to easily stitch inner rim nodes with box nodes. + :param nodeIds: List of box node ids to be rearranged. + :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. + """ + capBoundaryNodeIds, capBoundaryNodeToBoxId = [], [] + boxElementsCountRow = self._elementsCountCoreBoxMajor + 1 + boxElementsCountColumn = self._elementsCountCoreBoxMinor + 1 + for n3 in range(boxElementsCountRow): + if n3 == 0 or n3 == boxElementsCountRow - 1: + ids = nodeIds[n3] if n3 == 0 else nodeIds[n3][::-1] + n1List = list(range(boxElementsCountColumn)) if n3 == 0 else ( + list(range(boxElementsCountColumn - 1, -1, -1))) + capBoundaryNodeIds += [ids[c] for c in range(boxElementsCountColumn)] + for n1 in n1List: + capBoundaryNodeToBoxId.append([n3, n1]) + else: + for n1 in [-1, 0]: + capBoundaryNodeIds.append(nodeIds[n3][n1]) + capBoundaryNodeToBoxId.append([n3, n1]) + + start = self._elementsCountCoreBoxMajor - 2 + idx = self._elementsCountCoreBoxMinor + 2 + for n in range(int(start), -1, -1): + capBoundaryNodeIds.append(capBoundaryNodeIds.pop(idx + 2 * n)) + capBoundaryNodeToBoxId.append(capBoundaryNodeToBoxId.pop(idx + 2 * n)) + + nloop = self._elementsCountCoreBoxMinor // 2 + for _ in range(nloop): + capBoundaryNodeIds.insert(len(capBoundaryNodeIds), capBoundaryNodeIds.pop(0)) + capBoundaryNodeToBoxId.insert(len(capBoundaryNodeToBoxId), capBoundaryNodeToBoxId.pop(0)) + + return capBoundaryNodeIds, capBoundaryNodeToBoxId + + def _getBoxCoordinates(self, m, n): + """ + :param m: Index along the major axis. + :param n: Index along the minor axis. + :return: x[], d1[], d2[], and d3[] of box nodes in the cap mesh. + """ + idx = 0 if self._isStartCap else -1 + return [self._boxCoordinates[idx][0][m][n], + self._boxCoordinates[idx][1][m][n], + self._boxCoordinates[idx][2][m][n], + self._boxCoordinates[idx][3][m][n]] + + def _getBoxExtCoordinates(self, m, n): + """ + :param m: Index along the major axis. + :param n: Index along the minor axis. + :return: x[], d1[], d2[], and d3[] of box nodes extended from the tube segment. + """ + idx = 0 if self._isStartCap else -1 + return [self._boxExtCoordinates[idx][0][m][n], + self._boxExtCoordinates[idx][1][m][n], + self._boxExtCoordinates[idx][2][m][n], + self._boxExtCoordinates[idx][3][m][n]] + + def _getRimExtCoordinates(self, n1, n3): + """ + :param n1: Index around rim. + :param n3: Index from inner to outer rim. + :return: x[], d1[], d2[], and d3[] of rim nodes extended from the tube segment. + """ + idx = 0 if self._isStartCap else -1 + transitionNodeCount = (len(self._transitionExtCoordinates[idx][0]) + if (self._transitionExtCoordinates and self._transitionExtCoordinates[idx]) else 0) + + if n3 < transitionNodeCount: + return [self._transitionExtCoordinates[idx][0][n3][n1], + self._transitionExtCoordinates[idx][1][n3][n1], + self._transitionExtCoordinates[idx][2][n3][n1], + self._transitionExtCoordinates[idx][3][n3][n1]] + sn3 = n3 - transitionNodeCount + return [self._shellExtCoordinates[idx][0][sn3][n1], + self._shellExtCoordinates[idx][1][sn3][n1], + self._shellExtCoordinates[idx][2][sn3][n1], + self._shellExtCoordinates[idx][3][sn3][n1]] + + def _getRimExtCoordinatesAround(self, n3): + """ + :param n3: Index from inner to outer rim. + :return: x[], d1[], d2[], and d3[] of rim nodes extended from the tube segment. + """ + idx = 0 if self._isStartCap else -1 + transitionNodeCount = (len(self._transitionExtCoordinates[idx][0]) + if (self._transitionExtCoordinates and self._transitionExtCoordinates[idx]) else 0) + + if n3 < transitionNodeCount: + return [self._transitionExtCoordinates[idx][0][n3], + self._transitionExtCoordinates[idx][1][n3], + self._transitionExtCoordinates[idx][2][n3], + self._transitionExtCoordinates[idx][3][n3]] + sn3 = n3 - transitionNodeCount + return [self._shellExtCoordinates[idx][0][sn3], + self._shellExtCoordinates[idx][1][sn3], + self._shellExtCoordinates[idx][2][sn3], + self._shellExtCoordinates[idx][3][sn3]] + + def _getRimCoordinatesWithCore(self, m, n, n2): + """ + Get coordinates and derivatives for cap rim. Only applies when core option is active. + :param m: Index across major axis. + :param n: Index across minor axis. + :param n2: Index along the tube. + :return: cap rim coordinates and derivatives for points at n3, m and n. + """ + idx = 0 if self._isStartCap else -1 + return [self._shellCoordinates[idx][0][n2][m][n], + self._shellCoordinates[idx][1][n2][m][n], + self._shellCoordinates[idx][2][n2][m][n], + self._shellCoordinates[idx][3][n2][m][n]] + + def _getTubeBoxCoordinates(self, m, n): + """ + Get coordinates and derivatives for tube box. + :param m: Index across major axis. + :param n: Index across minor axis. + :return: Tube box coordinates and derivatives for points at n2, m and n. + """ + idx = 0 if self._isStartCap else -1 + return [self._tubeBoxCoordinates[0][idx][m][n], + self._tubeBoxCoordinates[1][idx][m][n], + self._tubeBoxCoordinates[2][idx][m][n], + self._tubeBoxCoordinates[3][idx][m][n]] + + def _getTubeRimCoordinates(self, n1, n2, n3): + """ + Get rim parameters at a point. + :param n1: Node index around. + :param n2: Node index along segment. + :param n3: Node index from inner to outer rim. + :return: Tube rim coordinates and derivatives for points at n1, n2 and n3. + """ + transitionNodeCount = (len(self._tubeTransitionCoordinates[0][0]) + if (self._tubeTransitionCoordinates and self._tubeTransitionCoordinates[0]) else 0) + if n3 < transitionNodeCount and self._isCore: + return [self._tubeTransitionCoordinates[0][n2][n3][n1], + self._tubeTransitionCoordinates[1][n2][n3][n1], + self._tubeTransitionCoordinates[2][n2][n3][n1], + self._tubeTransitionCoordinates[3][n2][n3][n1]] + sn3 = n3 - transitionNodeCount + return [self._tubeShellCoordinates[0][n2][sn3][n1], + self._tubeShellCoordinates[1][n2][sn3][n1], + self._tubeShellCoordinates[2][n2][sn3][n1], + self._tubeShellCoordinates[3][n2][sn3][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._getNodesCountCoreBoxMinor() // 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, isShell=False): + """ + 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. + :param isShell: True if the triple point is located on the shell layer, False if located on the core box. + :return: Location identifier. + """ + em = self._elementsCountCoreBoxMinor // 2 + eM = self._elementsCountCoreBoxMajor // 2 + 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)) + + ni = len(lftColumnElements) // 2 + if e1 == topRowElements[0] or e1 == lftColumnElements[ni - 1]: + location = 1 # "TopLeft" + elif e1 == topRowElements[-1] or e1 == rhtColumnElements[0]: + location = -1 # "TopRight" + elif e1 == btmRowElements[-1] or e1 == lftColumnElements[ni]: + location = 2 # "BottomLeft" + elif e1 == btmRowElements[0] or e1 == rhtColumnElements[-1]: + location = -2 # "BottomRight" + else: + location = 0 + + if isShell and not self._isStartCap and location != 0: + location = location + 1 if location > 0 else location - 1 + if abs(location) > 2: + location = location - 2 if location > 0 else location + 2 + + return location + + def _getBoxBoundaryLocation(self, m, n): + """ + Determines: 1. Where the node is located on the box boundary. There are four side locations: Top, bottom, left, + and right; and 2. the location of a triple point relative to the solid core box. There are four locations: + Top left, top right, bottom left, and bottom right. + Location is 0 if not located at any of the four specified locations. + :return: Side location identifier and triple point location identifier. + """ + mEnd = self._getNodesCountCoreBoxMajor() - 1 + nEnd = self._getNodesCountCoreBoxMinor() - 1 + + m = m + self._getNodesCountCoreBoxMajor() if m < 0 else m + n = n + self._getNodesCountCoreBoxMinor() if n < 0 else n + + if 0 < m < mEnd: + location = 1 if n == nEnd else -1 if n == 0 else 0 + elif 0 < n < nEnd: + location = 2 if m == 0 else -2 if m == mEnd else 0 + else: + location = 0 + + tpLocation = 0 + if location == 0: + if n == nEnd: + tpLocation = 1 if m == 0 else -1 # Top Left or Top Right + elif n == 0: + tpLocation = 2 if m == 0 else -2 # Bottom Left or Bottom Right + + return location, tpLocation + + def _getNodesCountCoreBoxMajor(self): + idx = -1 if self._boxCoordinates[0] is None else 0 + return len(self._boxCoordinates[idx][0]) + + def _getNodesCountCoreBoxMinor(self): + idx = -1 if self._boxCoordinates[0] is None else 0 + return len(self._boxCoordinates[idx][0][0]) + + def _getNodesCountRim(self): + return self._elementsCountThroughShell + self._elementsCountTransition + + def _getElementsCountRim(self): + elementsCountRim = max(1, self._elementsCountThroughShell) + self._elementsCountTransition + return elementsCountRim + + def _generateNodes(self): + """ + Blackbox function for generating cap nodes. + """ + generateData = self._generateData + coordinates = generateData.getCoordinates() + fieldcache = generateData.getFieldcache() + nodes = generateData.getNodes() + nodetemplate = generateData.getNodetemplate() + + nodesCountCoreBoxMajor = self._getNodesCountCoreBoxMajor() + nodesCountCoreBoxMinor = self._getNodesCountCoreBoxMinor() + nloop = self._getNodesCountRim() + 1 + capNodeIds = [] + idx = 0 if self._isStartCap else -1 + for n3 in range(nloop): + if n3 == 0: + capNodeIds.append([]) + if not self._isCore: + continue + + for m in range(nodesCountCoreBoxMajor): + capNodeIds[n3].append([]) + for n in range(nodesCountCoreBoxMinor): + rx, rd1, rd2, rd3 = (self._boxCoordinates[idx][i][m][n] for i in range(4)) + 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, rd1, rd2, rd3]): + coordinates.setNodeParameters(fieldcache, -1, nodeValue, 1, rValue) + capNodeIds[n3][m].append(nodeIdentifier) + else: + capNodeIds.append([]) + for m in range(nodesCountCoreBoxMajor): + n3p = n3 - 1 + capNodeIds[n3].append([]) + for n in range(nodesCountCoreBoxMinor): + rx, rd1, rd2, rd3 = (self._shellCoordinates[idx][i][n3p][m][n] for i in range(4)) + 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, rd1, rd2, rd3]): + coordinates.setNodeParameters(fieldcache, -1, nodeValue, 1, rValue) + capNodeIds[n3][m].append(nodeIdentifier) + + if self._isStartCap: + self._startCapNodeIds = capNodeIds + else: + self._endCapNodeIds = capNodeIds + + def _generateExtendedTubeNodes(self): + """ + Blackbox function for generating tube nodes extended from the original tube segment. + """ + idx = 0 if self._isStartCap else -1 + generateData = self._generateData + coordinates = generateData.getCoordinates() + fieldcache = generateData.getFieldcache() + nodes = generateData.getNodes() + nodetemplate = generateData.getNodetemplate() + + # create core box nodes + self._boxExtNodeIds = [None, None] if self._boxExtNodeIds is None else self._boxExtNodeIds + if self._isCore: + self._boxExtNodeIds[idx] = [] + nodesCountCoreBoxMajor = self._getNodesCountCoreBoxMajor() + nodesCountAcrossMinor = self._getNodesCountCoreBoxMinor() + for n3 in range(nodesCountCoreBoxMajor): + self._boxExtNodeIds[idx].append([]) + rx, rd1, rd2, rd3 = [self._boxExtCoordinates[idx][i][n3] for i in range(4)] + 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._boxExtNodeIds[idx][n3].append(nodeIdentifier) + + # create rim nodes and transition nodes (if there are more than 1 layer of transition) + nodesCountRim = self._getNodesCountRim() + elementsCountTransition = self._elementsCountTransition + self._rimExtNodeIds = [None, None] if self._rimExtNodeIds is None else self._rimExtNodeIds + self._rimExtNodeIds[idx] = [] + for n3 in range(nodesCountRim): + n3p = n3 - (elementsCountTransition - 1) + tx = self._transitionExtCoordinates[idx] if self._isCore and elementsCountTransition > 1 and n3 < ( + elementsCountTransition - 1) else self._shellExtCoordinates[idx] + rx, rd1, rd2, rd3 = [ + tx[i][n3 if elementsCountTransition > 1 and n3 < (elementsCountTransition - 1) else n3p] for i in + range(4)] + + ringNodeIds = [] + for n1 in range(self._elementsCountAround): + 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) + ringNodeIds.append(nodeIdentifier) + self._rimExtNodeIds[idx].append(ringNodeIds) + + def _generateElements(self, annotationMeshGroups): + """ + Blackbox function for generating cap elements. Used only when the tube segment has a core. + """ + isStartCap = self._isStartCap + idx = 0 if isStartCap else -1 + elementsCountAround = self._elementsCountAround + elementsCountCoreBoxMinor = self._elementsCountCoreBoxMinor + elementsCountCoreBoxMajor = self._elementsCountCoreBoxMajor + elementsCountTransition = self._elementsCountTransition + elementsCountRim = self._getElementsCountRim() + + generateData = self._generateData + coordinates = generateData.getCoordinates() + mesh = generateData.getMesh() + elementtemplateStd, eftStd = generateData.getStandardElementtemplate() + + nodeLayoutTransition = generateData.getNodeLayoutTransition() + nodeLayoutCapTransition = generateData.getNodeLayoutCapTransition() + nodeLayoutCapTransitionSpecial = generateData.getNodeLayoutCapTransition(isSpecial=True) + + if isStartCap: + capNodeIds = self._startCapNodeIds + self._startCapElementIds = [] if self._startCapElementIds is None else self._startCapElementIds + else: + capNodeIds = self._endCapNodeIds + self._endCapElementIds = [] if self._endCapElementIds is None else self._endCapElementIds + boxExtNodeIds = self._boxExtNodeIds[idx] + rimExtNodeIds = self._rimExtNodeIds[idx] + triplePointIndexesList = self._getTriplePointIndexes() + + capElementIds = [] + boxBoundaryNodeIds, boxBoundaryNodeToCapIndex = [], [] + rimBoundaryNodeIds, rimBoundaryNodeToCapIndex = [], [] + for n2 in range(elementsCountRim + 1): + if n2 == 0: + if not self._isCore: + continue + + boxBoundaryNodeIds.append(self._createBoundaryNodeIdsList(capNodeIds[n2])[0]) + boxBoundaryNodeToCapIndex.append(self._createBoundaryNodeIdsList(capNodeIds[n2])[1]) + else: + rimBoundaryNodeIds.append(self._createBoundaryNodeIdsList(capNodeIds[n2])[0]) + rimBoundaryNodeToCapIndex.append(self._createBoundaryNodeIdsList(capNodeIds[n2])[1]) + + # box & shield + if self._isCore: + # box + boxElementIds = [] + for e3 in range(elementsCountCoreBoxMajor): + boxElementIds.append([]) + e3p = e3 + 1 + for e1 in range(elementsCountCoreBoxMinor): + nids, nodeParameters, nodeLayouts = [], [], [] + for n1 in [e1, e1 + 1]: + nids += [capNodeIds[0][e3][n1], capNodeIds[0][e3p][n1], + boxExtNodeIds[e3][n1], boxExtNodeIds[e3p][n1]] + if not isStartCap: + for a in [nids]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplateStd) + element.setNodesByIdentifier(eftStd, nids) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + boxElementIds[e3].append(elementIdentifier) + capElementIds.append(boxElementIds) + + # box shield elements (elements joining the box and the shield elements) + boxshieldElementIds = [] + for e3 in range(elementsCountCoreBoxMajor): + boxshieldElementIds.append([]) + e3p = e3 + 1 + for e1 in range(elementsCountCoreBoxMinor): + nids, nodeParameters, nodeLayouts = [], [], [] + elementIdentifier = generateData.nextElementIdentifier() + for n1 in [e1, e1 + 1]: + for n3 in [e3, e3p]: + nids += [capNodeIds[1][n3][n1]] + nodeParameter = self._getRimCoordinatesWithCore(n3, n1, 0) + nodeParameters.append(nodeParameter) + if elementsCountTransition > 1: + nodeLayouts.append(nodeLayoutCapTransitionSpecial) + else: + nodeLayouts.append(nodeLayoutCapTransition) + for n3 in [e3, e3p]: + boxLocation, tpLocation = self._getBoxBoundaryLocation(n3, n1) + nid = capNodeIds[0][n3][n1] + nids += [nid] + nodeParameter = self._getBoxCoordinates(n3, n1) + nodeParameters.append(nodeParameter) + nodeLayoutCapBoxShield = generateData.getNodeLayoutCapBoxShield(boxLocation, isStartCap) + nodeLayoutCapBoxShieldTriplePoint = generateData.getNodeLayoutCapBoxShieldTriplePoint( + tpLocation, isStartCap) + if nid in boxBoundaryNodeIds[0]: + nodeLayouts.append( + nodeLayoutCapBoxShield if tpLocation == 0 else nodeLayoutCapBoxShieldTriplePoint) + else: + nodeLayouts.append(None) + if not isStartCap: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + boxshieldElementIds[e3].append(elementIdentifier) + capElementIds.append(boxshieldElementIds) + + # shield + for e3 in range(elementsCountRim - 1): + e3p = e3 + 1 + shieldElementIds = [] + for e2 in range(self._elementsCountCoreBoxMajor): + e2p = e2 + 1 + shieldElementIds.append([]) + for e1 in range(elementsCountCoreBoxMinor): + e1p = e1 + 1 + nids = [] + for n3 in [e3p, e3p + 1]: + nids += [capNodeIds[n3][e2][e1], capNodeIds[n3][e2p][e1], + capNodeIds[n3][e2][e1p], capNodeIds[n3][e2p][e1p]] + if not isStartCap: + for a in [nids]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplateStd) + element.setNodesByIdentifier(eftStd, nids) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + shieldElementIds[e2].append(elementIdentifier) + capElementIds.append(shieldElementIds) + + if isStartCap: + self._startCapElementIds.append(capElementIds) + else: + self._endCapElementIds.append(capElementIds) + + # rim + capElementIds = [] + if self._isCore: + # box shell (elements between the box core and the outer rim) + ringElementIds = [] + boxExtBoundaryNodeIds, boxExtBoundaryNodestoBoxIds = self._createBoundaryNodeIdsList(boxExtNodeIds) + for e1 in range(elementsCountAround): + nids, nodeParameters, nodeLayouts = [], [], [] + n1p = (e1 + 1) % self._elementsCountAround + boxLocation = self._getTriplePointLocation(e1) + shellLocation = self._getTriplePointLocation(e1, isShell=True) + nodeLayoutTransitionTriplePoint = generateData.getNodeLayoutTransitionTriplePoint(boxLocation) + nodeLayoutCapShellTransitionTriplePoint = generateData.getNodeLayoutCapShellTriplePoint(shellLocation) + for n3 in [0, 1]: + for n1 in [e1, n1p]: + nid = boxBoundaryNodeIds[n3][n1] if n3 == 0 else rimBoundaryNodeIds[n3 - 1][n1] + nids += [nid] + mi, ni = boxBoundaryNodeToCapIndex[n3][n1] if n3 == 0 else rimBoundaryNodeToCapIndex[n3 - 1][n1] + location, tpLocation = self._getBoxBoundaryLocation(mi, ni) + nodeLayoutCapBoxShield = generateData.getNodeLayoutCapBoxShield(location, isStartCap) + nodeLayoutCapBoxShieldTriplePoint = generateData.getNodeLayoutCapBoxShieldTriplePoint( + tpLocation, isStartCap) + if n3 == 0: + nodeParameter = self._getBoxCoordinates(mi, ni) + nodeLayout = nodeLayoutCapBoxShieldTriplePoint if n1 in triplePointIndexesList else ( + nodeLayoutCapBoxShield) + else: + nodeParameter = self._getRimCoordinatesWithCore(mi, ni, 0) + nodeLayout = nodeLayoutCapShellTransitionTriplePoint if n1 in triplePointIndexesList else ( + nodeLayoutCapTransition) + nodeParameters.append(nodeParameter) + nodeLayouts.append(nodeLayout) + for n1 in [e1, n1p]: + if n3 == 0: + nid = boxExtBoundaryNodeIds[n1] + mi, ni = boxExtBoundaryNodestoBoxIds[n1] + nodeParameter = self._getBoxExtCoordinates(mi, ni) + else: + nid = rimExtNodeIds[0][n1] + nodeParameter = self._getRimExtCoordinates(n1, 0) + nids += [nid] + nodeParameters.append(nodeParameter) + nodeLayouts.append(nodeLayoutTransitionTriplePoint if n1 in triplePointIndexesList and n3 == 0 + else nodeLayoutTransition) + if not isStartCap: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + 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) + capElementIds.append(ringElementIds) + + # shell + triplePointIndexesList = self._getTriplePointIndexes() + for e3 in range(elementsCountRim - 1): + rimElementIds = [] + for e1 in range(self._elementsCountAround): + nids, nodeParameters, nodeLayouts = [], [], [] + e1p = (e1 + 1) % self._elementsCountAround + location = self._getTriplePointLocation(e1, isShell=True) + nodeLayoutCapTransition = generateData.getNodeLayoutCapTransition() + nodeLayoutCapShellTriplePoint = generateData.getNodeLayoutCapShellTriplePoint(location) + for n3 in [e3, e3 + 1]: + for n1 in [e1, e1p]: + nids += [rimBoundaryNodeIds[n3][n1]] + mi, ni = rimBoundaryNodeToCapIndex[n3][n1] + nodeParameter = self._getRimCoordinatesWithCore(mi, ni, n3) + nodeParameters.append(nodeParameter) + nodeLayouts.append(nodeLayoutCapShellTriplePoint if n1 in triplePointIndexesList else + nodeLayoutCapTransition) + for n1 in [e1, e1p]: + nids += [rimExtNodeIds[n3][n1]] + nodeParameter = self._getRimExtCoordinates(n1, n3) + nodeParameters.append(nodeParameter) + nodeLayouts.append(None) + if not isStartCap: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + elementIdentifier = generateData.nextElementIdentifier() + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + rimElementIds.append(elementIdentifier) + capElementIds.append(rimElementIds) + + if isStartCap: + self._startCapElementIds.append(capElementIds) + else: + self._endCapElementIds.append(capElementIds) + + def _generateExtendedTubeElements(self, tubeBoxNodeIds, tubeRimNodeIds, annotationMeshGroups): + """ + Blackbox function for generating extended tube elements. + :param tubeBoxNodeIds: List of tube box nodes. + :param tubeRimNodeIds: List of tube rim nodes. + """ + isStartCap = self._isStartCap + idx = 0 if isStartCap else -1 + scalingMode = 3 if isStartCap else 4 + generateData = self._generateData + coordinates = generateData.getCoordinates() + mesh = generateData.getMesh() + elementtemplateStd, eftStd = generateData.getStandardElementtemplate() + boxExtElementIds, rimExtElementIds = [], [] + + boxExtNodeIds = self._boxExtNodeIds[idx] + rimExtNodeIds = self._rimExtNodeIds[idx] + + if isStartCap: + self._startExtElementIds = [] if self._startExtElementIds is None else self._startExtElementIds + else: + self._endExtElementIds = [] if self._endExtElementIds is None else self._endExtElementIds + + if self._isCore: + boxExtBoundaryNodeIds, boxExtBoundaryNodesToBoxIds = self._createBoundaryNodeIdsList(boxExtNodeIds) + tubeBoxBoundaryNodeIds, tubeBoxBoundaryNodesToBoxIds = self._createBoundaryNodeIdsList(tubeBoxNodeIds[idx]) + # create box elements + boxElementIds = [] + for e3 in range(self._elementsCountCoreBoxMajor): + boxElementIds.append([]) + e3p = e3 + 1 + for e1 in range(self._elementsCountCoreBoxMinor): + nids = [] + for n1 in [e1, e1 + 1]: + nids += [boxExtNodeIds[e3][n1], boxExtNodeIds[e3p][n1], + tubeBoxNodeIds[idx][e3][n1], tubeBoxNodeIds[idx][e3p][n1]] + if not isStartCap: + for a in [nids]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + elementIdentifier = generateData.nextElementIdentifier() + element = mesh.createElement(elementIdentifier, elementtemplateStd) + element.setNodesByIdentifier(eftStd, nids) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + boxElementIds[e3].append(elementIdentifier) + boxExtElementIds.append(boxElementIds) + + # create core transition elements first layer after box + triplePointIndexesList = self._getTriplePointIndexes() + boxElementIds = [] + 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 [0, 1]: + for n1 in [e1, n1p]: + if n2 == 0: + nid = boxExtBoundaryNodeIds[n1] + mi, ni = boxExtBoundaryNodesToBoxIds[n1] + nodeParameter = self._getBoxExtCoordinates(mi, ni) + else: + nid = tubeBoxBoundaryNodeIds[n1] + mi, ni = tubeBoxBoundaryNodesToBoxIds[n1] + nodeParameter = self._getTubeBoxCoordinates(mi, ni) + nids += [nid] + nodeParameters.append(nodeParameter) + nodeLayouts.append(nodeLayoutTransitionTriplePoint if n1 in triplePointIndexesList else + nodeLayoutTransition) + if not isStartCap: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + for n2 in [0, 1]: + for n1 in [e1, n1p]: + if n2 == 0: + nid = rimExtNodeIds[0][n1] + nodeParameter = self._getRimExtCoordinates(n1, 0) + else: + nid = tubeRimNodeIds[idx][0][n1] + nodeParameter = self._getTubeRimCoordinates(n1, idx, 0) + nids += [nid] + nodeParameters.append(nodeParameter) + nodeLayouts.append(None) + if not isStartCap: + for a in [nids, nodeParameters, nodeLayouts]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft, scalefactors = determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts) + if self._elementsCountTransition == 1: + eft, scalefactors = generateData.resolveEftCoreBoundaryScaling( + eft, scalefactors, nodeParameters, nids, scalingMode) + 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) + boxElementIds.append(elementIdentifier) + boxExtElementIds.append(boxElementIds) + + # create regular rim elements - all elements outside first transition layer + elementsCountRim = self._getElementsCountRim() + elementsCountRimRegular = elementsCountRim - 1 + nTransition = elementsCountRimRegular - self._elementsCountThroughShell + for e3 in range(elementsCountRimRegular): + if e3 < nTransition: + boxElementIds = [] + else: + rimExtElementIds.append([]) + lastTransition = self._isCore and (e3 == (self._elementsCountTransition - 2)) + for e1 in range(self._elementsCountAround): + elementtemplate, eft = elementtemplateStd, eftStd + n1p = (e1 + 1) % self._elementsCountAround + nids = [] + for n3 in [e3, e3 + 1]: + nids += [rimExtNodeIds[n3][e1], rimExtNodeIds[n3][n1p], + tubeRimNodeIds[idx][n3][e1], tubeRimNodeIds[idx][n3][n1p]] + if not isStartCap: + for a in [nids]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + elementIdentifier = generateData.nextElementIdentifier() + scalefactors = [] + if lastTransition: + # get node parameters for computing scale factors + nodeParameters = [] + for n3 in (e3, e3 + 1): + for n2 in (0, 1): + for n1 in (e1, n1p): + if n2 == 0: + nodeParameter = self._getRimExtCoordinates(n1, n3) + else: + nodeParameter = self._getTubeRimCoordinates(n1, idx, n3) + nodeParameters.append(nodeParameter) + if not isStartCap: + for a in [nodeParameters]: + a[-4], a[-2] = a[-2], a[-4] + a[-3], a[-1] = a[-1], a[-3] + eft = generateData.createElementfieldtemplate() + eft, scalefactors = generateData.resolveEftCoreBoundaryScaling( + eft, scalefactors, nodeParameters, nids, scalingMode) + elementtemplateTransition = mesh.createElementtemplate() + elementtemplateTransition.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplateTransition.defineField(coordinates, -1, eft) + elementtemplate = elementtemplateTransition + element = mesh.createElement(elementIdentifier, elementtemplate) + element.setNodesByIdentifier(eft, nids) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) + if e3 < nTransition: + boxElementIds.append(elementIdentifier) + else: + rimExtElementIds[-1].append(elementIdentifier) + if e3 < nTransition: + boxExtElementIds.append(boxElementIds) + + if isStartCap: + self._startExtElementIds.append(boxExtElementIds) + self._startExtElementIds.append(rimExtElementIds) + else: + self._endExtElementIds.append(boxExtElementIds) + self._endExtElementIds.append(rimExtElementIds) + + def sampleCoordinates(self, tubeBoxCoordinates, tubeTransitionCoordinates, tubeShellCoordinates): + """ + Sample cap coordinates. + """ + self._tubeBoxCoordinates = tubeBoxCoordinates # tube box coordinates + self._tubeTransitionCoordinates = tubeTransitionCoordinates # tube transition coordinates + self._tubeShellCoordinates = tubeShellCoordinates # tube rim coordinates + + self._createShellCoordinatesList() + for s in range(2): + self._isStartCap = True if self._isCap[0] and s == 0 else False + if self._isCap[s]: + self._sampleCapCoordinates(s) + else: + continue + + return self._boxExtCoordinates, self._transitionExtCoordinates, self._shellExtCoordinates + + def generateNodes(self, generateData, isStartCap=True): + """ + Blackbox function for generating cap and extended tube nodes. + :param generateData: Class object from TubeNetworkMeshGenerateData. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + """ + self._isStartCap = isStartCap + self._generateData = generateData + if isStartCap: + self._generateNodes() + self._generateExtendedTubeNodes() + else: + self._generateExtendedTubeNodes() + self._generateNodes() + + def generateElements(self, tubeBoxNodeIds, tubeRimNodeIds, annotationMeshGroups, isStartCap=True): + """ + Blackbox function for generating cap and extended tube elements. + :param tubeBoxNodeIds: List of tube box nodes. + :param tubeRimNodeIds: List of tube rim nodes. + :param annotationMeshGroups: List of all annotated mesh groups. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + """ + self._isStartCap = isStartCap + if isStartCap: + self._generateElements(annotationMeshGroups) + self._generateExtendedTubeElements(tubeBoxNodeIds, tubeRimNodeIds, annotationMeshGroups) + else: + self._generateExtendedTubeElements(tubeBoxNodeIds, tubeRimNodeIds, annotationMeshGroups) + self._generateElements(annotationMeshGroups) + + def addCapBoxElementsToMeshGroup(self, meshGroup, e1Range=None, e2Range=None, e3Range=None, mode=0): + """ + Add ranges of box elements to mesh group. + :param meshGroup: Zinc MeshGroup to add elements to. + :param e1Range: Range between start and limit element indexes in major / d2 direction. + :param e2Range: Range between start and limit element indexes around the tube. + :param e3Range: Range between start and limit element indexes in minor / d3 direction. + :param mode: + """ + elementsCountAround = self._elementsCountAround + elementsCountCoreBoxMajor = self._elementsCountCoreBoxMajor + elementsCountCoreBoxMinor = self._elementsCountCoreBoxMinor + elementsCountTransition = self._elementsCountTransition + + mesh = meshGroup.getMasterMesh() + e1Range = range(elementsCountCoreBoxMajor) if e1Range is None else e1Range + e2Range = range(elementsCountAround) if e2Range is None else e2Range + e3Range = range(elementsCountCoreBoxMinor) if e3Range is None else e3Range + + nloop = 2 if mode == 0 else 1 + for idx in range(nloop): + if mode > 0: + idx = 0 if mode == 1 else 1 + isCap = self._isCap[idx] + if not isCap: + continue + + isStart = (idx == 0) + capElementIds = self._startCapElementIds if isStart else self._endCapElementIds + extElementIds = self._startExtElementIds if isStart else self._endExtElementIds + + # Cap base block (structured m, t, n) + addElementsFromIdentifiers(mesh, meshGroup, capElementIds[0], tRange=range(elementsCountTransition + 1), mRange=e1Range, nRange=e3Range) + + # Cap tube elements (t, e2) + addElementsFromIdentifiers(mesh, meshGroup, capElementIds[1], tRange=range(elementsCountTransition), e2Range=e2Range) + + # Extension base block (just one layer at t=0) + addElementsFromIdentifiers(mesh, meshGroup, extElementIds[0][0:1], tRange=[0], mRange=e1Range, nRange=e3Range) + + # Extension tube (t > 0, e2) + addElementsFromIdentifiers(mesh, meshGroup, extElementIds[0], tRange=range(1, elementsCountTransition + 1), e2Range=e2Range) + + def addCapShellElementsToMeshGroup(self, meshGroup, e1Range=None, e2Range=None, e3Range=None, mode=0): + """ + Add ranges of shell elements to mesh group. + :param meshGroup: Zinc MeshGroup to add elements to. + :param e1Range: Range between start and limit element indexes in major / d2 direction. + :param e2Range: Range between start and limit element indexes around the tube. + :param e3Range: Range between start and limit element indexes in minor / d3 direction. + :mode: + """ + elementsCountAround = self._elementsCountAround + elementsCountCoreBoxMajor = self._elementsCountCoreBoxMajor + elementsCountCoreBoxMinor = self._elementsCountCoreBoxMinor + elementsCountTransition = self._elementsCountTransition + elementsCountThroughShell = self._elementsCountThroughShell + + mesh = meshGroup.getMasterMesh() + e1Range = range(elementsCountCoreBoxMajor) if e1Range is None else e1Range + e2Range = range(elementsCountAround) if e2Range is None else e2Range + e3Range = range(elementsCountCoreBoxMinor) if e3Range is None else e3Range + + nloop = 2 if mode == 0 else 1 + for idx in range(nloop): + if mode > 0: + idx = 0 if mode == 1 else 1 + isCap = self._isCap[idx] + if not isCap: + continue + + isStart = (idx == 0) + capElementIds = self._startCapElementIds if isStart else self._endCapElementIds + extElementIds = self._startExtElementIds if isStart else self._endExtElementIds + + tCapShellRange = range(elementsCountTransition + 1, elementsCountTransition + 1 + elementsCountThroughShell) + addElementsFromIdentifiers(mesh, meshGroup, capElementIds[0], tRange=tCapShellRange, mRange=e1Range, nRange=e3Range) + + tTubeShellRange = range(elementsCountTransition, elementsCountTransition + elementsCountThroughShell) + addElementsFromIdentifiers(mesh, meshGroup, capElementIds[1], tRange=tTubeShellRange, e2Range=e2Range) + + tExtShellRange = range(elementsCountThroughShell) + addElementsFromIdentifiers(mesh, meshGroup, extElementIds[1], tRange=tExtShellRange, e2Range=e2Range) + + + def addCapSideD1ElementsToMeshGroup(self, side: bool, meshGroup): + """ + Add elements to the mesh group on side of +d1 or -d1, often matching anterior and posterior. + :param side: False for +d1 direction, True for -d1 direction. + :param meshGroup: Zinc MeshGroup to add elements to. + """ + mode = 1 if side else 2 + if self._isCore: + self.addCapBoxElementsToMeshGroup(meshGroup, mode=mode) + self.addCapShellElementsToMeshGroup(meshGroup, mode=mode) + + def addCapSideD2ElementsToMeshGroup(self, side: bool, meshGroup): + """ + Add elements to the mesh group on side of +d2 or -d2, often matching left and right. + Only works with even numbers around and phase starting at +d2. + :param side: False for +d2 direction, True for -d2 direction. + :param meshGroup: Zinc MeshGroup to add elements to. + """ + e2Start = (self._elementsCountAround // 4) if side else -((self._elementsCountAround + 2) // 4) + e2Limit = e2Start + (self._elementsCountAround // 2) + e2Limit = e2Limit + 1 if (self._elementsCountAround % 4) == 2 else e2Limit + e2Range = range(e2Start, e2Limit) + if self._isCore: + e1Start = (self._elementsCountCoreBoxMajor // 2) if side else 0 + e1Limit = self._elementsCountCoreBoxMajor if side else ((self._elementsCountCoreBoxMajor + 1) // 2) + e1Range = range(e1Start, e1Limit) + self.addCapBoxElementsToMeshGroup(meshGroup, e1Range=e1Range, e2Range=e2Range) + else: + e1Range = None + self.addCapShellElementsToMeshGroup(meshGroup, e1Range=e1Range, e2Range=e2Range) + + def addCapSideD3ElementsToMeshGroup(self, side: bool, meshGroup): + """ + Add elements to the mesh group on side of +d3 or -d3, often matching anterior/ventral and posterior/dorsal. + Only works with even numbers around and phase starting at +d2. + :param side: False for +d3 direction, True for -d3 direction. + :param meshGroup: Zinc MeshGroup to add elements to. + """ + e2Start = (self._elementsCountAround // 2) if side else 0 + e2Limit = e2Start + (self._elementsCountAround // 2) + e2Range = range(e2Start, e2Limit) + if self._isCore: + e3Start = 0 if side else (self._elementsCountCoreBoxMinor // 2) + e3Limit = ((self._elementsCountCoreBoxMinor + 1) // 2) if side else self._elementsCountCoreBoxMinor + e3Range = range(e3Start, e3Limit) + self.addCapBoxElementsToMeshGroup(meshGroup, e2Range=e2Range, e3Range=e3Range) + else: + e3Range = None + self.addCapShellElementsToMeshGroup(meshGroup, e2Range=e2Range, e3Range=e3Range) + + +def sampleCurvesOnSphere(x1, x2, origin, elementsOut): + """ + Sample coordinates and d1 derivatives of + :param x1: Coordinates of point 1 on the spherical surface of cap mesh. + :param x2: Coordinates of point 2 on the spherical surface of cap mesh. + :param origin: Centre point coordinates. + :param elementsOut: The number of elements required between points 1 and 2. + :return: Lists of sampled x and d1 between points 1 and 2. + """ + r1, r2 = sub(x1, origin), sub(x2, origin) + deltax = sub(r2, r1) + normal = cross(r1, deltax) + theta = angle(r1, r2) + anglePerElement = theta / elementsOut + arcLengthPerElement = calculate_arc_length(x1, x2, origin) / elementsOut + + nx, nd1 = [], [] + for n1 in range(elementsOut + 1): + radiansAcross = n1 * anglePerElement + r = rotate_vector_around_vector(r1, normal, radiansAcross) + nx.append(add(r, origin)) + nd1.append(set_magnitude(cross(normal, r), arcLengthPerElement)) + + return nx, nd1 + + +def addElementsFromIdentifiers(mesh, meshGroup, identifiers, tRange=None, mRange=None, nRange=None, e2Range=None): + """ + Adds elements to the mesh group based on a structured list of element identifiers. + :param mesh: The master mesh object used to look up elements by identifier. + :param meshGroup: Zinc MeshGroup to add elements to. + :param identifiers: A nested list of element identifiers. + :param tRange: Range of transition indices. + :param mRange: Range along the major axis. + :param nRange: Range along the minor axis. + :param e2Range: Range around the tube, for ring/circular indexing. + """ + for t in tRange: + if mRange is not None and nRange is not None: + for m in mRange: + for n in nRange: + elementIdentifier = identifiers[t][m][n] + element = mesh.findElementByIdentifier(elementIdentifier) + meshGroup.addElement(element) + elif e2Range is not None: + for e2 in e2Range: + elementIdentifier = identifiers[t][e2] + element = mesh.findElementByIdentifier(elementIdentifier) + meshGroup.addElement(element) diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index 2594eaac..6b8222f4 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -22,12 +22,13 @@ def getEftTermScaling(eft, functionIndex, termIndex): scaleFactorCount, scaleFactorIndexes = eft.getTermScaling(functionIndex, termIndex, 1) if scaleFactorCount < 0: if eft.getNumberOfLocalScaleFactors() > 0: - print('getEftTermScaling function', functionIndex, ' term', termIndex, ' scaleFactorCount ', scaleFactorCount) + print('getEftTermScaling function', functionIndex, ' term', termIndex, ' scaleFactorCount ', + scaleFactorCount) return [] if scaleFactorCount == 0: return [] if scaleFactorCount == 1: - return [ scaleFactorIndexes ] + return [scaleFactorIndexes] scaleFactorCount, scaleFactorIndexes = eft.getTermScaling(functionIndex, termIndex, scaleFactorCount) return scaleFactorIndexes @@ -41,7 +42,8 @@ def mapEftFunction1Node1Term(eft, function, localNode, valueLabel, version, scal eft.setTermScaling(function, 1, scaleFactors) -def mapEftFunction1Node2Terms(eft, function, localNode, valueLabel1, version1, scaleFactors1, valueLabel2, version2, scaleFactors2): +def mapEftFunction1Node2Terms(eft, function, localNode, valueLabel1, version1, scaleFactors1, valueLabel2, version2, + scaleFactors2): ''' Set function of eft to sum 2 terms the respective valueLabels, versions from localNode with scaleFactors ''' @@ -62,7 +64,8 @@ def remapEftLocalNodes(eft, newNodeCount, localNodeIndexes): for f in range(1, functionCount + 1): termCount = eft.getFunctionNumberOfTerms(f) for t in range(1, termCount + 1): - eft.setTermNodeParameter(f, t, localNodeIndexes[eft.getTermLocalNodeIndex(f, t) - 1], eft.getTermNodeValueLabel(f, t), eft.getTermNodeVersion(f, t)) + eft.setTermNodeParameter(f, t, localNodeIndexes[eft.getTermLocalNodeIndex(f, t) - 1], + eft.getTermNodeValueLabel(f, t), eft.getTermNodeVersion(f, t)) eft.setNumberOfLocalNodes(newNodeCount) @@ -130,7 +133,7 @@ def remapEftNodeValueLabelsVersion(eft, localNodeIndexes, valueLabels, version): valueLabel = eft.getTermNodeValueLabel(f, t) if (localNodeIndex in localNodeIndexes) and (valueLabel in valueLabels): result = eft.setTermNodeParameter(f, t, localNodeIndex, valueLabel, version) - #print('remap result', result) + # print('remap result', result) def remapEftNodeValueLabelWithNodes(eft, localNodeIndex, fromValueLabel, expressionTerms): @@ -145,7 +148,8 @@ def remapEftNodeValueLabelWithNodes(eft, localNodeIndex, fromValueLabel, express functionCount = eft.getNumberOfFunctions() for f in range(1, functionCount + 1): if eft.getFunctionNumberOfTerms(f) == 1: - if (eft.getTermLocalNodeIndex(f, 1) == localNodeIndex) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and (not getEftTermScaling(eft, f, 1)): + if (eft.getTermLocalNodeIndex(f, 1) == localNodeIndex) and ( + eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and (not getEftTermScaling(eft, f, 1)): termCount = len(expressionTerms) eft.setFunctionNumberOfTerms(f, termCount) version = eft.getTermNodeVersion(f, 1) @@ -246,7 +250,7 @@ def createEftElementSurfaceLayer(elementIn, eftIn, eftfactory, eftStd, removeNod eft = eftStd scalefactors = None scaleFactorCountIn = eftIn.getNumberOfLocalScaleFactors() - scaleFactorsUsed = [False]*scaleFactorCountIn # flag scale factors used on top surface of element + scaleFactorsUsed = [False] * scaleFactorCountIn # flag scale factors used on top surface of element functionCountIn = eftIn.getNumberOfFunctions() halfFunctionCountIn = functionCountIn // 2 elementBasis = eftfactory.getElementbasis() @@ -272,7 +276,7 @@ def createEftElementSurfaceLayer(elementIn, eftIn, eftfactory, eftStd, removeNod scalingCount, scalingIndexes = eftIn.getTermScaling(fIn, t, 0) if scalingCount > 0: scalingIndexes = eftIn.getTermScaling(fIn, t, scalingCount)[1] - #if (scalingCount > 1) or (scalingIndexes != 1): + # if (scalingCount > 1) or (scalingIndexes != 1): # print("Scaling", scalingIndexes) if scalingCount == 1: scalingIndexes = [scalingIndexes] @@ -306,7 +310,7 @@ def createEftElementSurfaceLayer(elementIn, eftIn, eftfactory, eftStd, removeNod else: newt = t assert noRemap or (newt > 0), "Element " + str(elementIn.getIdentifier()) + " f " + str(f) + \ - " terms " + str(termCountIn) + " terms " + str(termCountIn) if True in scaleFactorsUsed: # get used scale factors, reorder indexes result, scalefactorsIn = elementIn.getScaleFactors(eftIn, scaleFactorCountIn) @@ -478,9 +482,9 @@ def determineTricubicHermiteEft(mesh, nodeParameters, nodeDerivativeFixedWeights functionNumber = n * 8 + crossDerivativeLabels[cd] derivativeIndexes = \ [0, 1] if (cd == 0) else \ - [0, 2] if (cd == 1) else \ - [1, 2] if (cd == 2) else \ - [0, 1, 2] + [0, 2] if (cd == 1) else \ + [1, 2] if (cd == 2) else \ + [0, 1, 2] sign = 1.0 crossDerivativeLabel = Node.VALUE_LABEL_VALUE for ed in derivativeIndexes: @@ -541,6 +545,9 @@ def _determinePermutations(self): # permutations.append((directions[j], directions[i])) else: normDir = [normalize(dir) for dir in self._directions] + # print("directions", self._directions) + # print("normDir", normDir) + # print("directionsCount", directionsCount) for i in range(directionsCount - 2): for j in range(i + 1, directionsCount - 1): if dot(normDir[i], normDir[j]) < -0.9: @@ -570,6 +577,7 @@ def _determinePermutations(self): permutations.append((self._directions[ii], self._directions[jj], self._directions[kk])) permutations.append((self._directions[jj], self._directions[kk], self._directions[ii])) permutations.append((self._directions[kk], self._directions[ii], self._directions[jj])) + # print("permutations", permutations) return permutations def getComplexity(self): @@ -610,11 +618,13 @@ def getDerivativeWeightsList(self, nodeDeltas, nodeDerivatives, localNodeIndex): inwardNodeDeltas = [[-d for d in nodeDeltas[i]] if flips[i] else nodeDeltas[i] for i in swizzleIndexes] derivativeWeightsList = None greatestSimilarity = -1.0 + # print("permutations", self._permutations) for permutation in self._permutations: # skip permutations using directions not in any supplied limitDirections skipPermutation = False limitIndex = 0 for limitDirections in self._limitDirections: + # print("limitDirections", limitDirections) if limitDirections: weights = permutation[swizzleIndexes[limitIndex]] if flips[limitIndex]: @@ -645,9 +655,12 @@ def getDerivativeWeightsList(self, nodeDeltas, nodeDerivatives, localNodeIndex): cosineSimilarity = dot(derivative, delta) / (magDerivative * magDelta) # magnitudeSimilarity = math.exp(-math.fabs((magDerivative - magDelta) / magDelta)) similarity += cosineSimilarity # * magnitudeSimilarity + # print("similarity", similarity) if similarity > greatestSimilarity: greatestSimilarity = similarity derivativeWeightsList = permutation + # print("swizzleIndexes", swizzleIndexes) + # print("derivativeWeightsList", derivativeWeightsList) finalWeightsList = [ [-w for w in derivativeWeightsList[swizzleIndexes[i]]] if flips[i] else derivativeWeightsList[swizzleIndexes[i]] for i in range(derivativesPerNode) @@ -719,19 +732,23 @@ def __init__(self): HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]), HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]), HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.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]]) + self._nodeLayoutTriplePoint = [ + 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]]), + 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]]), + 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]]), + 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._nodeLayoutTriplePoint23Front = HermiteNodeLayout( [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 1.0]]) self._nodeLayoutTriplePoint23Back = HermiteNodeLayout( [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, -1.0]]) self.nodeLayoutsBifurcation6WayTriplePoint = {} + self._nodeLayoutCapShellTriplePoint = [ + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [-1.0, 1.0, 0.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [1.0, 1.0, 0.0]]), + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0], [-1.0, -1.0, 0.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0], [1.0, -1.0, 0.0]])] + self._nodeLayoutCapBoxShield = None + self._nodeLayoutCapBoxShieldTriplePoint = None def getNodeLayoutRegularPermuted(self, d3Defined, limitDirections=None): """ @@ -847,9 +864,7 @@ def getNodeLayoutTriplePoint(self): Bottom Left, and Bottom Right) each with its specific node layout. :return: List of 4 HermiteNodeLayout. """ - nodeLayouts = [self._nodeLayoutTriplePointTopLeft, self._nodeLayoutTriplePointTopRight, - self._nodeLayoutTriplePointBottomLeft, self._nodeLayoutTriplePointBottomRight] - return nodeLayouts + return self._nodeLayoutTriplePoint def getNodeLayoutBifurcation6WayTriplePoint(self, segmentsIn, sequence, maxMajorSegment, top): """ @@ -927,7 +942,7 @@ def getNodeLayoutBifurcation6WayTriplePoint(self, segmentsIn, sequence, maxMajor [0.0, 0.0, 1.0], [0.0, 1.0, -1.0]]) nodeLayoutBifurcationCoreTransitionMMM = HermiteNodeLayout( [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [1.0, 1.0, 0.0], - [0.0, 0.0, 11.0], [-1.0, -1.0, -1.0]]) + [0.0, 0.0, 1.0], [-1.0, -1.0, -1.0]]) nodeLayouts = ( nodeLayoutBifurcationCoreTransitionPOP, nodeLayoutBifurcationCoreTransitionOPP, @@ -953,6 +968,76 @@ def getNodeLayoutBifurcation6WayTriplePoint(self, segmentsIn, sequence, maxMajor return self._nodeLayoutBifurcationCoreTransitionBottomGeneral return nodeLayouts[layoutIndex] + def getNodeLayoutCapTransition(self, isSpecial=False): + """ + Get node layout for transition elements between the core box elements and the shell elements in the cap mesh. + :param isSpecial: True if cap transition is a special case where the number of transition elements > 1 for + box shield elements. False for all other cases. + :return: HermiteNodeLayout. + """ + if isSpecial: + nodeLayoutCapTransition = HermiteNodeLayout( + [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]]) + else: + nodeLayoutCapTransition = self._nodeLayoutRegularPermuted_d3Defined + return nodeLayoutCapTransition + + def getNodeLayoutCapShellTriplePoint(self): + """ + Get node layout for triple-point corners of shell elements in the cap mesh. There are four corners + (Top Left, Top Right, Bottom Left, and Bottom Right) each with its specific node layout. + :return: HermiteNodeLayout. + """ + return self._nodeLayoutCapShellTriplePoint + + def getNodeLayoutCapBoxShield(self, isStartCap=True): + """ + Get node layout for elements between the core box and the shield in the cap mesh. There are four types + (Top, Bottom, Left, and Right) each with its specific node layout. The node layout also changes direction in + the y-axis depending on whether the cap is at the start or the end of a tube segment. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + :return: HermiteNodeLayout. + """ + if isStartCap: + self._nodeLayoutCapBoxShield = [ + HermiteNodeLayout([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0], [0.0, -1.0, 1.0]]), + HermiteNodeLayout([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.0, -1.0, -1.0]]), + HermiteNodeLayout([[1.0, 0.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]]), + HermiteNodeLayout([[-1.0, 0.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]])] + else: + self._nodeLayoutCapBoxShield = [ + HermiteNodeLayout([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 1.0]]), + HermiteNodeLayout([[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 1.0, -1.0]]), + HermiteNodeLayout([[1.0, 0.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]]), + HermiteNodeLayout([[-1.0, 0.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]])] + + return self._nodeLayoutCapBoxShield + + def getNodeLayoutCapBoxShieldTriplePoint(self, isStartCap=True): + """ + Get node layout for triple-point corners of core box elements in the cap mesh. There are four corners + (Top Left, Top Right, Bottom Left, and Bottom Right) each with its specific node layout. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + :return: HermiteNodeLayout. + """ + if isStartCap: + self._nodeLayoutCapBoxShieldTriplePoint = [ + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, -1.0, 1.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]), + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [-1.0, -1.0, -1.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, -1.0, -1.0]])] + else: + self._nodeLayoutCapBoxShieldTriplePoint = [ + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]]), + HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [-1.0, 1.0, -1.0]]), + HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 1.0, -1.0]])] + + return self._nodeLayoutCapBoxShieldTriplePoint + + def determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts): """ Determine the bicubic or tricubic Hermite serendipity element field template for @@ -1034,7 +1119,7 @@ def determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts): nodeParameters[n][1], nodeParameters[n][2], nodeParameters[n][3] if d3Defined else None] - derivativeWeightsList =\ + derivativeWeightsList = \ nodeLayout.getDerivativeWeightsList(deltas[n], nodeDerivatives, n) if nodeLayout else None for ed in range(derivativesPerNode): if nodeLayout: @@ -1128,7 +1213,8 @@ def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodePara :param localNodeIndexes: Local node indexes to scale value label at. Currently must be in [5, 6, 7, 8]. :param valueLabel: Single value label to scale. Currently only implemened for D_DS3. :param version: Set to a number > 1 to use that derivative version instead of scale factors, and client is - responsible for assigning that derivative version in proportion to the returned scale factors. + responsible for assigning that derivative version in proportion to the returned scale factors. Versions 3 and 4 + are used when the cap mesh is active, each representing the version used for start cap and end cap, respectively. :return: Modified eft, final scalefactors, add scalefactors (multiplying the affected derivatives; these are included in final scalefactors if version == 1, otherwise client must use these to scale versioned derivatives). """ @@ -1150,5 +1236,11 @@ def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodePara newScalefactors = scalefactors + addScalefactors if scalefactors else addScalefactors addScaleEftNodesValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS3, 3) return eft, newScalefactors, addScalefactors + elif version in [3, 4]: + addScalefactors = addScalefactors[-2:] if version == 3 else addScalefactors[:2] + localNodeIndexes = [7, 8] if version == 3 else [5, 6] + newScalefactors = scalefactors + addScalefactors if scalefactors else addScalefactors + addScaleEftNodesValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS3, 3) + return eft, newScalefactors, addScalefactors remapEftNodeValueLabelsVersion(eft, localNodeIndexes, [valueLabel], version) return eft, scalefactors, addScalefactors diff --git a/src/scaffoldmaker/utils/networkmesh.py b/src/scaffoldmaker/utils/networkmesh.py index 039c7476..da43e223 100644 --- a/src/scaffoldmaker/utils/networkmesh.py +++ b/src/scaffoldmaker/utils/networkmesh.py @@ -102,15 +102,17 @@ class NetworkSegment: Describes a segment of a network between junctions as a sequence of nodes with node derivative versions. """ - def __init__(self, networkNodes: list, nodeVersions: list, isPatch): + def __init__(self, networkNodes: list, nodeVersions: list, isCap: list, isPatch): """ :param networkNodes: List of NetworkNodes from start to end. Must be at least 2. :param nodeVersions: List of node versions to use for derivatives at network nodes. + :param isCap: List of boolean true if segment requires a cap at either ends. [Start, End] :param isPatch: True if segment at the other end of the junction requires a patch. """ assert isinstance(networkNodes, list) and (len(networkNodes) > 1) and (len(nodeVersions) == len(networkNodes)) self._networkNodes = networkNodes self._nodeVersions = nodeVersions + self._isCap = isCap self._isPatch = isPatch self._elementIdentifiers = [None] * (len(networkNodes) - 1) for networkNode in networkNodes[1:-1]: @@ -162,6 +164,12 @@ def isCyclic(self): """ return False # not implemented, assume not cyclic + def isCap(self): + """ + :return: True if the segment requires a cap mesh, False if not. + """ + return self._isCap + def isPatch(self): """ :return: True if the segment is a patch, False if not. @@ -216,16 +224,24 @@ def build(self, structureString): self._networkSegments = [] sequenceStrings = structureString.split(",") for sequenceString in sequenceStrings: - # check if segment is a patch - if not sequenceString[0].isnumeric(): - try: - isPatch = True if sequenceString[0] == "#" else False - sequenceString = sequenceString[2:] if isPatch else sequenceString - except ValueError: - print("Network mesh: Skipping invalid cap sequence", sequenceString, file=sys.stderr) - continue - else: - isPatch = False + # check if the segment requires a cap or a patch + isStartCap = isEndCap = isPatch = False + try: + # Check and handle caps + if sequenceString.startswith("("): + isStartCap = True + sequenceString = sequenceString[1:] + if sequenceString.endswith(")"): + isEndCap = True + sequenceString = sequenceString[:-1] + # Check and handle patch + if sequenceString.startswith("#"): + isPatch = True + sequenceString = sequenceString[2:] + except (ValueError, IndexError): + print("Network mesh: Skipping invalid cap sequence", sequenceString, file=sys.stderr) + continue + isCap = [isStartCap, isEndCap] nodeIdentifiers = [] nodeVersions = [] @@ -260,7 +276,7 @@ def build(self, structureString): sequenceNodes.append(networkNode) sequenceVersions.append(nodeVersion) if (len(sequenceNodes) > 1) and (existingNetworkNode or (nodeIdentifier == nodeIdentifiers[-1])): - networkSegment = NetworkSegment(sequenceNodes, sequenceVersions, isPatch) + networkSegment = NetworkSegment(sequenceNodes, sequenceVersions, isCap, isPatch) self._networkSegments.append(networkSegment) sequenceNodes = sequenceNodes[-1:] sequenceVersions = sequenceVersions[-1:] diff --git a/src/scaffoldmaker/utils/tubenetworkmesh.py b/src/scaffoldmaker/utils/tubenetworkmesh.py index 224b503e..24998f6a 100644 --- a/src/scaffoldmaker/utils/tubenetworkmesh.py +++ b/src/scaffoldmaker/utils/tubenetworkmesh.py @@ -4,6 +4,8 @@ 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.capmesh import CapMesh from scaffoldmaker.utils.eft_utils import ( addTricubicHermiteSerendipityEftParameterScaling, determineCubicHermiteSerendipityEft, HermiteNodeLayoutManager) from scaffoldmaker.utils.interpolation import ( @@ -80,6 +82,12 @@ def __init__(self, region, meshDimension, coordinateFieldName="coordinates", d3Defined, limitDirections=[None, [[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], None]) self._nodeLayoutTransitionTriplePoint = None + self._nodeLayoutCapTransition = None + self._nodeLayoutCapShellTriplePoint = None + self._nodeLayoutCapBoxShield = None + self._nodeLayoutCapBoxShieldTriplePoint = None + + # annotation groups created if core: self._coreGroup = None self._shellGroup = None @@ -88,6 +96,14 @@ def __init__(self, region, meshDimension, coordinateFieldName="coordinates", self._rightGroup = None self._dorsalGroup = None self._ventralGroup = None + # annotations for the kidney scaffold: + self._renalCapsuleGroup = None + self._lateralGroup = None + self._medialGroup = None + self._anteriorGroup = None + self._posteriorGroup = None + self._openingGroup = None + def getStandardElementtemplate(self): return self._standardElementtemplate, self._standardEft @@ -98,6 +114,9 @@ def createElementfieldtemplate(self): """ return self._mesh.createElementfieldtemplate(self._elementbasis) + def getCapElementtemplate(self): + return self._capElementtemplate + def getNodeLayout5Way(self): return self._nodeLayout5Way @@ -168,6 +187,92 @@ def getNodeLayoutBifurcation6WayTriplePoint(self, segmentsIn, sequence, maxMajor return self._nodeLayoutManager.getNodeLayoutBifurcation6WayTriplePoint( segmentsIn, sequence, maxMajorSegment, top) + def getNodeLayoutCapTransition(self, isSpecial=False): + """ + Node layout for generating cap transition elements, excluding at triple points. + :param isSpecial:True if cap transition is a special case where the number of transition elements > 1 for + box shield elements. False for all other cases. + :return: Node layout. + """ + self._nodeLayoutCapTransition = self._nodeLayoutManager.getNodeLayoutCapTransition(isSpecial=isSpecial) + return self._nodeLayoutCapTransition + + def getNodeLayoutCapShellTriplePoint(self, location): + """ + Special node layout for generating cap shell 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.getNodeLayoutCapShellTriplePoint() + 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._nodeLayoutCapTransition + + self._nodeLayoutCapShellTriplePoint = nodeLayout + return self._nodeLayoutCapShellTriplePoint + + def getNodeLayoutCapBoxShield(self, location, isStartCap=True): + """ + Special node layout for generating cap box-shield transition elements. + There are four layouts relative to the core box: Top (location = 1); bottom (location = -1); + left (location = 2); and right (location = -2). + :param location: Location identifier identifying four corners of solid core box. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + :return: Node layout. + """ + nodeLayouts = self._nodeLayoutManager.getNodeLayoutCapBoxShield(isStartCap) + assert location in [1, -1, 2, -2, 0] + if location == 1: # "Top" + nodeLayout = nodeLayouts[0] + elif location == -1: # "Bottom" + nodeLayout = nodeLayouts[1] + elif location == 2: # "Left" + nodeLayout = nodeLayouts[2] + elif location == -2: # "Right" + nodeLayout = nodeLayouts[3] + else: + nodeLayout = None + + self._nodeLayoutCapBoxShield = nodeLayout + return self._nodeLayoutCapBoxShield + + def getNodeLayoutCapBoxShieldTriplePoint(self, location, isStartCap=True): + """ + Special node layout for generating cap box-shield transition elements at triple points. + There are four layouts relative to 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. + :param isStartCap: True if generating a cap mesh at the start of a tube segment, False if generating at the end + of a tube segment. + :return: Node layout. + """ + nodeLayouts = self._nodeLayoutManager.getNodeLayoutCapBoxShieldTriplePoint(isStartCap) + 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._nodeLayoutCapBoxShield + + self._nodeLayoutCapBoxShieldTriplePoint = nodeLayout + return self._nodeLayoutCapBoxShieldTriplePoint + def getNodetemplate(self): return self._nodetemplate @@ -207,6 +312,36 @@ def getVentralMeshGroup(self): self._ventralGroup = self.getOrCreateAnnotationGroup(("ventral", "")) return self._ventralGroup.getMeshGroup(self._mesh) + def getRenalCapsuleMeshGroup(self): + if not self._renalCapsuleGroup: + self._renalCapsuleGroup = self.getOrCreateAnnotationGroup(("renal capsule", "")) + return self._renalCapsuleGroup.getMeshGroup(self._mesh) + + def getAnteriorMeshGroup(self): + if not self._anteriorGroup: + self._anteriorGroup = self.getOrCreateAnnotationGroup(("anterior", "")) + return self._anteriorGroup.getMeshGroup(self._mesh) + + def getPosteriorMeshGroup(self): + if not self._posteriorGroup: + self._posteriorGroup = self.getOrCreateAnnotationGroup(("posterior", "")) + return self._posteriorGroup.getMeshGroup(self._mesh) + + def getLateralMeshGroup(self): + if not self._lateralGroup: + self._lateralGroup = self.getOrCreateAnnotationGroup(("lateral", "")) + return self._lateralGroup.getMeshGroup(self._mesh) + + def getMedialMeshGroup(self): + if not self._medialGroup: + self._medialGroup = self.getOrCreateAnnotationGroup(("medial", "")) + return self._medialGroup.getMeshGroup(self._mesh) + + def getOpeningMeshGroup(self): + if not self._openingGroup: + self._openingGroup = self.getOrCreateAnnotationGroup(("opening", "")) + return self._openingGroup.getMeshGroup(self._mesh) + def getNewTrimAnnotationGroup(self): self._trimAnnotationGroupCount += 1 return self.getOrCreateAnnotationGroup(("trim surface " + "{:03d}".format(self._trimAnnotationGroupCount), "")) @@ -221,10 +356,12 @@ def resolveEftCoreBoundaryScaling(self, eft, scalefactors, nodeParameters, nodeI x, d1, d2, d3 each with 3 components. :param nodeIdentifiers: List over 8 3-D local nodes giving global node identifiers. :param mode: 1 to set scale factors, 2 to add version 2 to d3 for the boundary nodes and assigning - values to that version equal to the scale factors x version 1. + values to that version equal to the scale factors x version 1. Modes 3 and 4 are used for cap mesh, which are + similar to mode 1, but the local node indexes are limited to [7, 8] for mode 3 (start cap) and [5, 6] for mode + 4 (end cap). :return: New eft, new scalefactors. """ - assert mode in (1, 2) + assert mode in (1, 2, 3, 4) eft, scalefactors, addScalefactors = addTricubicHermiteSerendipityEftParameterScaling( eft, scalefactors, nodeParameters, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS3, version=mode) if mode == 2: @@ -304,10 +441,27 @@ def __init__(self, networkSegment, pathParametersList, elementsCountAround, elem # [nAlong][nAcrossMajor][nAcrossMinor] format. self._boxElementIds = None # [along][major][minor] + self._boxExtCoordinates = None + self._transitionExtCoordinates = None + self.rimExtCoordinates = None + + self._networkPathParameters = pathParametersList + self._isCap = networkSegment.isCap() + if self._isCap: + if not self._isCore: + self._getNumberOfBoxElements() + self._elementsCountTransition = 1 + self._capMesh = CapMesh(self._elementsCountAround, self._elementsCountCoreBoxMajor, + self._elementsCountCoreBoxMinor, self._elementsCountThroughShell, + self._elementsCountTransition, self._networkPathParameters, + self._isCap, self._isCore) + def getNetworkSegment(self): return self._networkSegment def getCoreBoundaryScalingMode(self): + modes = (1, 2, 3, 4) if self._isCap else (1, 2) + assert self._coreBoundaryScalingMode in modes return self._coreBoundaryScalingMode def getElementsCountAround(self): @@ -319,6 +473,12 @@ def getRawTubeCoordinates(self, pathIndex=0): def getIsCore(self): return self._isCore + def getIsCap(self): + return self._isCap + + def getCapMesh(self): + return self._capMesh + def getElementsCountCoreBoxMajor(self): return self._elementsCountCoreBoxMajor @@ -534,10 +694,30 @@ def sample(self, fixedElementsCountAlong, targetElementLength): self._rimElementIds = [None] * elementsCountAlong self._boxElementIds = [None] * elementsCountAlong - if self._isCore: - # sample coordinates for the solid core + # sample coordinates for the solid core + if self._dimension == 3: self._sampleCoreCoordinates(elementsCountAlong) + if self._isCap: + # sample coordinates for the cap mesh at the ends of a tube segment + self._boxExtCoordinates, self._transitionExtCoordinates, self.rimExtCoordinates = ( + self._capMesh.sampleCoordinates(self._boxCoordinates, self._transitionCoordinates, self._rimCoordinates)) + # if self._isCore: + self._smoothD2DerivativesAtCapTubeJoint() + + def _getNumberOfBoxElements(self): + """ + Calculates the number of core box elements required to form a square-looking shield mesh. + Only used when the tube is without a core. + """ + # assert + half = self._elementsCountAround // 4 + if half % 2 == 0: + self._elementsCountCoreBoxMajor = self._elementsCountCoreBoxMinor = half + else: + self._elementsCountCoreBoxMinor = half - 1 if (half - 1) % 2 == 0 else half - 3 + self._elementsCountCoreBoxMajor = half + 1 if (half + 1) % 2 == 0 else half + 3 + def _sampleCoreCoordinates(self, elementsCountAlong): """ Black box function for sampling coordinates for the solid core. @@ -1053,6 +1233,57 @@ def _createBoxBoundaryNodeIdsList(self, startSkipCount=None, endSkipCount=None): return boxBoundaryNodeIds, boxBoundaryNodeToBoxId + def _smoothD2DerivativesAtCapTubeJoint(self): + """ + Smooths D2 derivatives at the joint where the cap and the tube surfaces join to eliminate zero Jacobian contours. + """ + def smoothBoxDerivatives(): + capCoordinates = self._boxExtCoordinates[i] + for m in range(self._elementsCountCoreBoxMajor + 1): + for n in range(self._elementsCountCoreBoxMinor + 1): + nx = [capCoordinates[0][m][n], self._boxCoordinates[0][i][m][n]] + nd2 = [capCoordinates[2][m][n], self._boxCoordinates[2][i][m][n]] + if i == -1: + nx.reverse() + nd2.reverse() + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2) + self._boxCoordinates[2][i][m][n] = sd2[1] + + def smoothRimDerivatives(): + capCoordinates = self.rimExtCoordinates[i] + for n3 in range(self._elementsCountThroughShell + 1): + for n1 in range(self._elementsCountAround): + nx = [capCoordinates[0][n3][n1], self._rimCoordinates[0][i][n3][n1]] + nd2 = [capCoordinates[2][n3][n1], self._rimCoordinates[2][i][n3][n1]] + if i == -1: + nx.reverse() + nd2.reverse() + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2) + self._rimCoordinates[2][i][n3][n1] = sd2[1] + + if self._transitionExtCoordinates is not None: + capCoordinates = self._transitionExtCoordinates[i] + for n3 in range(self._elementsCountTransition - 1): + for n1 in range(self._elementsCountAround): + nx = [capCoordinates[0][n3][n1], self._transitionCoordinates[0][i][n3][n1]] + nd2 = [capCoordinates[2][n3][n1], self._transitionCoordinates[2][i][n3][n1]] + if i == -1: + nx.reverse() + nd2.reverse() + sd2 = smoothCubicHermiteDerivativesLine(nx, nd2) + self._transitionCoordinates[2][i][n3][n1] = sd2[1] + + for i in [0, -1]: + if not self._isCap[i]: + continue + + if self._isCore: + smoothBoxDerivatives() + smoothRimDerivatives() + else: + smoothRimDerivatives() + + @classmethod def blendSampledCoordinates(cls, segment1, nodeIndexAlong1, segment2, nodeIndexAlong2): nodesCountAround = segment1._elementsCountAround @@ -1363,19 +1594,23 @@ def setBoxElementId(self, e1, e2, e3, elementIdentifier): [None] * self._elementsCountCoreBoxMinor for _ in range(self._elementsCountCoreBoxMajor)] self._boxElementIds[e2][e3][e1] = elementIdentifier - def _addBoxElementsToMeshGroup(self, e1Start, e1Limit, e3Start, e3Limit, meshGroup): + def _addBoxElementsToMeshGroup(self, e1Start, e1Limit, e3Start, e3Limit, meshGroup, e2Start=None, e2Limit=None): """ Add ranges of box elements to mesh group. :param e1Start: Start element index in major / d2 direction. :param e1Limit: Limit element index in major / d2 direction. :param e3Start: Start element index in minor / d3 direction. :param e3Limit: Limit element index in minor / d3 direction. + :param e2Start: Start element index in d1 direction. If None, use 0 as default. + :param e2Limit: Limit element index in d1 direction. If None, use elementsCountAlong as default. :param meshGroup: Zinc MeshGroup to add elements to. """ # print("Add box elements", e1Start, e1Limit, e3Start, e3Limit, meshGroup.getName()) elementsCountAlong = self.getSampledElementsCountAlong() mesh = meshGroup.getMasterMesh() - for e2 in range(elementsCountAlong): + e2Start = 0 if e2Start is None else e2Start + e2Limit = elementsCountAlong if e2Limit is None else e2Limit + for e2 in range(e2Start, e2Limit): boxSlice = self._boxElementIds[e2] if boxSlice: # print(boxSlice[e1Start:e1Limit]) @@ -1384,19 +1619,23 @@ def _addBoxElementsToMeshGroup(self, e1Start, e1Limit, e3Start, e3Limit, meshGro element = mesh.findElementByIdentifier(elementIdentifier) meshGroup.addElement(element) - def _addRimElementsToMeshGroup(self, e1Start, e1Limit, e3Start, e3Limit, meshGroup): + def _addRimElementsToMeshGroup(self, e1Start, e1Limit, e3Start, e3Limit, meshGroup, e2Start=None, e2Limit=None): """ Add ranges of rim elements to mesh group. :param e1Start: Start element index around. Can be negative which supports wrapping. :param e1Limit: Limit element index around. :param e3Start: Start element index rim. :param e3Limit: Limit element index rim. + :param e2Start: Start element index in d1 direction. If None, use 0 as default. + :param e2Limit: Limit element index in d1 direction. If None, use elementsCountAlong as default. :param meshGroup: Zinc MeshGroup to add elements to. """ # print("Add rim elements", e1Start, e1Limit, e3Start, e3Limit, meshGroup.getName()) elementsCountAlong = self.getSampledElementsCountAlong() mesh = meshGroup.getMasterMesh() - for e2 in range(elementsCountAlong): + e2Start = 0 if e2Start is None else e2Start + e2Limit = elementsCountAlong if e2Limit is None else e2Limit + for e2 in range(e2Start, e2Limit): rimSlice = self._rimElementIds[e2] if rimSlice: for elementIdentifiersList in rimSlice[e3Start:e3Limit]: @@ -1438,6 +1677,23 @@ def addAllElementsToMeshGroup(self, meshGroup): self.addCoreElementsToMeshGroup(meshGroup) self.addShellElementsToMeshGroup(meshGroup) + def addSideD1ElementsToMeshGroup(self, side: bool, meshGroup): + """ + Add elements to the mesh group on side of +d1 or -d1, often matching anterior and posterior. + :param side: False for +d1 direction, True for -d1 direction. + :param meshGroup: Zinc MeshGroup to add elements to. + """ + elementsCountAlong = self.getSampledElementsCountAlong() + e2Start = 0 if side else elementsCountAlong // 2 + e2Limit = elementsCountAlong // 2 if side else elementsCountAlong + if self._isCore: + self._addBoxElementsToMeshGroup(0, self._elementsCountCoreBoxMajor, + 0, self._elementsCountCoreBoxMinor, meshGroup, + e2Start=e2Start, e2Limit=e2Limit) + self._addRimElementsToMeshGroup(0, self._elementsCountAround, + 0, self.getElementsCountRim(), meshGroup, + e2Start=e2Start, e2Limit=e2Limit) + def addSideD2ElementsToMeshGroup(self, side: bool, meshGroup): """ Add elements to the mesh group on side of +d2 or -d2, often matching left and right. @@ -1470,6 +1726,19 @@ def addSideD3ElementsToMeshGroup(self, side: bool, meshGroup): e1Limit = e1Start + (self._elementsCountAround // 2) self._addRimElementsToMeshGroup(e1Start, e1Limit, 0, self.getElementsCountRim(), meshGroup) + def addShellOpeningElementsToMeshGroup(self, e1Start, e1Limit, meshGroup): + """ + Add specific elements in the shell to mesh group. + :param e1Start: Start element index around. + :param e1Limit: Limit element index around. + :param meshGroup: Zinc MeshGroup to add elements to. + """ + elementsCountRim = self.getElementsCountRim() + elementsCountShell = self._elementsCountThroughShell + e3ShellStart = elementsCountRim - elementsCountShell + + self._addRimElementsToMeshGroup(e1Start, e1Limit, e3ShellStart, elementsCountRim, meshGroup) + def getRimNodeIdsSlice(self, n2): """ Get slice of rim node IDs. @@ -1547,8 +1816,13 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): self._rimNodeIds[n2] = rimNodeIds continue + capMesh = self._capMesh + # create cap nodes at the start section of a tube segment + if self._isCap[0] and n2 == 0: + capMesh.generateNodes(generateData, isStartCap=True) + # create core box nodes - if self._boxCoordinates: + if self._isCore: self._boxNodeIds[n2] = [] if self._boxNodeIds[n2] is None else self._boxNodeIds[n2] coreBoxMajorNodesCount = self.getCoreBoxMajorNodesCount() coreBoxMinorNodesCount = self.getCoreBoxMinorNodesCount() @@ -1598,6 +1872,10 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): ringNodeIds.append(nodeIdentifier) self._rimNodeIds[n2].append(ringNodeIds) + # create cap nodes at the end section of a tube segment + if self._isCap[-1] and n2 == elementsCountAlong: + self._endCapNodeIds = capMesh.generateNodes(generateData, isStartCap=False) + # create a new list containing box node ids are located at the boundary if self._isCore: self._boxBoundaryNodeIds, self._boxBoundaryNodeToBoxId = ( @@ -1614,6 +1892,10 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): self._boxElementIds[e2] = [] self._rimElementIds[e2] = [] e2p = e2 + 1 + # create cap elements at the start of the tube + if self._isCap[0] and e2 == 0: + capMesh.generateElements(self._boxNodeIds, self._rimNodeIds, annotationMeshGroups, isStartCap=True) + if self._isCore: # create box elements elementsCountAcrossMinor = self.getCoreBoxMinorNodesCount() - 1 @@ -1710,6 +1992,10 @@ def generateMesh(self, generateData: TubeNetworkMeshGenerateData, n2Only=None): ringElementIds.append(elementIdentifier) self._rimElementIds[e2].append(ringElementIds) + # create cap elements at the end of the tube + if self._isCap[-1] and e2 == (elementsCountAlong - endSkipCount - 1): + capMesh.generateElements(self._boxNodeIds, self._rimNodeIds, annotationMeshGroups, isStartCap=False) + def generateJunctionRimElements(self, junction, generateData): """ Generates rim elements for junction part of the segment after segment nodes and elements have been made. @@ -4052,6 +4338,12 @@ def createSegment(self, networkSegment): self._layoutNodes, self._layoutInnerCoordinates, pathValueLabels, networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions())) elementsCountAround = self._defaultElementsCountAround + if any(networkSegment.isCap()) and not self._isCore: + if elementsCountAround < 8: + elementsCountAround = 8 + elif elementsCountAround % 4: + elementsCountAround += 4 - self._defaultElementsCountAround % 4 + elementsCountCoreBoxMinor = self._defaultElementsCountCoreBoxMinor coreBoundaryScalingMode = self._defaultCoreBoundaryScalingMode @@ -4106,8 +4398,14 @@ def generateMesh(self, generateData): shellMeshGroup = generateData.getShellMeshGroup() for networkSegment in self._networkMesh.getNetworkSegments(): segment = self._segments[networkSegment] + segmentCaps = segment.getIsCap() segment.addCoreElementsToMeshGroup(coreMeshGroup) segment.addShellElementsToMeshGroup(shellMeshGroup) + for isCap in segmentCaps: + if isCap: + capMesh = segment.getCapMesh() + capMesh.addCapBoxElementsToMeshGroup(coreMeshGroup) + capMesh.addCapShellElementsToMeshGroup(shellMeshGroup) class BodyTubeNetworkMeshBuilder(TubeNetworkMeshBuilder): @@ -4128,6 +4426,11 @@ def generateMesh(self, generateData): ventralMeshGroup = generateData.getVentralMeshGroup() for networkSegment in self._networkMesh.getNetworkSegments(): segment = self._segments[networkSegment] + segmentCaps = segment.getIsCap() + if True in segmentCaps: + capMesh = segment.getCapMesh() + else: + capMesh = None annotationTerms = segment.getAnnotationTerms() for annotationTerm in annotationTerms: if "left" in annotationTerm[0]: @@ -4140,8 +4443,127 @@ def generateMesh(self, generateData): # segment on main axis segment.addSideD2ElementsToMeshGroup(False, leftMeshGroup) segment.addSideD2ElementsToMeshGroup(True, rightMeshGroup) + if capMesh: + capMesh.addCapSideD2ElementsToMeshGroup(False, leftMeshGroup) + capMesh.addCapSideD2ElementsToMeshGroup(True, rightMeshGroup) segment.addSideD3ElementsToMeshGroup(False, ventralMeshGroup) segment.addSideD3ElementsToMeshGroup(True, dorsalMeshGroup) + if capMesh: + capMesh.addCapSideD3ElementsToMeshGroup(False, ventralMeshGroup) + capMesh.addCapSideD3ElementsToMeshGroup(True, dorsalMeshGroup) + + +class KidneyTubeNetworkMeshBuilder(TubeNetworkMeshBuilder): + """ + Specialization of TubeNetworkMeshBuilder adding annotations for anterior, posterior, lateral, medial, and hilum regions. + """ + + def __init__(self, networkMesh: NetworkMesh, targetElementDensityAlongLongestSegment: float, + layoutAnnotationGroups: list=[], annotationElementsCountsAlong: list=[], + defaultElementsCountAround: int=8, annotationElementsCountsAround: list=[], + elementsCountThroughShell: int=1, isCore=False, elementsCountTransition: int=1, + defaultElementsCountCoreBoxMinor: int=2, annotationElementsCountsCoreBoxMinor: list=[], + defaultCoreBoundaryScalingMode=1, annotationCoreBoundaryScalingMode=[], + useOuterTrimSurfaces=True, showKidneys=[]): + """ + Builds specialized continuous tube network meshes for kidney scaffold. + :param showKidneys: List of flags for showing left and/or right kidneys. + """ + super(KidneyTubeNetworkMeshBuilder, self).__init__( + networkMesh, targetElementDensityAlongLongestSegment, layoutAnnotationGroups, annotationElementsCountsAlong, + defaultElementsCountAround, annotationElementsCountsAround, elementsCountThroughShell, isCore, + elementsCountTransition, defaultElementsCountCoreBoxMinor, annotationElementsCountsCoreBoxMinor, + defaultCoreBoundaryScalingMode, annotationCoreBoundaryScalingMode, useOuterTrimSurfaces) + + self._showKidneys = showKidneys + + + def generateMesh(self, generateData): + super(KidneyTubeNetworkMeshBuilder, self).generateMesh(generateData) + # build anterior, posterior, lateral, medial annotation groups + anteriorMeshGroup = generateData.getAnteriorMeshGroup() + posteriorMeshGroup = generateData.getPosteriorMeshGroup() + lateralMeshGroup = generateData.getLateralMeshGroup() + medialMeshGroup = generateData.getMedialMeshGroup() + dorsalMeshGroup = generateData.getDorsalMeshGroup() + ventralMeshGroup = generateData.getVentralMeshGroup() + openingMeshGroup = generateData.getOpeningMeshGroup() + + elementsCountAround = self._defaultElementsCountAround + halfElementsCountAround = elementsCountAround // 2 + increment = max(1, elementsCountAround // 8) + + leftKidney, rightKidney = 0, 1 + # Kidney configuration mapping + kidney_configs = { + leftKidney: { + 'lateral_flag': False, + 'medial_flag': True, + 'e1_start': halfElementsCountAround - increment, + 'e1_end': halfElementsCountAround + increment + }, + rightKidney: { + 'lateral_flag': True, + 'medial_flag': False, + 'e1_start': -increment, + 'e1_end': increment + } + } + + def add_common_elements(mesh_obj, method_prefix=""): + """ + Add D1 and D3 elements that are common to all kidneys. + :param mesh_obj: The mesh object (segment or capMesh) to add elements to + :param method_prefix: Prefix for method names (e.g., "addCap" for capMesh methods) + """ + d1_method = f"{method_prefix}SideD1ElementsToMeshGroup" + d3_method = f"{method_prefix}SideD3ElementsToMeshGroup" + + getattr(mesh_obj, d1_method)(False, anteriorMeshGroup) + getattr(mesh_obj, d1_method)(True, posteriorMeshGroup) + getattr(mesh_obj, d3_method)(False, ventralMeshGroup) + getattr(mesh_obj, d3_method)(True, dorsalMeshGroup) + + def add_kidney_specific_elements(mesh_obj, kidney, method_prefix=""): + """ + Add kidney-specific D2 elements. + :param mesh_obj: The mesh object (segment or capMesh) to add elements to + :param kidney: Kidney identifier (leftKidney=0, rightKidney=1) + :param method_prefix: Prefix for method names (e.g., "addCap" for capMesh methods) + """ + if kidney not in kidney_configs: + return + + config = kidney_configs[kidney] + d2_method = f"{method_prefix}SideD2ElementsToMeshGroup" + + getattr(mesh_obj, d2_method)(config['lateral_flag'], lateralMeshGroup) + getattr(mesh_obj, d2_method)(config['medial_flag'], medialMeshGroup) + + # Shell opening elements only for segment (not capMesh) + if "Cap" not in method_prefix: + mesh_obj.addShellOpeningElementsToMeshGroup(config['e1_start'], config['e1_end'], openingMeshGroup) + + for kidney in [leftKidney, rightKidney]: + if not self._showKidneys[kidney]: + continue + + idx = 0 if False in self._showKidneys else kidney + networkSegment = self._networkMesh.getNetworkSegments()[idx] + segment = self._segments[networkSegment] + segmentCaps = segment.getIsCap() + capMesh = segment.getCapMesh() if True in segmentCaps else None + + # Apply common elements to segment + add_common_elements(segment, "add") + + # Apply kidney-specific elements to segment + add_kidney_specific_elements(segment, kidney, "add") + + # Apply elements to capMesh if it exists + if capMesh: + add_common_elements(capMesh, "addCap") + add_kidney_specific_elements(capMesh, kidney, "addCap") class TubeEllipseGenerator: diff --git a/tests/test_capmesh.py b/tests/test_capmesh.py new file mode 100644 index 00000000..b407414d --- /dev/null +++ b/tests/test_capmesh.py @@ -0,0 +1,400 @@ +import unittest + +from cmlibs.utils.zinc.finiteelement import evaluateFieldNodesetRange +from cmlibs.utils.zinc.general import ChangeManager +from cmlibs.utils.zinc.group import identifier_ranges_to_string, mesh_group_add_identifier_ranges, \ + mesh_group_to_identifier_ranges +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 findAnnotationGroupByName +from scaffoldmaker.meshtypes.meshtype_3d_tubenetwork1 import MeshType_3d_tubenetwork1 +from scaffoldmaker.scaffoldpackage import ScaffoldPackage + +from testutils import assertAlmostEqualList + +class CapScaffoldTestCase(unittest.TestCase): + + def test_3d_cap_tube_network_default(self): + """ + Test default 3-D tube network with cap at both ends without core is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Default") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + # change the network layout to have cap at both ends of the tube + networkLayoutSettings["Structure"] = "(1-2)" + self.assertEqual("(1-2)", networkLayoutSettings["Structure"]) + + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(13, len(settings)) + self.assertEqual(8, settings["Number of elements around"]) + self.assertEqual(1, settings["Number of elements through shell"]) + self.assertEqual([0], settings["Annotation numbers of elements around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertEqual([0], settings["Annotation numbers of elements along"]) + self.assertFalse(settings["Use linear through shell"]) + self.assertTrue(settings["Use outer trim surfaces"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(2, settings["Number of elements across core box minor"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation numbers of elements across core box minor"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(8 * 4 + (8 * 2 + 4) * 2, mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual((8 * 5 + 8 * 2 * 2 + 2) * 2, 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.150, -0.100, -0.100], X_TOL) + assertAlmostEqualList(self, maximums, [1.150, 0.100, 0.100], X_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(4) + 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, 0.014392962237699036, delta=X_TOL) + self.assertAlmostEqual(outerSurfaceArea, 0.8148830981004552, delta=X_TOL) + self.assertAlmostEqual(innerSurfaceArea, 0.6329726542126899, delta=X_TOL) + + def test_3d_cap_tube_network_default_core(self): + """ + Test default 3-D tube network with cap at both ends and with a solid core is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Default") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + # change the network layout to have cap at both ends of the tube + networkLayoutSettings["Structure"] = "(1-2)" + self.assertEqual("(1-2)",networkLayoutSettings["Structure"]) + + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(13, len(settings)) + self.assertEqual(8, settings["Number of elements around"]) + self.assertEqual(1, settings["Number of elements through shell"]) + self.assertEqual([0], settings["Annotation numbers of elements around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertEqual([0], settings["Annotation numbers of elements along"]) + self.assertFalse(settings["Use linear through shell"]) + self.assertTrue(settings["Use outer trim surfaces"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(2, settings["Number of elements across core box minor"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation numbers of elements across core box minor"]) + 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 + 8 + 4) * 4 + ((8 + 8 + 4) * 2 + 4 + 4) * 2, mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual((8 + 8 + 9) * 5 + ((8 + 8 + 9) + 9 * 3) * 2, 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.150, -0.100, -0.100], X_TOL) + assertAlmostEqualList(self, maximums, [1.150, 0.100, 0.100], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + 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(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.03862588577889281, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 0.8148830981004554, delta=X_TOL) + + def test_3d_cap_tube_network_bifurcation(self): + """ + Test bifurcation 3-D tube network with cap at both ends is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Bifurcation") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + # change the network layout to have cap at both ends of the tube + networkLayoutSettings["Structure"] = "(1-2.1,2.2-3),2.3-4)" + self.assertEqual("(1-2.1,2.2-3),2.3-4)",networkLayoutSettings["Structure"]) + + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(13, len(settings)) + self.assertEqual(8, settings["Number of elements around"]) + self.assertEqual(1, settings["Number of elements through shell"]) + self.assertEqual([0], settings["Annotation numbers of elements around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertEqual([0], settings["Annotation numbers of elements along"]) + self.assertFalse(settings["Use linear through shell"]) + self.assertTrue(settings["Use outer trim surfaces"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(2, settings["Number of elements across core box minor"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation numbers of elements across core box minor"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual((8 * 4 + 8 * 2 + 4) * 3 , mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual((8 * 4 * 3 + 3 * 3 + 2) * 2 + (8 * 2 * 2 + 2) * 3, 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.1500000000000000, -0.6172290095800493, -0.1000000000000000], X_TOL) + assertAlmostEqualList(self, maximums, [2.1395896893550477, 0.6172290095800493, 0.1000000000000000], X_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(4) + 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, 0.039564164189831115, delta=X_TOL) + self.assertAlmostEqual(outerSurfaceArea, 2.2086099452749344, delta=X_TOL) + self.assertAlmostEqual(innerSurfaceArea, 1.755307334843727, delta=X_TOL) + + def test_3d_tube_network_bifurcation_core(self): + """ + Test bifurcation 3-D tube network with solid core and cap at both ends is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Bifurcation") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + # change the network layout to have cap at both ends of the tube + networkLayoutSettings["Structure"] = "(1-2.1,2.2-3),2.3-4)" + self.assertEqual("(1-2.1,2.2-3),2.3-4)",networkLayoutSettings["Structure"]) + + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(13, len(settings)) + self.assertEqual(8, settings["Number of elements around"]) + self.assertEqual(1, settings["Number of elements through shell"]) + self.assertEqual([0], settings["Annotation numbers of elements around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertEqual([0], settings["Annotation numbers of elements along"]) + self.assertFalse(settings["Use linear through shell"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(2, settings["Number of elements across core box minor"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation numbers of elements across core box minor"]) + 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) +((8 + 8 + 4) * 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) + ((8 + 8 + 9) + 9 * 3) * 3, 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.1500000000000000, -0.6172290095800493, -0.1000000000000000], X_TOL) + assertAlmostEqualList(self, maximums, [2.1395896893550477, 0.6172290095800493, 0.1000000000000000], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + 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(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.1096893209534004, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 2.2086099452749344, delta=X_TOL) + + def test_3d_tube_network_converging_bifurcation_core(self): + """ + Test converging bifurcation 3-D tube network with solid core and 12, 12, 8 elements around. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Bifurcation") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + # change the network layout to have cap at both ends of the tube + networkLayoutSettings["Structure"] = "(1-3.1,(2-3.2,3.3-4)" + self.assertEqual("(1-3.1,(2-3.2,3.3-4)",networkLayoutSettings["Structure"]) + + self.assertTrue(networkLayoutSettings["Define inner coordinates"]) + self.assertEqual(13, len(settings)) + self.assertEqual(8, settings["Number of elements around"]) + self.assertEqual(1, settings["Number of elements through shell"]) + self.assertEqual([0], settings["Annotation numbers of elements around"]) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + self.assertEqual([0], settings["Annotation numbers of elements along"]) + self.assertFalse(settings["Use linear through shell"]) + self.assertFalse(settings["Show trim surfaces"]) + self.assertFalse(settings["Core"]) + self.assertEqual(2, settings["Number of elements across core box minor"]) + self.assertEqual(1, settings["Number of elements across core transition"]) + self.assertEqual([0], settings["Annotation numbers of elements across core box minor"]) + settings["Core"] = True + settings["Number of elements around"] = 12 + settings["Annotation numbers of elements around"] = [8] + + 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(("segment 3", "SEGMENT:3")) + group = annotationGroup1.getGroup() + mesh1d = tmpFieldmodule.findMeshByDimension(1) + meshGroup = group.createMeshGroup(mesh1d) + mesh_group_add_identifier_ranges(meshGroup, [[3, 3]]) + self.assertEqual(1, meshGroup.getSize()) + self.assertEqual(1, annotationGroup1.getDimension()) + identifier_ranges_string = identifier_ranges_to_string(mesh_group_to_identifier_ranges(meshGroup)) + self.assertEqual("3", identifier_ranges_string) + networkLayoutScaffoldPackage.updateUserAnnotationGroups() + + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + annotationGroups = scaffoldPackage.getAnnotationGroups() + self.assertEqual(3, len(annotationGroups)) + self.assertTrue(findAnnotationGroupByName(annotationGroups, "core") is not None) + self.assertTrue(findAnnotationGroupByName(annotationGroups, "shell") is not None) + + fieldmodule = region.getFieldmodule() + + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(544, mesh3d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(680, nodes.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + # check annotation group transferred to 3D tube + annotationGroup = findAnnotationGroupByName(annotationGroups, "segment 3") + self.assertTrue(annotationGroup is not None) + self.assertEqual("SEGMENT:3", annotationGroup.getId()) + self.assertEqual(128, annotationGroup.getMeshGroup(fieldmodule.findMeshByDimension(3)).getSize()) + + X_TOL = 1.0E-6 + + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.14453769545049167, -0.6221797157655231, -0.1000000000000000], X_TOL) + assertAlmostEqualList(self, maximums, [2.15, 0.6221795584199485, 0.1000000000000000], X_TOL) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + 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(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.10923734449979247, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 2.205681406661578, delta=X_TOL) + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_kidney.py b/tests/test_kidney.py new file mode 100644 index 00000000..35d5e067 --- /dev/null +++ b/tests/test_kidney.py @@ -0,0 +1,124 @@ +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.field import Field +from cmlibs.zinc.result import RESULT_OK + +from scaffoldmaker.annotation.annotationgroup import getAnnotationGroupForTerm +from scaffoldmaker.annotation.kidney_terms import get_kidney_term +from scaffoldmaker.meshtypes.meshtype_3d_kidney1 import MeshType_3d_kidney1 + + +from testutils import assertAlmostEqualList + + +class KidneyScaffoldTestCase(unittest.TestCase): + + def test_kidney(self): + """ + Test creation of renal capsule scaffold. + """ + scaffold = MeshType_3d_kidney1 + parameterSetNames = scaffold.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default", "Human 1"]) + options = scaffold.getDefaultOptions("Human 1") + + self.assertEqual(13, len(options)) + self.assertEqual(8, options["Number of elements around"]) + self.assertEqual(1, options["Number of elements through cortex"]) + self.assertEqual([0], options["Annotation numbers of elements around"]) + self.assertEqual(2.0, options["Target element density along longest segment"]) + self.assertEqual(2, options["Number of elements across core box minor"]) + self.assertEqual(1, options["Number of elements across core transition"]) + self.assertEqual([0], options["Annotation numbers of elements across core box minor"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold.generateMesh(region, options)[0] + self.assertEqual(52, len(annotationGroups)) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(136 * 2, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(436 * 2, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(478 * 2, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(181 * 2, 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-4 + assertAlmostEqualList(self, minimums, [-0.5137479110210048, -0.75, -0.2], tol) + assertAlmostEqualList(self, maximums, [0.5137479110210048, 0.75, 0.2], tol) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + 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(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(4) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + + self.assertAlmostEqual(volume, 0.27576984106019536, delta=tol) + self.assertAlmostEqual(surfaceArea, 2.9158676929090253, delta=tol) + + # check some annotation groups: + + expectedSizes3d = { + "renal medulla": (80 * 2, 0.08481313137381906), + "cortex of kidney": (52 * 2, 0.1770664248818641), + "kidney": (136 * 2, 0.2757746910245236) + } + for name in expectedSizes3d: + term = get_kidney_term(name) + annotationGroup = getAnnotationGroupForTerm(annotationGroups, term) + size = annotationGroup.getMeshGroup(mesh3d).getSize() + self.assertEqual(expectedSizes3d[name][0], size, name) + volumeMeshGroup = annotationGroup.getMeshGroup(mesh3d) + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, volumeMeshGroup) + volumeField.setNumbersOfPoints(4) + fieldcache = fieldmodule.createFieldcache() + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, expectedSizes3d[name][1], delta=tol) + + expectedSizes2d = { + "kidney capsule": (56 * 2, 2.9158676929090253), + "anterior surface of kidney": (28 * 2, 1.4579338464549303), + "posterior surface of kidney": (28 * 2, 1.4579338464540916) + } + for name in expectedSizes2d: + term = get_kidney_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()