Merge pull request #246 from rchristie/box_network

Box network

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",
 ]

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)

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

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

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]

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.

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).

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):
         """