diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 083d0f67..67b13e47 100644 --- 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 index c26b54ce..c8f62142 100644 --- 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 new file mode 100644 index 00000000..04693be9 --- /dev/null +++ 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 new file mode 100644 index 00000000..35132cb5 --- /dev/null +++ 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() diff --git a/tests/test_heart.py b/tests/test_heart.py index 1c1cdfc4..59ca64db 100644 --- a/tests/test_heart.py +++ b/tests/test_heart.py @@ -5,11 +5,7 @@ from scaffoldmaker.annotation.annotationgroup import AnnotationGroup from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1 from scaffoldmaker.utils.zinc_utils import evaluateFieldRange - -def assertAlmostEqualList(testcase, actualList, expectedList, delta): - assert len(actualList) == len(expectedList) - for actual, expected in zip(actualList, expectedList): - testcase.assertAlmostEqual(actual, expected, delta=delta) +from tests.testutils import assertAlmostEqualList class HeartScaffoldTestCase(unittest.TestCase): diff --git a/tests/testutils.py b/tests/testutils.py new file mode 100644 index 00000000..ff8f89bf --- /dev/null +++ b/tests/testutils.py @@ -0,0 +1,8 @@ +''' +Utility function for tests. +''' + +def assertAlmostEqualList(testcase, actualList, expectedList, delta): + assert len(actualList) == len(expectedList) + for actual, expected in zip(actualList, expectedList): + testcase.assertAlmostEqual(actual, expected, delta=delta) \ No newline at end of file