diff --git a/setup.py b/setup.py index 53429ad0..9e221082 100644 --- a/setup.py +++ b/setup.py @@ -27,7 +27,7 @@ def readfile(filename, split=False): # minimal requirements listing "cmlibs.maths >= 0.3", "cmlibs.utils >= 0.6", - "cmlibs.zinc >= 4.0", + "cmlibs.zinc >= 4.1", "scipy", "numpy", ] diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py index fa958f21..eac2d920 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py @@ -143,7 +143,7 @@ def generateBaseMesh(cls, region, options): cd2[nodeIndexes[ln]].append(d2) cd13[nodeIndexes[ln]].append(mult(d1Unit, edgeAngle * tubeRadius)) # fix the one node out of order: - for d in [cd1[4], cd2[4]]: + for d in [cd1[4], cd2[4], cd13[4]]: d[0:2] = [d[1], d[0]] for n in range(8): node = nodes.findNodeByIdentifier(n + 1) diff --git a/src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py b/src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py index 7875a1c0..3d6e2315 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py @@ -81,7 +81,7 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non @classmethod def checkOptions(cls, options): if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"): - options["Network layout"] = cls.getOptionScaffoldPackage("Network layout") + options["Network layout"] = cls.getOptionScaffoldPackage("Network layout", MeshType_1d_network_layout1) elementsCountAround = options["Elements count around"] if options["Elements count around"] < 4: options["Elements count around"] = 4 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_boxnetwork1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_boxnetwork1.py index 381fdc29..b097a5c9 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_boxnetwork1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_boxnetwork1.py @@ -1,9 +1,6 @@ """ Generates a hermite x bilinear 3-D box network mesh from a 1-D network layout. """ - -import copy - from cmlibs.maths.vectorops import add, sub from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates from cmlibs.zinc.element import Element, Elementbasis @@ -13,7 +10,111 @@ 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.networkmesh import NetworkMesh +from scaffoldmaker.utils.eft_utils import remapEftLocalNodes, remapEftNodeValueLabelVersion, setEftScaleFactorIds +from scaffoldmaker.utils.interpolation import getCubicHermiteCurvesLength, interpolateSampleCubicHermite, \ + sampleCubicHermiteCurvesSmooth +from scaffoldmaker.utils.zinc_utils import get_nodeset_path_ordered_field_parameters +import math + + +class BoxSegmentData: + """ + Definition of box mesh for network segment. + """ + + def __init__(self, pathParameters): + """ + :param pathParameters: (sx[], sd1[], sd2[], sd12[], sd3[], sd13[]) + """ + self._pathParameters = pathParameters # original parameters from network layout + self._segmentLength = getCubicHermiteCurvesLength(pathParameters[0], pathParameters[1]) + self._sampledBoxCoordinates = None + self._sampledNodeIdentifiers = [] + self._annotationMeshGroups = [] + + def getPathParameters(self): + return self._pathParameters + + def getSegmentLength(self): + return self._segmentLength + + def getSampledBoxCoordinates(self): + return self._sampledBoxCoordinates + + def setSampledBoxCoordinates(self, sampledBoxCoordinates): + """ + :param sampledBoxCoordinates: (sx[], sd1[], sd2[], sd12[], sd3[], sd13[]) + """ + self._sampledBoxCoordinates = sampledBoxCoordinates + self._sampledNodeIdentifiers = [None] * len(self._sampledBoxCoordinates[0]) + + def getSampledElementsCountAlong(self): + """ + Must have previously called setSampledBoxCoordinates + :return: Number of elements along sampled box. + """ + return len(self._sampledBoxCoordinates[0]) - 1 + + def getSampledNodeIdentifier(self, nodeIndex): + """ + Get start node identifier for supplying to adjacent straight or junction. + :param nodeIndex: Index of sampled node from 0 to sampledElementsCountAlong, or -1 for last. + """ + return self._sampledNodeIdentifiers[nodeIndex] + + def setSampledNodeIdentifier(self, nodeIndex, nodeIdentifier): + """ + :param nodeIndex: Index of sampled node from 0 to sampledElementsCountAlong, or -1 for last. + :param nodeIdentifier: Identifier of node to set at nodeIndex. + """ + self._sampledNodeIdentifiers[nodeIndex] = nodeIdentifier + + def addAnnotationMeshGroup(self, annotationMeshGroup): + """ + Add an annotation mesh group for segment elements to be added to. + :param annotationMeshGroup: Mesh group to add. + """ + self._annotationMeshGroups.append(annotationMeshGroup) + + def getAnnotationMeshGroups(self): + return self._annotationMeshGroups + + +class BoxJunctionData: + """ + Describes junction between multiple segments, some in, some out. + """ + + def __init__(self, networkSegmentsIn: list, networkSegmentsOut: list, boxSegmentData): + """ + :param networkSegmentsIn: List of input segments. + :param networkSegmentsOut: List of output segments. + :param boxSegmentData: dict NetworkSegment -> BoxSegmentData. + """ + self._networkSegmentsIn = networkSegmentsIn + self._networkSegmentsOut = networkSegmentsOut + self._networkSegments = networkSegmentsIn + networkSegmentsOut + segmentsCount = len(self._networkSegments) + assert segmentsCount > 0 + self._boxData = [boxSegmentData[networkSegment] for networkSegment in self._networkSegments] + self._segmentsIn = [self._networkSegments[s] in self._networkSegmentsIn for s in range(segmentsCount)] + self._nodeIdentifier = None + + def getSegmentsIn(self): + return self._segmentsIn + + def getBoxData(self): + return self._boxData + + def getNodeIdentifier(self): + return self._nodeIdentifier + + def setNodeIdentifier(self, nodeIdentifier): + self._nodeIdentifier = nodeIdentifier + segmentsCount = len(self._networkSegments) + for s in range(segmentsCount): + nodeIndex = -1 if self._segmentsIn[s] else 0 + self._boxData[s].setSampledNodeIdentifier(nodeIndex, nodeIdentifier) class MeshType_3d_boxnetwork1(Scaffold_base): @@ -21,21 +122,27 @@ class MeshType_3d_boxnetwork1(Scaffold_base): Generates a hermite x bilinear 3-D box network mesh from a 1-D network layout. """ - @staticmethod - def getName(): + @classmethod + def getName(cls): return "3D Box Network 1" + @classmethod + def getParameterSetNames(cls): + return MeshType_1d_network_layout1.getParameterSetNames() + @classmethod def getDefaultOptions(cls, parameterSetName="Default"): options = { - "Network layout": ScaffoldPackage(MeshType_1d_network_layout1) + "Network layout": ScaffoldPackage(MeshType_1d_network_layout1, defaultParameterSetName=parameterSetName), + "Target element density along longest segment": 4.0 } return options - @staticmethod - def getOrderedOptionNames(): + @classmethod + def getOrderedOptionNames(cls): return [ - "Network layout" + "Network layout", + "Target element density along longest segment" ] @classmethod @@ -70,12 +177,14 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non @classmethod def checkOptions(cls, options): if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"): - options["Network layout"] = cls.getOptionScaffoldPackage("Network layout") + options["Network layout"] = cls.getOptionScaffoldPackage("Network layout", MeshType_1d_network_layout1) + if options["Target element density along longest segment"] < 1.0: + options["Target element density along longest segment"] = 1.0 dependentChanges = False return dependentChanges - @staticmethod - def generateBaseMesh(region, options): + @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. @@ -83,6 +192,7 @@ def generateBaseMesh(region, options): :return: list of AnnotationGroup, None """ networkLayout = options["Network layout"] + targetElementDensityAlongLongestSegment = options["Target element density along longest segment"] layoutRegion = region.createRegion() layoutFieldmodule = layoutRegion.getFieldmodule() @@ -91,7 +201,6 @@ def generateBaseMesh(region, options): networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters layoutAnnotationGroups = networkLayout.getAnnotationGroups() layoutCoordinates = findOrCreateFieldCoordinates(layoutFieldmodule) - layoutFieldcache = layoutFieldmodule.createFieldcache() networkMesh = networkLayout.getConstructionObject() @@ -101,95 +210,168 @@ def generateBaseMesh(region, options): nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) nodeIdentifier = 1 - nodetemplate = nodes.createNodetemplate() - nodetemplate.defineField(coordinates) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplates = {} # map from derivativeVersionsCount to nodetemplate mesh = fieldmodule.findMeshByDimension(3) elementIdentifier = 1 hermiteBilinearBasis = fieldmodule.createElementbasis(3, Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE) hermiteBilinearBasis.setFunctionType(1, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) - eft = mesh.createElementfieldtemplate(hermiteBilinearBasis) - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate.defineField(coordinates, -1, eft) + elementtemplates = {} # map from (startVersion, endVersion) to (elementtemplate, eft) # make box annotation groups from network layout annotations annotationGroups = [] - layoutBoxMeshGroups = {} # map from group name + layoutAnnotationMeshGroupMap = [] # List of tuples of (layout annotation mesh group, box mesh group) for layoutAnnotationGroup in layoutAnnotationGroups: if layoutAnnotationGroup.getDimension() == 1: annotationGroup = AnnotationGroup(region, layoutAnnotationGroup.getTerm()) annotationGroups.append(annotationGroup) - layoutBoxMeshGroups[layoutAnnotationGroup.getName()] = \ - (layoutAnnotationGroup.getMeshGroup(layoutMesh), annotationGroup.getMeshGroup(mesh)) + layoutAnnotationMeshGroupMap.append( + (layoutAnnotationGroup.getMeshGroup(layoutMesh), annotationGroup.getMeshGroup(mesh))) + + valueLabels = [ + Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3] networkSegments = networkMesh.getNetworkSegments() - slices = {} # map from network layout node identifier to list of 4 box node identifiers on slice + boxSegmentData = {} # map from NetworkSegment to BoxSegmentData + longestSegmentLength = 0.0 for networkSegment in networkSegments: + pathParameters = get_nodeset_path_ordered_field_parameters( + layoutNodes, layoutCoordinates, valueLabels, + networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions()) + boxSegmentData[networkSegment] = boxData = BoxSegmentData(pathParameters) + segmentLength = boxData.getSegmentLength() + if segmentLength > longestSegmentLength: + longestSegmentLength = segmentLength + for layoutAnnotationMeshGroup, annotationMeshGroup in layoutAnnotationMeshGroupMap: + if networkSegment.hasLayoutElementsInMeshGroup(layoutAnnotationMeshGroup): + boxData.addAnnotationMeshGroup(annotationMeshGroup) + if longestSegmentLength > 0.0: + targetElementLength = longestSegmentLength / targetElementDensityAlongLongestSegment + else: + targetElementLength = 1.0 + + # map from NetworkNodes to BoxJunctionData + boxNodeJunctionDataMap = {} + for networkSegment in networkSegments: + boxData = boxSegmentData[networkSegment] + segmentNodes = networkSegment.getNetworkNodes() + startSegmentNode = segmentNodes[0] + startBoxJunctionData = boxNodeJunctionDataMap.get(startSegmentNode) + if not startBoxJunctionData: + startInSegments = startSegmentNode.getInSegments() + startOutSegments = startSegmentNode.getOutSegments() + startBoxJunctionData = BoxJunctionData(startInSegments, startOutSegments, boxSegmentData) + boxNodeJunctionDataMap[startSegmentNode] = startBoxJunctionData + endSegmentNode = segmentNodes[-1] + endBoxJunctionData = boxNodeJunctionDataMap.get(endSegmentNode) + if not endBoxJunctionData: + endInSegments = endSegmentNode.getInSegments() + endOutSegments = endSegmentNode.getOutSegments() + endBoxJunctionData = BoxJunctionData(endInSegments, endOutSegments, boxSegmentData) + boxNodeJunctionDataMap[endSegmentNode] = endBoxJunctionData + segmentLength = boxData.getSegmentLength() + elementsCountAlong = max(1, math.ceil(segmentLength / targetElementLength)) + loop = (len(startSegmentNode.getInSegments()) == 1) and \ + (startSegmentNode.getInSegments()[0] is networkSegment) and \ + (networkSegment.getNodeVersions()[0] == networkSegment.getNodeVersions()[-1]) + if (elementsCountAlong < 2) and loop: + elementsCountAlong = 2 # at least 2 segments around loop + pathParameters = boxData.getPathParameters() + sx, sd1, pe, pxi, psf = sampleCubicHermiteCurvesSmooth( + pathParameters[0], pathParameters[1], elementsCountAlong) + sd2, sd12 = interpolateSampleCubicHermite(pathParameters[2], pathParameters[3], pe, pxi, psf) + sd3, sd13 = interpolateSampleCubicHermite(pathParameters[4], pathParameters[5], pe, pxi, psf) + boxData.setSampledBoxCoordinates((sx, sd1, sd2, sd12, sd3, sd13)) + + lastNodeIdentifier = None + for networkSegment in networkSegments: + boxData = boxSegmentData[networkSegment] segmentNodes = networkSegment.getNetworkNodes() - # segmentVersions = networkSegment.getNodeVersions() - segmentElementIdentifiers = networkSegment.getElementIdentifiers() - segmentSlices = [] - segmentNodeCount = len(segmentNodes) - lastSlice = None - for n in range(segmentNodeCount): - segmentNode = segmentNodes[n] - layoutNodeIdentifier = segmentNode.getNodeIdentifier() - slice = slices.get(layoutNodeIdentifier) - if slice: - segmentSlices.append(slice) - lastSlice = slice - continue - layoutNode = layoutNodes.findNodeByIdentifier(layoutNodeIdentifier) - layoutFieldcache.setNode(layoutNode) - # currently only supports node version 1 - _, lx = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_VALUE, 1, 3) - _, ld1 = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) - _, ld2 = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) - _, ld12 = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, 3) - _, ld3 = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) - _, ld13 = layoutCoordinates.getNodeParameters(layoutFieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, 3) - - slice = [] - for n3 in range(2): - for n2 in range(2): - node = nodes.createNode(nodeIdentifier, nodetemplate) - fieldcache.setNode(node) - x = lx - d1 = ld1 - if n2 == 0: - x = sub(x, ld2) - d1 = sub(d1, ld12) - else: - x = add(x, ld2) - d1 = add(d1, ld12) - if n3 == 0: - x = sub(x, ld3) - d1 = sub(d1, ld13) - else: - x = add(x, ld3) - d1 = add(d1, ld13) - coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, x) - coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) - slice.append(nodeIdentifier) - nodeIdentifier += 1 - slices[layoutNodeIdentifier] = slice - - if lastSlice: - element = mesh.createElement(elementIdentifier, elementtemplate); - nids = [lastSlice[0], slice[0], - lastSlice[1], slice[1], - lastSlice[2], slice[2], - lastSlice[3], slice[3]] + layoutNodeVersions = networkSegment.getNodeVersions() + sx, sd1, sd2, sd12, sd3, sd13 = boxData.getSampledBoxCoordinates() + elementsCountAlong = boxData.getSampledElementsCountAlong() + annotationMeshGroups = boxData.getAnnotationMeshGroups() + for n in range(elementsCountAlong + 1): + thisNodeIdentifier = None + versionsCount = 1 + version = 1 + boxJunctionData = None + if n in [0, elementsCountAlong]: + nLayout = 0 if (n == 0) else -1 + segmentNode = segmentNodes[nLayout] + boxJunctionData = boxNodeJunctionDataMap[segmentNode] + thisNodeIdentifier = boxJunctionData.getNodeIdentifier() + version = layoutNodeVersions[nLayout] + versionsCount = segmentNode.getVersionsCount() + if thisNodeIdentifier is None: + nodetemplate = nodetemplates.get(versionsCount) + if not nodetemplate: + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + for valueLabel in valueLabels[1:]: + nodetemplate.setValueNumberOfVersions(coordinates, -1, valueLabel, versionsCount) + nodetemplates[versionsCount] = nodetemplate + node = nodes.createNode(nodeIdentifier, nodetemplate) + if boxJunctionData: + boxJunctionData.setNodeIdentifier(nodeIdentifier) + thisNodeIdentifier = nodeIdentifier + nodeIdentifier += 1 + else: + node = nodes.findNodeByIdentifier(thisNodeIdentifier) + # note following will set shared versions of coordinates or derivatives multiple times + # future: average derivatives from all adjoining segments + fieldcache.setNode(node) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, sx[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, version, sd1[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS2, version, sd2[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, sd12[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS3, version, sd3[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, sd13[n]) + + if n > 0: + startVersion = layoutNodeVersions[0] if (n == 1) else 1 + templateKey = (startVersion, version) + elementtemplate_and_eft = elementtemplates.get(templateKey) + if elementtemplate_and_eft: + elementtemplate, eft = elementtemplate_and_eft + else: + eft = mesh.createElementfieldtemplate(hermiteBilinearBasis) + setEftScaleFactorIds(eft, [1], []) + ln = 1 + for n3 in range(2): + s3 = [1] if (n3 == 0) else [] + for n2 in range(2): + s2 = [1] if (n2 == 0) else [] + for n1 in range(2): + v = startVersion if (n1 == 0) else version + remapEftNodeValueLabelVersion( + eft, [ln], Node.VALUE_LABEL_VALUE, + [(Node.VALUE_LABEL_VALUE, 1, []), + (Node.VALUE_LABEL_D_DS2, v, s2), + (Node.VALUE_LABEL_D_DS3, v, s3)]) + remapEftNodeValueLabelVersion( + eft, [ln], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, v, []), + (Node.VALUE_LABEL_D2_DS1DS2, v, s2), + (Node.VALUE_LABEL_D2_DS1DS3, v, s3)]) + ln += 1 + ln_map = [1, 2, 1, 2, 1, 2, 1, 2] + remapEftLocalNodes(eft, 2, ln_map) + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + elementtemplates[templateKey] = (elementtemplate, eft) + + element = mesh.createElement(elementIdentifier, elementtemplate) + nids = [lastNodeIdentifier, thisNodeIdentifier] element.setNodesByIdentifier(eft, nids) - layoutElementIdentifier = segmentElementIdentifiers[n - 1] - layoutElement = layoutMesh.findElementByIdentifier(layoutElementIdentifier) - for layoutBoxMeshGroup in layoutBoxMeshGroups.values(): - if layoutBoxMeshGroup[0].containsElement(layoutElement): - layoutBoxMeshGroup[1].addElement(element) + element.setScaleFactors(eft, [-1.0]) + for annotationMeshGroup in annotationMeshGroups: + annotationMeshGroup.addElement(element) elementIdentifier += 1 - lastSlice = slice + lastNodeIdentifier = thisNodeIdentifier return annotationGroups, None diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py index 0a320bc5..c2a2dcb6 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_tubenetwork1.py @@ -87,7 +87,7 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non @classmethod def checkOptions(cls, options): if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"): - options["Network layout"] = cls.getOptionScaffoldPackage("Network layout") + options["Network layout"] = cls.getOptionScaffoldPackage("Network layout", MeshType_1d_network_layout1) elementsCountAround = options["Elements count around"] if options["Elements count around"] < 4: options["Elements count around"] = 4 diff --git a/src/scaffoldmaker/utils/bifurcation.py b/src/scaffoldmaker/utils/bifurcation.py index 637b3698..84752e64 100644 --- a/src/scaffoldmaker/utils/bifurcation.py +++ b/src/scaffoldmaker/utils/bifurcation.py @@ -1683,7 +1683,7 @@ def generateTubeBifurcationTree(networkMesh: NetworkMesh, region, coordinates, n mesh = fieldmodule.findMeshByDimension(dimension) fieldcache = fieldmodule.createFieldcache() - # make 2D annotation groups from 1D network layout annotation groups + # make tube mesh annotation groups from 1D network layout annotation groups annotationGroups = [] layoutAnnotationMeshGroupMap = [] # List of tuples of layout annotation mesh group to final mesh group for layoutAnnotationGroup in layoutAnnotationGroups: @@ -1764,7 +1764,7 @@ def generateTubeBifurcationTree(networkMesh: NetworkMesh, region, coordinates, n if not startTubeBifurcationData: startInSegments = startSegmentNode.getInSegments() startOutSegments = startSegmentNode.getOutSegments() - if ((len(startInSegments) + len(startOutSegments)) == 3): + if (len(startInSegments) + len(startOutSegments)) == 3: # print("create start", networkSegment, startSegmentNode) startTubeBifurcationData = TubeBifurcationData(startInSegments, startOutSegments, segmentTubeData) nodeTubeBifurcationData[startSegmentNode] = startTubeBifurcationData @@ -1774,11 +1774,10 @@ def generateTubeBifurcationTree(networkMesh: NetworkMesh, region, coordinates, n endSegmentNode = segmentNodes[-1] endTubeBifurcationData = nodeTubeBifurcationData.get(endSegmentNode) endSurface = None - createEndBifurcationData = not endTubeBifurcationData - if createEndBifurcationData: + if not endTubeBifurcationData: endInSegments = endSegmentNode.getInSegments() endOutSegments = endSegmentNode.getOutSegments() - if ((len(endInSegments) + len(endOutSegments)) == 3): + if (len(endInSegments) + len(endOutSegments)) == 3: # print("create end", networkSegment, endSegmentNode) endTubeBifurcationData = TubeBifurcationData(endInSegments, endOutSegments, segmentTubeData) nodeTubeBifurcationData[endSegmentNode] = endTubeBifurcationData @@ -1802,9 +1801,9 @@ def generateTubeBifurcationTree(networkMesh: NetworkMesh, region, coordinates, n if (elementsCountAlong == 1) and startTubeBifurcationData and endTubeBifurcationData: # at least 2 segments if bifurcating at both ends elementsCountAlong = 2 - elif (elementsCountAlong < 3) and loop: - # at least 3 segments around loop; 2 should work, but zinc currently makes incorrect faces - elementsCountAlong = 3 + elif (elementsCountAlong < 2) and loop: + # at least 2 segments around loop + elementsCountAlong = 2 else: # must match count from outer surface! outerTubeData = outerSegmentTubeData[networkSegment] diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index be47b5ac..05856d5c 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -74,7 +74,8 @@ def remapEftNodeValueLabel(eft, localNodeIndexes, fromValueLabel, expressionTerm for f in range(1, functionCount + 1): if eft.getFunctionNumberOfTerms(f) == 1: localNodeIndex = eft.getTermLocalNodeIndex(f, 1) - if (localNodeIndex in localNodeIndexes) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and (not getEftTermScaling(eft, f, 1)): + if ((localNodeIndex in localNodeIndexes) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and + (not getEftTermScaling(eft, f, 1))): termCount = len(expressionTerms) eft.setFunctionNumberOfTerms(f, termCount) version = eft.getTermNodeVersion(f, 1) @@ -85,6 +86,30 @@ def remapEftNodeValueLabel(eft, localNodeIndexes, fromValueLabel, expressionTerm eft.setTermScaling(f, t, expressionTerm[1]) +def remapEftNodeValueLabelVersion(eft, localNodeIndexes, fromValueLabel, expressionTerms): + ''' + Remap all uses of the given valueLabels to the expressionTerms. + Note: Assumes valueLabel is currently single term and unscaled! + :param localNodeIndexes: List of local node indexes >= 1 to remap at. + :param fromValueLabel: Node value label to be remapped. + :param expressionTerms: List of (valueLabel, version, scaleFactorIndexesList) to remap to. + e.g. [ (Node.VALUE_LABEL_D_DS2, 1, []), (Node.VALUE_LABEL_D_DS3, 2, [5, 6]) ] + ''' + functionCount = eft.getNumberOfFunctions() + for f in range(1, functionCount + 1): + if eft.getFunctionNumberOfTerms(f) == 1: + localNodeIndex = eft.getTermLocalNodeIndex(f, 1) + if ((localNodeIndex in localNodeIndexes) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and + (not getEftTermScaling(eft, f, 1))): + termCount = len(expressionTerms) + eft.setFunctionNumberOfTerms(f, termCount) + for t in range(1, termCount + 1): + expressionTerm = expressionTerms[t - 1] + eft.setTermNodeParameter(f, t, localNodeIndex, expressionTerm[0], expressionTerm[1]) + if expressionTerm[2]: + eft.setTermScaling(f, t, expressionTerm[2]) + + def remapEftNodeValueLabelsVersion(eft, localNodeIndexes, valueLabels, version): ''' Remap all uses of the given valueLabels to use the version. diff --git a/src/scaffoldmaker/utils/interpolation.py b/src/scaffoldmaker/utils/interpolation.py index 967426dd..724f6c67 100644 --- a/src/scaffoldmaker/utils/interpolation.py +++ b/src/scaffoldmaker/utils/interpolation.py @@ -445,7 +445,7 @@ def sampleCubicHermiteCurvesSmooth(nx, nd1, elementsCountOut, derivatives appropriate for elementsCountOut. If unspecified these are calculated from the other end or set to be equal for even spaced elements. 0.0 is a valid derivative magnitude. :param startLocation: Optional tuple of 'in' (element, xi) to start curve at. - :param endLocation: Optional tuple of 'in' (element, xi) to end curve at. + :param endLocation: Optional tuple of 'out' (element, xi) to end curve at. :return: px[], pd1[], pe[], pxi[], psf[], where pe[] and pxi[] are lists of element indices and xi locations in the 'in' elements to pass to partner interpolateSample functions. psf[] is a list of scale factors for converting derivatives from old to new xi coordinates: dxi(old)/dxi(new). diff --git a/tests/test_network.py b/tests/test_network.py index 61d40353..9ebb8683 100644 --- a/tests/test_network.py +++ b/tests/test_network.py @@ -1,3 +1,4 @@ +import math import unittest from cmlibs.utils.zinc.finiteelement import evaluateFieldNodesetRange @@ -11,6 +12,7 @@ from cmlibs.zinc.result import RESULT_OK from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1 from scaffoldmaker.meshtypes.meshtype_2d_tubenetwork1 import MeshType_2d_tubenetwork1 +from scaffoldmaker.meshtypes.meshtype_3d_boxnetwork1 import MeshType_3d_boxnetwork1 from scaffoldmaker.meshtypes.meshtype_3d_tubenetwork1 import MeshType_3d_tubenetwork1 from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.utils.zinc_utils import get_nodeset_path_ordered_field_parameters @@ -183,8 +185,8 @@ def test_2d_tube_network_sphere_cube(self): X_TOL = 1.0E-6 minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) - assertAlmostEqualList(self, minimums, [-0.5663822833834603, -0.5836195039860382, -0.5981791955829937], X_TOL) - assertAlmostEqualList(self, maximums, [0.5664184183783737, 0.5965021612010702, 0.5981698817825564], X_TOL) + assertAlmostEqualList(self, minimums, [-0.5663822833834603, -0.5965021612010701, -0.598179822484941], X_TOL) + assertAlmostEqualList(self, maximums, [0.5663822151894911, 0.5965021612010702, 0.5981798465355403], X_TOL) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -193,7 +195,7 @@ def test_2d_tube_network_sphere_cube(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 4.033577495198774, delta=X_TOL) + self.assertAlmostEqual(surfaceArea, 4.057905325323945, delta=X_TOL) def test_3d_tube_network_sphere_cube(self): """ @@ -253,8 +255,8 @@ def test_3d_tube_network_sphere_cube(self): X_TOL = 1.0E-6 minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) - assertAlmostEqualList(self, minimums, [-0.5663822833834602, -0.5836195039860383, -0.5981791955829938], X_TOL) - assertAlmostEqualList(self, maximums, [0.5664184183783736, 0.5965021612010702, 0.5981698817825564], X_TOL) + assertAlmostEqualList(self, minimums, [-0.5663822833834603, -0.5965021612010702, -0.598179822484941], X_TOL) + assertAlmostEqualList(self, maximums, [0.5663822151894911, 0.5965021612010702, 0.5981798465355403], X_TOL) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -270,19 +272,73 @@ def test_3d_tube_network_sphere_cube(self): volumeField.setNumbersOfPoints(4) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 0.07344892961686832, delta=X_TOL) + self.assertAlmostEqual(volume, 0.074451650669961, delta=X_TOL) outerSurfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_1, coordinates, mesh2d) outerSurfaceAreaField.setNumbersOfPoints(4) result, outerSurfaceArea = outerSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(outerSurfaceArea, 4.0335774951987755, delta=X_TOL) + self.assertAlmostEqual(outerSurfaceArea, 4.057905325323947, delta=X_TOL) innerSurfaceAreaField = fieldmodule.createFieldMeshIntegral(isExteriorXi3_0, coordinates, mesh2d) innerSurfaceAreaField.setNumbersOfPoints(4) result, innerSurfaceArea = innerSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(innerSurfaceArea, 3.3276607694171063, delta=X_TOL) + self.assertAlmostEqual(innerSurfaceArea, 3.347440907189292, delta=X_TOL) + + def test_3d_box_network_bifurcation(self): + """ + Test 3-D box network bifurcation is generated correctly. + """ + scaffoldPackage = ScaffoldPackage(MeshType_3d_boxnetwork1, defaultParameterSetName="Bifurcation") + settings = scaffoldPackage.getScaffoldSettings() + networkLayoutScaffoldPackage = settings["Network layout"] + networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings() + self.assertEqual(2, len(settings)) + self.assertEqual(4.0, settings["Target element density along longest segment"]) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + scaffoldPackage.generate(region) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(12, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(63, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(108, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(13, nodes.getSize()) + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + + X_TOL = 1.0E-6 + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [0.0, -0.5, 0.0], X_TOL) + assertAlmostEqualList(self, maximums, [2.0, 0.5, 0.0], X_TOL) + + L2 = math.sqrt(1.25) + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + isExterior = fieldmodule.createFieldIsExterior() + fieldcache = fieldmodule.createFieldcache() + + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(1) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + expectedVolume = 0.2 * 0.2 * (1.0 + 2 * L2) + self.assertAlmostEqual(volume, expectedVolume, delta=X_TOL) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExterior, coordinates, mesh2d) + surfaceAreaField.setNumbersOfPoints(1) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + expectedSurfaceArea = 6 * 0.2 * 0.2 + 4 * 0.2 * (1.0 + 2 * L2) + self.assertAlmostEqual(surfaceArea, expectedSurfaceArea, delta=X_TOL) def test_3d_tube_network_loop(self): """