Merge pull request #61 from mlin865/unittest

Tests for colon and colon segment meshtypes
ABI-Software · Nov 27, 2019 · 3e32f0e · 3e32f0e
2 parents 65fa7cb + 222a331
commit 3e32f0e
Show file tree

Hide file tree

Showing 6 changed files with 342 additions and 16 deletions.
diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py
@@ -138,17 +138,6 @@ class MeshType_3d_colon1(Scaffold_base):
                 [ [   32.0, -49.3,  24.4 ], [   -51.1, -50.1,  -1.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ],
                 [ [  -51.1, -45.3,  15.7 ], [   -38.0,  52.2, -16.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ],
                 [ [ -122.9, 124.2, -36.0 ], [   -21.3,  64.5, -18.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ] ] )
-            } ),
-        'Test Line' : ScaffoldPackage(MeshType_1d_path1, {
-            'scaffoldSettings' : {
-                'Coordinate dimensions' : 3,
-                'Length' : 1.0,
-                'Number of elements' : 1
-                },
-            'meshEdits' : zinc_utils.exnodeStringFromNodeValues(
-                [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [
-                [ [ -40.0, 10.0, 30.0 ], [ 50.0, 10.0, -30.0 ], [ 8.0, -40.0, 0.0 ], [0.0, 0.0, 0.0 ] ],
-                [ [  10.0, 20.0,  0.0 ], [ 50.0, 10.0, -30.0 ], [ 8.0, -40.0, 0.0 ], [0.0, 0.0, 0.0 ] ] ])
             } )
         }
 

