From a7783ceb80475e61ce28cdbefd147798e173987c Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 29 May 2026 13:30:49 +1200 Subject: [PATCH 1/3] Add define grid field on mesh utility --- src/cmlibs/utils/zinc/finiteelement.py | 62 +++++++ tests/resources/square.exf | 81 ++++++++ tests/resources/two_lines.exf | 55 ++++++ tests/resources/warped_two_element_cube.exf | 2 +- tests/test_zinc_finiteelement.py | 196 ++++++++++++++++++++ 5 files changed, 395 insertions(+), 1 deletion(-) create mode 100644 tests/resources/square.exf create mode 100644 tests/resources/two_lines.exf create mode 100644 tests/test_zinc_finiteelement.py diff --git a/src/cmlibs/utils/zinc/finiteelement.py b/src/cmlibs/utils/zinc/finiteelement.py index e22b91b..4302324 100644 --- a/src/cmlibs/utils/zinc/finiteelement.py +++ b/src/cmlibs/utils/zinc/finiteelement.py @@ -660,6 +660,68 @@ def get_scalar_field_minimum_in_mesh(real_field, mesh=None): return minimum_element_id, minimum_element_value +def define_grid_field_on_mesh(mesh, field_name, number_of_cells, number_of_components=1, value_type=float): + """ + Define a grid field over all elements of a mesh. This is an embedded structured mesh of sub-cells over + the element space. Supports cube, square or line shape. Supports only trilinear grid cells. + This can add to an existing grid field provided the field is defined identically. + New field has zero values. + :param mesh: Mesh to define grid field on. + :param field_name: Name of grid field to define. + :param number_of_cells: List giving number of sub-cells in each element direction. If fewer values + than the mesh dimension are supplied, the last number of cells is used for remaining directions. + :param number_of_components: Number of components of grid field. + :param value_type: float or int. + """ + number_count = len(number_of_cells) + assert number_count > 0 + assert all((number >= 1) for number in number_of_cells) + assert number_of_components >= 1 + assert value_type in (float, int) + mesh_dimension = mesh.getDimension() + face_mesh_string = "" if (mesh_dimension == 1) else f", face mesh=mesh{mesh_dimension - 1}d" + number_in_xi = [] + number_in_xi_strings = [] + for d in range(mesh_dimension): + number_in_xi.append(number_of_cells[d if (d < number_count) else -1]) + number_in_xi_strings.append(f" #xi{d + 1}={number_in_xi[d]}") + number_in_xi_string = ",".join(number_in_xi_strings) + value_type_string = "integer" if (value_type is int) else "real" + basis_string = "*".join(["l.Lagrange"] * mesh_dimension) + components_string = "" + for c in range(number_of_components): + components_string += f" {c + 1}. {basis_string}, no modify, grid based.\n" + number_in_xi_string + "\n" + header_string = \ +f"""EX Version: 3 +Region: / +!#mesh mesh{mesh_dimension}d, dimension={mesh_dimension}{face_mesh_string}, nodeset=nodes +Define element template: element1 +Shape. Dimension={mesh_dimension}, {"*".join(["line"] * mesh_dimension)} +#Scale factor sets=0 +#Nodes=0 +#Fields=1 +1) {field_name}, field, {value_type_string}, #Components={number_of_components} +{components_string}Element template: element1 +""" + higher_number_in_xi = number_of_components + for d in range(1, mesh_dimension): + higher_number_in_xi *= (number_in_xi[d] + 1) + values_string = "Values:\n" + (" 0" * (number_in_xi[0] + 1) + "\n") * higher_number_in_xi + ex_strings = [header_string] + elementiterator = mesh.createElementiterator() + element = elementiterator.next() + while element.isValid(): + element_string = "Element: " + str(element.getIdentifier()) + "\n" + values_string + ex_strings.append(element_string) + element = elementiterator.next() + ex_string = "".join(ex_strings) + region = mesh.getFieldmodule().getRegion() + sir = region.createStreaminformationRegion() + sir.createStreamresourceMemoryBuffer(ex_string) + result = region.read(sir) + assert result == RESULT_OK + + createCubeElement = create_cube_element createSquareElement = create_square_element createLineElement = create_line_element diff --git a/tests/resources/square.exf b/tests/resources/square.exf new file mode 100644 index 0000000..f7921f0 --- /dev/null +++ b/tests/resources/square.exf @@ -0,0 +1,81 @@ +EX Version: 3 +Region: / +!#nodeset nodes +Define node template: node1 +Shape. Dimension=0 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +Node template: node1 +Node: 1 + 0.000000000000000e+00 + 0.000000000000000e+00 + 0.000000000000000e+00 +Node: 2 + 1.000000000000000e+00 + 0.000000000000000e+00 + 0.000000000000000e+00 +Node: 3 + 0.000000000000000e+00 + 1.000000000000000e+00 + 0.000000000000000e+00 +Node: 4 + 1.000000000000000e+00 + 1.000000000000000e+00 + 0.000000000000000e+00 +!#mesh mesh1d, dimension=1, nodeset=nodes +Define element template: element1 +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element template: element1 +Element: 1 +Element: 2 +Element: 3 +Element: 4 +!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes +Define element template: element2 +Shape. Dimension=2, line*line +#Scale factor sets=0 +#Nodes=4 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. l.Lagrange*l.Lagrange, no modify, standard node based. + #Nodes=4 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value + 3. #Values=1 + Value labels: value + 4. #Values=1 + Value labels: value + y. l.Lagrange*l.Lagrange, no modify, standard node based. + #Nodes=4 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value + 3. #Values=1 + Value labels: value + 4. #Values=1 + Value labels: value + z. l.Lagrange*l.Lagrange, no modify, standard node based. + #Nodes=4 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value + 3. #Values=1 + Value labels: value + 4. #Values=1 + Value labels: value +Element template: element2 +Element: 1 + Faces: + 1 2 3 4 + Nodes: + 1 2 3 4 diff --git a/tests/resources/two_lines.exf b/tests/resources/two_lines.exf new file mode 100644 index 0000000..5638d5d --- /dev/null +++ b/tests/resources/two_lines.exf @@ -0,0 +1,55 @@ +EX Version: 3 +Region: / +!#nodeset nodes +Define node template: node1 +Shape. Dimension=0 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +Node template: node1 +Node: 1 + 0.0 + 0.0 + 0.0 +Node: 2 + 1.0 + 0.0 + 0.0 +Node: 3 + 2.0 + 1.0 + 1.0 +!#mesh mesh1d, dimension=1, nodeset=nodes +Define element template: element1 +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=2 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. l.Lagrange, no modify, standard node based. + #Nodes=2 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value + y. l.Lagrange, no modify, standard node based. + #Nodes=2 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value + z. l.Lagrange, no modify, standard node based. + #Nodes=2 + 1. #Values=1 + Value labels: value + 2. #Values=1 + Value labels: value +Element template: element1 +Element: 1 + Nodes: + 1 2 +Element: 2 + Nodes: + 2 3 diff --git a/tests/resources/warped_two_element_cube.exf b/tests/resources/warped_two_element_cube.exf index e127038..91935e2 100644 --- a/tests/resources/warped_two_element_cube.exf +++ b/tests/resources/warped_two_element_cube.exf @@ -340,4 +340,4 @@ Element: 2 7 8 4 9 10 11 Nodes: 3 4 5 6 9 10 11 12 -Group name: .scene_selection + diff --git a/tests/test_zinc_finiteelement.py b/tests/test_zinc_finiteelement.py new file mode 100644 index 0000000..5bf7537 --- /dev/null +++ b/tests/test_zinc_finiteelement.py @@ -0,0 +1,196 @@ +import math +import unittest + +from cmlibs.utils.zinc.finiteelement import define_grid_field_on_mesh +from cmlibs.zinc.context import Context +from cmlibs.zinc.result import RESULT_OK +from utilities import assert_almost_equal_list, get_test_resource_name + + +class ZincFiniteElementTestCase(unittest.TestCase): + + def test_mesh1d_grid_vector(self): + exf_file = get_test_resource_name('two_lines.exf') + + context = Context("test") + region = context.createRegion() + result = region.readFile(exf_file) + self.assertTrue(result == RESULT_OK) + + fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(1) + + define_grid_field_on_mesh(mesh, "grid_coordinates", [4], 3, value_type=float) + + coordinates = fm.findFieldByName("coordinates") + grid_coordinates = fm.findFieldByName("grid_coordinates") + self.assertTrue(grid_coordinates.isValid()) + self.assertEqual(grid_coordinates.getNumberOfComponents(), 3) + + # Fieldassignment doesn't work for grid, so assign to all points + # fieldassignment = grid_coordinates.createFieldassignment(coordinates) + # result = fieldassignment.assign() + # self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi1 in (0.0, 0.25, 0.5, 0.75, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, xi1) + result, x = coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + result = grid_coordinates.assignReal(fieldcache, x) + self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi1 in (0.0, 0.25, 0.5, 0.75, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, xi1) + result, x = coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + result, grid_x = grid_coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + assert_almost_equal_list(self, x, grid_x, delta=1.0E-8) + + def test_mesh2d_grid_scalar(self): + exf_file = get_test_resource_name('square.exf') + + context = Context("test") + region = context.createRegion() + result = region.readFile(exf_file) + self.assertTrue(result == RESULT_OK) + + fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(2) + + define_grid_field_on_mesh(mesh, "grid_mag", [2, 4], 1, value_type=float) + + coordinates = fm.findFieldByName("coordinates") + mag_coordinates = fm.createFieldMagnitude(coordinates) + grid_mag = fm.findFieldByName("grid_mag") + self.assertTrue(grid_mag.isValid()) + self.assertEqual(grid_mag.getNumberOfComponents(), 1) + + # Fieldassignment doesn't work for grid, so assign to all points + # fieldassignment = grid_mag.createFieldassignment(mag_coordinates) + # result = fieldassignment.assign() + # self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi2 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + result = grid_mag.assignReal(fieldcache, mag_x) + self.assertTrue(result == RESULT_OK) + # new field cache so not recalling from cache + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi2 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + result, grid_mag_x = grid_mag.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + self.assertAlmostEqual(mag_x, grid_mag_x, delta=1.0E-8) + + def test_mesh3d_grid_scalar(self): + exf_file = get_test_resource_name('warped_two_element_cube.exf') + + context = Context("test") + region = context.createRegion() + result = region.readFile(exf_file) + self.assertTrue(result == RESULT_OK) + + fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(3) + + define_grid_field_on_mesh(mesh, "grid_mag", [2, 3, 4], 1, value_type=float) + + coordinates = fm.findFieldByName("coordinates") + mag_coordinates = fm.createFieldMagnitude(coordinates) + grid_mag = fm.findFieldByName("grid_mag") + self.assertTrue(grid_mag.isValid()) + self.assertEqual(grid_mag.getNumberOfComponents(), 1) + + # Fieldassignment doesn't work for grid, so assign to all points + # fieldassignment = grid_mag.createFieldassignment(mag_coordinates) + # result = fieldassignment.assign() + # self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + result = grid_mag.assignReal(fieldcache, mag_x) + self.assertTrue(result == RESULT_OK) + # new field cache so not recalling from cache + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + result, grid_mag_x = grid_mag.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + self.assertAlmostEqual(mag_x, grid_mag_x, delta=1.0E-8) + + def test_mesh3d_grid_vector(self): + exf_file = get_test_resource_name('warped_two_element_cube.exf') + + context = Context("test") + region = context.createRegion() + result = region.readFile(exf_file) + self.assertTrue(result == RESULT_OK) + + fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(3) + + define_grid_field_on_mesh(mesh, "grid_coordinates", [4], 3, value_type=float) + + coordinates = fm.findFieldByName("coordinates") + grid_coordinates = fm.findFieldByName("grid_coordinates") + self.assertTrue(grid_coordinates.isValid()) + self.assertEqual(grid_coordinates.getNumberOfComponents(), 3) + + # Fieldassignment doesn't work for grid, so assign to all points + # fieldassignment = grid_coordinates.createFieldassignment(coordinates) + # result = fieldassignment.assign() + # self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi1 in (0.0, 0.25, 0.5, 0.75, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, x = coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + result = grid_coordinates.assignReal(fieldcache, x) + self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi1 in (0.0, 0.25, 0.5, 0.75, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, x = coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + result, grid_x = grid_coordinates.evaluateReal(fieldcache, 3) + self.assertTrue(result == RESULT_OK) + assert_almost_equal_list(self, x, grid_x, delta=1.0E-8) + + +if __name__ == "__main__": + unittest.main() From 649e82f06955183263936c1be7cdc2b0ec649340 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 29 May 2026 13:41:15 +1200 Subject: [PATCH 2/3] Add integer grid test --- tests/test_zinc_finiteelement.py | 53 +++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/tests/test_zinc_finiteelement.py b/tests/test_zinc_finiteelement.py index 5bf7537..cd4b4b1 100644 --- a/tests/test_zinc_finiteelement.py +++ b/tests/test_zinc_finiteelement.py @@ -4,7 +4,7 @@ from cmlibs.utils.zinc.finiteelement import define_grid_field_on_mesh from cmlibs.zinc.context import Context from cmlibs.zinc.result import RESULT_OK -from utilities import assert_almost_equal_list, get_test_resource_name +from .utilities import assert_almost_equal_list, get_test_resource_name class ZincFiniteElementTestCase(unittest.TestCase): @@ -145,6 +145,57 @@ def test_mesh3d_grid_scalar(self): self.assertTrue(result == RESULT_OK) self.assertAlmostEqual(mag_x, grid_mag_x, delta=1.0E-8) + def test_mesh3d_grid_scalar_int(self): + exf_file = get_test_resource_name('warped_two_element_cube.exf') + + context = Context("test") + region = context.createRegion() + result = region.readFile(exf_file) + self.assertTrue(result == RESULT_OK) + + fm = region.getFieldmodule() + mesh = fm.findMeshByDimension(3) + + define_grid_field_on_mesh(mesh, "grid_mag", [2, 3, 4], 1, value_type=int) + + coordinates = fm.findFieldByName("coordinates") + mag_coordinates = fm.createFieldMagnitude(coordinates) + grid_mag = fm.findFieldByName("grid_mag") + self.assertTrue(grid_mag.isValid()) + self.assertEqual(grid_mag.getNumberOfComponents(), 1) + + # Fieldassignment doesn't work for grid, so assign to all points + # fieldassignment = grid_mag.createFieldassignment(mag_coordinates) + # result = fieldassignment.assign() + # self.assertTrue(result == RESULT_OK) + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + int_mag_x = round(mag_x * 1000.0) + result = grid_mag.assignReal(fieldcache, int_mag_x) + self.assertTrue(result == RESULT_OK) + # new field cache so not recalling from cache + fieldcache = fm.createFieldcache() + for element_identifier in (1, 2): + for xi3 in (0.0, 0.25, 0.5, 0.75, 1.0): + for xi2 in (0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0): + for xi1 in (0.0, 0.5, 1.0): + element = mesh.findElementByIdentifier(element_identifier) + fieldcache.setMeshLocation(element, [xi1, xi2, xi3]) + result, mag_x = mag_coordinates.evaluateReal(fieldcache, 1) + self.assertTrue(result == RESULT_OK) + int_mag_x = round(mag_x * 1000.0) + result, grid_mag_x = grid_mag.evaluateReal(fieldcache, 1) + int_grid_mag_x = int(grid_mag_x) + self.assertTrue(result == RESULT_OK) + self.assertEqual(int_mag_x, int_grid_mag_x) + def test_mesh3d_grid_vector(self): exf_file = get_test_resource_name('warped_two_element_cube.exf') From baba8e375c1951b99100c5e584b73c0208bf6d36 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 29 May 2026 13:42:37 +1200 Subject: [PATCH 3/3] Fix import --- tests/test_zinc_finiteelement.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_zinc_finiteelement.py b/tests/test_zinc_finiteelement.py index cd4b4b1..37c0fcc 100644 --- a/tests/test_zinc_finiteelement.py +++ b/tests/test_zinc_finiteelement.py @@ -4,7 +4,7 @@ from cmlibs.utils.zinc.finiteelement import define_grid_field_on_mesh from cmlibs.zinc.context import Context from cmlibs.zinc.result import RESULT_OK -from .utilities import assert_almost_equal_list, get_test_resource_name +from utilities import assert_almost_equal_list, get_test_resource_name class ZincFiniteElementTestCase(unittest.TestCase):