diff --git a/src/cmlibs/utils/zinc/field.py b/src/cmlibs/utils/zinc/field.py index 360bf35..483b328 100644 --- a/src/cmlibs/utils/zinc/field.py +++ b/src/cmlibs/utils/zinc/field.py @@ -2,6 +2,7 @@ Utilities for creating and working with Zinc Fields. """ from cmlibs.utils.zinc.general import ChangeManager +from cmlibs.utils.zinc.regionbufferio import region_read_from_buffer, region_write_to_buffer from cmlibs.zinc.element import Mesh from cmlibs.zinc.field import Field, FieldFiniteElement, FieldGroup, \ FieldStoredMeshLocation @@ -9,8 +10,6 @@ from cmlibs.zinc.node import Node from cmlibs.zinc.result import RESULT_OK, RESULT_WARNING_PART_DONE -from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer - def field_is_managed_coordinates(field_in: Field): """ @@ -382,12 +381,12 @@ def create_field_finite_element_clone(source_field: Field, name: str, managed=Fa # Zinc needs a function to do this efficiently; currently serialise to string, replace field name and reload! source_name = source_field.getName() region = fieldmodule.getRegion() - buffer = write_to_buffer(region, field_names=[source_name]) + buffer = region_write_to_buffer(region, field_names=[source_name]) # small risk of modifying other text here: source_bytes = bytes(") " + source_name + ",", "utf-8") target_bytes = bytes(") " + name + ",", "utf-8") buffer = buffer.replace(source_bytes, target_bytes) - result = read_from_buffer(region, buffer) + result = region_read_from_buffer(region, buffer) assert result == RESULT_OK # note currently must have called endChange before field can be found field = fieldmodule.findFieldByName(name).castFiniteElement() diff --git a/src/cmlibs/utils/zinc/finiteelement.py b/src/cmlibs/utils/zinc/finiteelement.py index f1eb591..c275978 100644 --- a/src/cmlibs/utils/zinc/finiteelement.py +++ b/src/cmlibs/utils/zinc/finiteelement.py @@ -4,6 +4,9 @@ import math from cmlibs.maths import vectorops +from cmlibs.utils.zinc.field import ( + find_or_create_field_coordinates, find_or_create_field_finite_element, find_or_create_field_group, + find_or_create_field_stored_string) from cmlibs.utils.zinc.general import ChangeManager from cmlibs.zinc.element import Element, Elementbasis, Elementfieldtemplate, Mesh from cmlibs.zinc.field import Field, FieldFindMeshLocation @@ -761,6 +764,118 @@ def define_grid_field_on_mesh(mesh, field_name, number_of_cells, number_of_compo assert result == RESULT_OK +def create_hermite_curve_elements(mesh, nx, nd1, loop=False, start_node_identifier=None, start_element_identifier=None, + coordinate_field_name="coordinates"): + """ + Generate a continuous 1-D mesh interpolated with cubic Hermite basis from the supplied coordinates and derivatives. + :param mesh: 1-D mesh or mesh group to define 1-D elements in. Nodes are added to the associated nodeset or + nodeset group. + :param nx: Coordinates along curve. + :param nd1: Derivatives along curve. + :param loop: True if curve loops back to first point, False if not. + :param start_node_identifier: Optional first node identifier to use. + :param start_element_identifier: Optional first 1D element identifier to use. + :param coordinate_field_name: Optional name of coordinate field to define, if omitted use "coordinates". + :return: next node identifier, next 1D element identifier + """ + assert mesh.getDimension() == 1 + fieldmodule = mesh.getFieldmodule() + # future, when supported: nodes = mesh.getMasterNodeset(), or mesh.getOrCreateNodeset() to replace following + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + mesh_group = mesh.castGroup() + if mesh_group.isValid(): + group = mesh_group.getFieldGroup() + nodes = group.getOrCreateNodesetGroup(nodes) + else: + group = None + with ChangeManager(fieldmodule): + coordinates = find_or_create_field_coordinates(fieldmodule, name=coordinate_field_name) + node_identifier = start_node_identifier if start_node_identifier is not None else \ + max(get_maximum_node_identifier(nodes), 0) + 1 + + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + + element_identifier = start_element_identifier if start_element_identifier is not None else \ + max(get_maximum_element_identifier(mesh), 0) + 1 + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_LINE) + cubicHermiteBasis = fieldmodule.createElementbasis( + 1, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + eft = mesh.createElementfieldtemplate(cubicHermiteBasis) + elementtemplate.defineField(coordinates, -1, eft) + + fieldcache = fieldmodule.createFieldcache() + nids = [] + nCount = len(nx) + for n in range(nCount): + node = nodes.createNode(node_identifier, nodetemplate) + nids.append(node_identifier) + fieldcache.setNode(node) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, nx[n]) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, nd1[n]) + node_identifier += 1 + eCount = nCount if loop else nCount - 1 + if loop: + nids.append(nids[0]) + for e in range(eCount): + element = mesh.createElement(element_identifier, elementtemplate) + element.setNodesByIdentifier(eft, nids[e:e + 2]) + element_identifier += 1 + + return node_identifier, element_identifier + + +def create_node_points(nodeset, px, start_node_identifier=None, + coordinate_field_name="coordinates", field_names_and_values=[]): + """ + Create a set of nodes or datapoints. + :param nodeset: Nodeset or nodeset group to create nodes in. + :param px: Coordinates of points. + :param start_data_identifier: Optional first node identifier to use. + :param coordinate_field_name: Optional name of coordinate field to define, if omitted use "coordinates". + :param field_names_and_values: Optional lists of (field_name, list of values) for additional fields to + define on the datapoints. Values may be scalar or vector (list of lists) real, or string. + Must be the same number of values as the number of points. + :return: next available node identifier + """ + fieldmodule = nodeset.getFieldmodule() + with ChangeManager(fieldmodule): + coordinates = find_or_create_field_coordinates(fieldmodule, name=coordinate_field_name) + + master_nodeset = nodeset.getMasterNodeset() + node_identifier = start_node_identifier if start_node_identifier is not None else \ + max(get_maximum_node_identifier(master_nodeset), 0) + 1 + + nodetemplate = master_nodeset.createNodetemplate() + nodetemplate.defineField(coordinates) + fields_values = [] # (field, is_string, values) + for field_name, field_values in field_names_and_values: + is_string = isinstance(field_values[0], str) + if is_string: + field = find_or_create_field_stored_string(fieldmodule, field_name, managed=True) + else: + components_count = len(field_values[0]) if isinstance(field_values[0], list) else 1 + field = find_or_create_field_finite_element(fieldmodule, field_name, components_count, managed=True) + nodetemplate.defineField(field) + fields_values.append((field, is_string, field_values)) + + fieldcache = fieldmodule.createFieldcache() + for n, x in enumerate(px): + node = nodeset.createNode(node_identifier, nodetemplate) + fieldcache.setNode(node) + coordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, x) + for field, is_string, values in fields_values: + if is_string: + field.assignString(fieldcache, values[n]) + else: + field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, values[n]) + node_identifier += 1 + + return node_identifier + + createCubeElement = create_cube_element createSquareElement = create_square_element createLineElement = create_line_element diff --git a/src/cmlibs/utils/zinc/region.py b/src/cmlibs/utils/zinc/region.py index 98c9e8f..16ca4da 100644 --- a/src/cmlibs/utils/zinc/region.py +++ b/src/cmlibs/utils/zinc/region.py @@ -1,4 +1,6 @@ from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, get_identifiers, get_maximum_node_identifier +from cmlibs.utils.zinc.regionbufferio import ( + region_read_from_buffer as read_from_buffer, region_write_to_buffer as write_to_buffer) from cmlibs.zinc.field import Field from cmlibs.zinc.result import RESULT_OK @@ -136,37 +138,3 @@ def determine_appropriate_glyph_size(region, coordinates): del fieldcache return glyph_width - - -def write_to_buffer(region, resource_domain_type=None, field_names=None): - """ - Write the contents of the given region to a buffer. - The content written to the buffer can be controlled with resource domain type and field names. - - :param region: The Zinc Region to write the content from. - :param resource_domain_type: Should be either Field.DOMAIN_TYPE_DATAPOINTS or Field.DOMAIN_TYPE_NODES. - :param field_names: A list of field names to output. - :return: The buffer in bytes written from the region. - """ - sir = region.createStreaminformationRegion() - srm = sir.createStreamresourceMemory() - if resource_domain_type is not None: - sir.setResourceDomainTypes(srm, resource_domain_type) - if field_names is not None: - sir.setFieldNames(field_names) - region.write(sir) - result, buffer = srm.getBuffer() - return buffer if result == RESULT_OK else None - - -def read_from_buffer(region, buffer): - """ - Read the given buffer into the given region. - - :param region: A Zinc Region to read the buffer into. - :param buffer: A buffer of bytes. - :return: The result of the region read operation. - """ - temp_sir = region.createStreaminformationRegion() - temp_sir.createStreamresourceMemoryBuffer(buffer) - return region.read(temp_sir) diff --git a/src/cmlibs/utils/zinc/regionbufferio.py b/src/cmlibs/utils/zinc/regionbufferio.py new file mode 100644 index 0000000..5f9828c --- /dev/null +++ b/src/cmlibs/utils/zinc/regionbufferio.py @@ -0,0 +1,45 @@ +""" +Utilities for reading and writing regions/field to/from a buffer. +""" +from cmlibs.zinc.result import RESULT_OK + + +def region_read_from_buffer(region, buffer): + """ + Read the given buffer into the given region. + + :param region: A Zinc Region to read the buffer into. + :param buffer: A buffer of bytes. + :return: The result of the region read operation. + """ + temp_sir = region.createStreaminformationRegion() + temp_sir.createStreamresourceMemoryBuffer(buffer) + return region.read(temp_sir) + + +def region_write_to_buffer(region, resource_domain_type=None, field_names=None, + recursion_mode=None): + """ + Write the contents of the given region to a buffer. + The content written to the buffer can be controlled with resource domain type and field names. + + :param region: The Zinc Region to write the content from. + :param resource_domain_type: A logical OR of any values in Field.DOMAIN_TYPE_* e.g. + Field.DOMAIN_TYPE_DATAPOINTS or Field.DOMAIN_TYPE_NODES or Field.DOMAIN_TYPE_MESH1D | Field.DOMAIN_TYPE_NODES. + :param field_names: Optional list of field names to output, if omitted writes all fields. + :param recursion_mode: Optionally specify mode controlling whether subregions are also output, currently either + StreaminformationRegion.RECURSION_MODE_ON (default) to write all subregions, or + StreaminformationRegion.RECURSION_MODE_OFF to write only this region and omit all groups. + :return: The buffer in bytes written from the region, or None if failed. + """ + sir = region.createStreaminformationRegion() + srm = sir.createStreamresourceMemory() + if resource_domain_type is not None: + sir.setResourceDomainTypes(srm, resource_domain_type) + if field_names is not None: + sir.setFieldNames(field_names) + if recursion_mode is not None: + sir.setRecursionMode(recursion_mode) + region.write(sir) + result, buffer = srm.getBuffer() + return buffer if result == RESULT_OK else None diff --git a/tests/test_zinc_finiteelement.py b/tests/test_zinc_finiteelement.py index 1641582..5b49a22 100644 --- a/tests/test_zinc_finiteelement.py +++ b/tests/test_zinc_finiteelement.py @@ -1,8 +1,14 @@ import math import unittest -from cmlibs.utils.zinc.finiteelement import define_grid_field_on_mesh, create_element_from_node_identifiers +from cmlibs.utils.zinc.field import find_or_create_field_coordinates +from cmlibs.utils.zinc.finiteelement import ( + create_element_from_node_identifiers, create_hermite_curve_elements, create_node_points, define_grid_field_on_mesh) +from cmlibs.utils.zinc.general import ChangeManager + from cmlibs.zinc.context import Context +from cmlibs.zinc.field import Field +from cmlibs.zinc.node import Node from cmlibs.zinc.result import RESULT_OK try: @@ -264,6 +270,92 @@ def test_create_element_from_node_identifiers(self): element_identifier = create_element_from_node_identifiers(mesh, coordinates_field, [134, 516]) self.assertEqual(22, element_identifier) + def test_create_hermite_curve(self): + """ + Test finite element utility to create a hermite curve mesh. + """ + ax = [ + [20614.63786414724, -2939.873650361369, -605.4150310438267], + [19866.346579952475, -4452.963739972778, -662.9986048963924], + [18688.55297975595, -5664.535791350198, -711.1634865902475]] + ad1 = [ + [-490.70985052735557, -1622.8873929937577, -60.65382802572525], + [-984.464558558513, -1380.6496536715067, -53.615308968309854], + [-1342.5638510703225, -1036.3223530813373, -42.38892426110936]] + context = Context("create hermite curve") + region = context.getDefaultRegion() + fieldmodule = region.getFieldmodule() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + with ChangeManager(fieldmodule): + group = fieldmodule.createFieldGroup() + group_name = "vagus_branch" + group.setName(group_name) + mesh1d = fieldmodule.findMeshByDimension(1) + mesh_group = group.getOrCreateMeshGroup(mesh1d) + next_node_identifier, next_element_identifier = create_hermite_curve_elements(mesh_group, ax, ad1) + self.assertEqual(next_node_identifier, 4) + self.assertEqual(next_element_identifier, 3) + nodeset_group = group.getNodesetGroup(nodes) + self.assertTrue(nodeset_group.isValid()) + + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + self.assertEqual(nodeset_group.getSize(), 3) + self.assertEqual(mesh_group.getSize(), 2) + fieldcache = fieldmodule.createFieldcache() + for n in range(3): + node = nodes.findNodeByIdentifier(n + 1) + self.assertTrue(node.isValid()) + fieldcache.setNode(node) + result, x = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, 3) + self.assertEqual(result, RESULT_OK) + result, d1 = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) + self.assertEqual(result, RESULT_OK) + for c in range(3): + self.assertAlmostEqual(x[c], ax[n][c], delta=1.0E-8) + self.assertAlmostEqual(d1[c], ad1[n][c], delta=1.0E-8) + + def test_create_node_points(self): + """ + Test creating node points with additional fields using utility. + """ + px = [[0.0, 1.1, 2.2], [1.0, 2.1, 3.2], [2.0, 3.1, 4.2]] + rv = [1.5, 3.0, 4.5] + sv = ["bob", "fred", "joe"] + context = Context("create node points") + region = context.getDefaultRegion() + fieldmodule = region.getFieldmodule() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + with ChangeManager(fieldmodule): + group = fieldmodule.createFieldGroup() + group_name = "group" + group.setName(group_name) + nodeset_group = group.createNodesetGroup(nodes) + next_node_identifier = create_node_points( + nodeset_group, px, field_names_and_values=[("reals", rv), ("strings", sv)]) + + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + reals = fieldmodule.findFieldByName("reals").castFiniteElement() + self.assertTrue(reals.isValid()) + strings = fieldmodule.findFieldByName("strings").castStoredString() + self.assertTrue(strings.isValid()) + self.assertEqual(nodeset_group.getSize(), 3) + fieldcache = fieldmodule.createFieldcache() + for n in range(3): + node = nodes.findNodeByIdentifier(n + 1) + self.assertTrue(node.isValid()) + fieldcache.setNode(node) + result, x = coordinates.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, 3) + self.assertEqual(result, RESULT_OK) + result, r = reals.getNodeParameters(fieldcache, -1, Node.VALUE_LABEL_VALUE, 1, 1) + self.assertEqual(result, RESULT_OK) + s = strings.evaluateString(fieldcache) + for c in range(3): + self.assertAlmostEqual(x[c], px[n][c], delta=1.0E-12) + self.assertAlmostEqual(r, rv[n], delta=1.0E-12) + self.assertEqual(s, sv[n]) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_zinc_region.py b/tests/test_zinc_region.py index 558fc9b..57d943a 100644 --- a/tests/test_zinc_region.py +++ b/tests/test_zinc_region.py @@ -4,8 +4,11 @@ from cmlibs.utils.zinc.field import find_or_create_field_coordinates from cmlibs.utils.zinc.finiteelement import create_nodes from cmlibs.utils.zinc.region import convert_nodes_to_datapoints, copy_nodeset +from cmlibs.utils.zinc.regionbufferio import region_read_from_buffer, region_write_to_buffer from cmlibs.zinc.context import Context from cmlibs.zinc.field import Field +from cmlibs.zinc.result import RESULT_OK +from cmlibs.zinc.streamregion import StreaminformationRegion class ZincRegionTestCase(unittest.TestCase): @@ -361,6 +364,35 @@ def test_transfer_nodes_VI(self): # This means we only lose the 'alt_coordinates' field from the source region. self.assertEqual(5, field_count) + def test_read_write_buffer_recursion_mode(self): + context = Context("test") + region = context.createRegion() + buffer = b"EX Version: 3\nRegion: /\nGroup name: bob\nRegion: /fred\nRegion: /joe\n" + result = region_read_from_buffer(region, buffer) + self.assertEqual(result, RESULT_OK) + fieldmodule = region.getFieldmodule() + bob = fieldmodule.findFieldByName("bob").castGroup() + self.assertTrue(bob.isValid()) + fred = region.findChildByName("fred") + self.assertTrue(fred.isValid()) + joe = region.findChildByName("joe") + self.assertTrue(joe.isValid()) + buffer2 = region_write_to_buffer(region, recursion_mode=StreaminformationRegion.RECURSION_MODE_ON) + self.assertTrue(buffer2 is not None) + self.assertEqual(buffer2, buffer) + buffer3 = region_write_to_buffer(region, recursion_mode=StreaminformationRegion.RECURSION_MODE_OFF) + self.assertTrue(buffer3 is not None) + if b"EX Version 3" in buffer3: + self.assertEqual(buffer3, b"EX Version: 3\nRegion: /\n") + result = region_read_from_buffer(fred, buffer3) + self.assertEqual(result, RESULT_OK) + fieldmodule2 = fred.getFieldmodule() + bob2 = fieldmodule2.findFieldByName("bob").castGroup() + self.assertFalse(bob2.isValid()) + fred2 = fred.findChildByName("fred") + self.assertFalse(fred2.isValid()) + joe2 = fred.findChildByName("joe") + self.assertFalse(joe2.isValid()) if __name__ == "__main__":