diff --git a/src/scaffoldmaker/utils/zinc_utils.py b/src/scaffoldmaker/utils/zinc_utils.py
@@ -3,7 +3,9 @@
 '''
 
 from opencmiss.zinc.context import Context
+from opencmiss.zinc.element import MeshGroup
 from opencmiss.zinc.field import Field
+from opencmiss.zinc.fieldmodule import Fieldmodule
 from opencmiss.zinc.node import Node, Nodeset
 from opencmiss.zinc.result import RESULT_OK
 from scaffoldmaker.utils import interpolation as interp
@@ -416,3 +418,19 @@ def exnodeStringFromNodeValues(
         region.write(sir)
         result, exString = srm.getBuffer()
     return exString
+
+def createFaceMeshGroupExteriorOnFace(fieldmodule : Fieldmodule, elementFaceType) -> MeshGroup:
+    """
+    Returns mesh group for the exterior surface on the face described
+    by elementFaceType.
+    """
+    with ZincCacheChanges(fieldmodule):
+        isExterior = fieldmodule.createFieldIsExterior()
+        isOnFace = fieldmodule.createFieldIsOnFace(elementFaceType)
+        mesh2d = fieldmodule.findMeshByDimension(2)
+        faceElementGroup = fieldmodule.createFieldElementGroup(mesh2d)
+        faceMeshGroup = faceElementGroup.getMeshGroup()
+        faceMeshGroup.addElementsConditional(fieldmodule.createFieldAnd(isExterior, isOnFace))
+        del isExterior
+        del isOnFace
+    return faceMeshGroup
diff --git a/tests/test_colon.py b/tests/test_colon.py
@@ -0,0 +1,182 @@
+import copy
+import unittest
+from opencmiss.zinc.context import Context
+from opencmiss.zinc.element import Element
+from opencmiss.zinc.field import Field
+from opencmiss.zinc.node import Node
+from opencmiss.zinc.result import RESULT_OK
+from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion
+from scaffoldmaker.meshtypes.meshtype_3d_colon1 import MeshType_3d_colon1
+from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1
+from scaffoldmaker.scaffoldpackage import ScaffoldPackage
+from scaffoldmaker.utils import zinc_utils
+from tests.testutils import assertAlmostEqualList
+
+class ColonScaffoldTestCase(unittest.TestCase):
+
+    def test_colon1(self):
+        """
+        Test creation of colon scaffold.
+        """
+        parameterSetNames = MeshType_3d_colon1.getParameterSetNames()
+        self.assertEqual(parameterSetNames, ["Default", "Human 1", "Human 2", "Mouse 1", "Mouse 2", "Pig 1"])
+        centralPathDefaultScaffoldPackages = {
+            'Test line': ScaffoldPackage(MeshType_1d_path1, {
+                'scaffoldSettings': {
+                    'Coordinate dimensions': 3,
+                    'Length': 1.0,
+                    'Number of elements': 1
+                },
+                'meshEdits': zinc_utils.exnodeStringFromNodeValues(
+                    [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2,
+                     Node.VALUE_LABEL_D2_DS1DS2], [
+                        [[163.7, -25.2, 12.2], [-21.7, 50.1, -18.1], [0.0, 0.0, 5.0], [0.0, 0.0, 0.5]],
+                        [[117.2, 32.8, -2.6], [-64.3, 34.4, -3.9], [0.0, 0.0, 5.0], [0.0, 0.0, 0.5]]])
+            })
+        }
+        centralPathOption = centralPathDefaultScaffoldPackages['Test line']
+        segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName='Human 1')
+        options = {
+            'Central path': copy.deepcopy(centralPathOption),
+            'Segment profile': segmentProfileOption,
+            'Number of segments': 3,
+            'Proximal length': 25.0,
+            'Transverse length': 25.0,
+            'Distal length': 25.0,
+            'Proximal inner radius': 20.0,
+            'Proximal tenia coli width': 8.0,
+            'Proximal-transverse inner radius': 18.0,
+            'Proximal-transverse tenia coli width': 6.0,
+            'Transverse-distal inner radius': 16.0,
+            'Transverse-distal tenia coli width': 5.0,
+            'Distal inner radius': 15.0,
+            'Distal tenia coli width': 5.0,
+            'Refine': False,
+            'Refine number of elements around': 1,
+            'Refine number of elements along': 1,
+            'Refine number of elements through wall': 1
+        }
+        self.assertEqual(18, len(options))
+        centralPath = options['Central path']
+        segmentProfile = options.get("Segment profile")
+        segmentSettings = segmentProfile.getScaffoldSettings()
+        self.assertEqual(8, segmentSettings.get("Number of elements around haustrum"))
+        self.assertEqual(0.5, segmentSettings.get("Corner inner radius factor"))
+        self.assertEqual(0.5, segmentSettings.get("Haustrum inner radius factor"))
+        self.assertEqual(0.5, segmentSettings.get("Segment length end derivative factor"))
+        self.assertEqual(3, segmentSettings.get("Number of tenia coli"))
+        self.assertEqual(1.6, segmentSettings.get("Tenia coli thickness"))
+        self.assertEqual(3, options.get("Number of segments"))
+        self.assertEqual(25.0, options.get("Transverse length"))
+        self.assertEqual(20.0, options.get("Proximal inner radius"))
+        self.assertEqual(6.0, options.get("Proximal-transverse tenia coli width"))
+        self.assertEqual(16.0, options.get("Transverse-distal inner radius"))
+        self.assertEqual(5.0, options.get("Distal tenia coli width"))
+
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+
+        tmpRegion = region.createRegion()
+        centralPath.generate(tmpRegion)
+        cx = extractPathParametersFromRegion(tmpRegion)[0]
+        self.assertEqual(2, len(cx))
+        assertAlmostEqualList(self, cx[0], [ 163.7, -25.2, 12.2 ], 1.0E-6)
+        assertAlmostEqualList(self, cx[1], [ 117.2, 32.8, -2.6 ], 1.0E-6)
+        del tmpRegion
+
+        annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options)
+        self.assertEqual(3, len(annotationGroups))
+
+        fieldmodule = region.getFieldmodule()
+        self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
+        if annotationGroups is not None:
+            for annotationGroup in annotationGroups:
+                annotationGroup.addSubelements()
+        mesh3d = fieldmodule.findMeshByDimension(3)
+        self.assertEqual(432, mesh3d.getSize())
+        mesh2d = fieldmodule.findMeshByDimension(2)
+        self.assertEqual(1656, mesh2d.getSize())
+        mesh1d = fieldmodule.findMeshByDimension(1)
+        self.assertEqual(2043, mesh1d.getSize())
+        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        self.assertEqual(819, nodes.getSize())
+        datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        self.assertEqual(0, datapoints.getSize())
+
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        minimums, maximums = zinc_utils.evaluateFieldRange(coordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ 108.05453644074798, -36.659788201178515, -25.896225462626255 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [ 185.4128545433126, 48.1410866643588, 34.90780743659052 ], 1.0E-6)
+
+        flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement()
+        self.assertTrue(flatCoordinates.isValid())
+        minimums, maximums = zinc_utils.evaluateFieldRange(flatCoordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [ 186.72664370397405, 77.41890571041102, 3.2000000000000006 ], 1.0E-6)
+
+        textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement()
+        minimums, maximums = zinc_utils.evaluateFieldRange(textureCoordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [ 0.9812487204616481, 1.0, 2.0 ], 1.0E-6)
+
+        with zinc_utils.ZincCacheChanges(fieldmodule):
+            one = fieldmodule.createFieldConstant(1.0)
+            faceMeshGroup = zinc_utils.createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1)
+            surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, faceMeshGroup)
+            surfaceAreaField.setNumbersOfPoints(4)
+            volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d)
+            volumeField.setNumbersOfPoints(3)
+        fieldcache = fieldmodule.createFieldcache()
+        result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(surfaceArea, 14244.909601959167, delta=1.0E-6)
+        result, volume = volumeField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(volume, 24014.934316646588, delta=1.0E-6)
+
+    def test_mousecolon1(self):
+        """
+        Test creation of mouse colon scaffold.
+        """
+        options = MeshType_3d_colon1.getDefaultOptions("Mouse 2")
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+        annotationGroups = MeshType_3d_colon1.generateBaseMesh(region, options)
+        self.assertEqual(2, len(annotationGroups))
+
+        fieldmodule = region.getFieldmodule()
+        self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
+        mesh3d = fieldmodule.findMeshByDimension(3)
+        self.assertEqual(400, mesh3d.getSize())
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement()
+        self.assertTrue(flatCoordinates.isValid())
+        textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement()
+        self.assertTrue(textureCoordinates.isValid())
+
+        with zinc_utils.ZincCacheChanges(fieldmodule):
+            one = fieldmodule.createFieldConstant(1.0)
+            faceMeshGroup = zinc_utils.createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1)
+            flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup)
+            flatSurfaceAreaField.setNumbersOfPoints(4)
+            textureSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, faceMeshGroup)
+            textureSurfaceAreaField.setNumbersOfPoints(4)
+            textureVolumeField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, mesh3d)
+            textureVolumeField.setNumbersOfPoints(3)
+        fieldcache = fieldmodule.createFieldcache()
+        result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(flatSurfaceArea, 641.9768900558535, delta=1.0E-3)
+        result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6)
+        result, textureVolume = textureVolumeField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(textureVolume, 1.0, delta=1.0E-6)
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_colonsegment.py b/tests/test_colonsegment.py
@@ -0,0 +1,133 @@
+import unittest
+from opencmiss.zinc.context import Context
+from opencmiss.zinc.element import Element
+from opencmiss.zinc.field import Field
+from opencmiss.zinc.result import RESULT_OK
+from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1
+from scaffoldmaker.utils.zinc_utils import evaluateFieldRange, ZincCacheChanges, createFaceMeshGroupExteriorOnFace
+from tests.testutils import assertAlmostEqualList
+
+class ColonSegmentScaffoldTestCase(unittest.TestCase):
+
+    def test_humancolonsegment1(self):
+        """
+        Test creation of human colon segment scaffold.
+        """
+        parameterSetNames = MeshType_3d_colonsegment1.getParameterSetNames()
+        self.assertEqual(parameterSetNames, [ "Default", "Human 1", "Mouse 1", "Pig 1" ])
+        options = MeshType_3d_colonsegment1.getDefaultOptions("Human 1")
+        self.assertEqual(26, len(options))
+        self.assertEqual(2, options.get("Number of elements around tenia coli"))
+        self.assertEqual(4, options.get("Number of elements along segment"))
+        self.assertEqual(1, options.get("Number of elements through wall"))
+        self.assertEqual(43.5, options.get("Start inner radius"))
+        self.assertEqual(0.0, options.get("Start inner radius derivative"))
+        self.assertEqual(33.0, options.get("End inner radius"))
+        self.assertEqual(0.0, options.get("End inner radius derivative"))
+        self.assertEqual(3.0, options.get("Segment length mid derivative factor"))
+        self.assertEqual(50.0, options.get("Segment length"))
+        self.assertEqual(3, options.get("Number of tenia coli"))
+        self.assertEqual(10.0, options.get("Start tenia coli width"))
+        self.assertEqual(0.0, options.get("End tenia coli width derivative"))
+        self.assertEqual(1.6, options.get("Wall thickness"))
+
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+        annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options)
+        self.assertEqual(3, len(annotationGroups))
+
+        fieldmodule = region.getFieldmodule()
+        self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
+        mesh3d = fieldmodule.findMeshByDimension(3)
+        self.assertEqual(144, mesh3d.getSize())
+        mesh2d = fieldmodule.findMeshByDimension(2)
+        self.assertEqual(576, mesh2d.getSize())
+        mesh1d = fieldmodule.findMeshByDimension(1)
+        self.assertEqual(747, mesh1d.getSize())
+        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        self.assertEqual(315, nodes.getSize())
+        datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
+        self.assertEqual(0, datapoints.getSize())
+
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        minimums, maximums = evaluateFieldRange(coordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ -2.172286248499807e-15, -58.70795100094912, -55.299531343601124 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [ 50.0, 48.620086017242905, 55.29953084925164 ], 1.0E-6)
+
+        flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement()
+        self.assertTrue(flatCoordinates.isValid())
+        minimums, maximums = evaluateFieldRange(flatCoordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [397.26513766571264, 50.0, 3.2000000000000006], 1.0E-6)
+
+        textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement()
+        minimums, maximums = evaluateFieldRange(textureCoordinates, nodes)
+        assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [ 0.9887754554800083, 1.0, 2.0 ], 1.0E-6)
+
+        with ZincCacheChanges(fieldmodule):
+            one = fieldmodule.createFieldConstant(1.0)
+            faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1)
+            surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, faceMeshGroup)
+            surfaceAreaField.setNumbersOfPoints(4)
+            volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d)
+            volumeField.setNumbersOfPoints(3)
+        fieldcache = fieldmodule.createFieldcache()
+        result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(surfaceArea, 20791.077901965793, delta=1.0E-6)
+        result, volume = volumeField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(volume, 36020.03580218517, delta=1.0E-6)
+
+    def test_mousecolonsegment1(self):
+        """
+        Test creation of mouse colon segment scaffold.
+        """
+        options = MeshType_3d_colonsegment1.getDefaultOptions("Mouse 1")
+        context = Context("Test")
+        region = context.getDefaultRegion()
+        self.assertTrue(region.isValid())
+        annotationGroups = MeshType_3d_colonsegment1.generateBaseMesh(region, options)
+        self.assertEqual(2, len(annotationGroups))
+
+        fieldmodule = region.getFieldmodule()
+        self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
+        mesh3d = fieldmodule.findMeshByDimension(3)
+        self.assertEqual(40, mesh3d.getSize())
+        coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
+        self.assertTrue(coordinates.isValid())
+        flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement()
+        self.assertTrue(flatCoordinates.isValid())
+        textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement()
+        self.assertTrue(textureCoordinates.isValid())
+
+        with ZincCacheChanges(fieldmodule):
+            one = fieldmodule.createFieldConstant(1.0)
+            faceMeshGroup = createFaceMeshGroupExteriorOnFace(fieldmodule, Element.FACE_TYPE_XI3_1)
+            surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, faceMeshGroup)
+            surfaceAreaField.setNumbersOfPoints(4)
+            flatSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, flatCoordinates, faceMeshGroup)
+            flatSurfaceAreaField.setNumbersOfPoints(4)
+            textureSurfaceAreaField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, faceMeshGroup)
+            textureSurfaceAreaField.setNumbersOfPoints(4)
+            textureVolumeField = fieldmodule.createFieldMeshIntegral(one, textureCoordinates, mesh3d)
+            textureVolumeField.setNumbersOfPoints(3)
+        fieldcache = fieldmodule.createFieldcache()
+        result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(surfaceArea, 468.17062489996886, delta=1.0E-6)
+        result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(flatSurfaceArea, surfaceArea, delta=1.0E-3)
+        result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6)
+        result, textureVolume = textureVolumeField.evaluateReal(fieldcache, 1)
+        self.assertEqual(result, RESULT_OK)
+        self.assertAlmostEqual(textureVolume, 1.0, delta=1.0E-6)
+
+if __name__ == "__main__":
+    unittest.main()