Add function to create hermite curve

src/cmlibs/utils/zinc/field.py

@@ -2,15 +2,14 @@
 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
 from cmlibs.zinc.fieldmodule import Fieldmodule
 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()

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

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)

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

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