Merge pull request #50 from hsorby/main

Add calculate_jacobian and report_on_lowest_value functions.

Author: hsorby

src/cmlibs/utils/zinc/field.py

@@ -723,6 +723,51 @@ def find_coordinate_fields(region):
     return field_list
 
 
+def _is_3_component_real_valued_field(field):
+    return field is not None and field.isValid() and (field.getNumberOfComponents() == 3) and (field.getValueType() == Field.VALUE_TYPE_REAL)
+
+
+def create_jacobian_determinant_field(coordinates, reference_coordinates, name=None):
+    """
+    Create the Jacobian determinant of a 3-component coordinate field
+    w.r.t. 3-D element reference coordinates.
+    This value should always be positive for valid right-handed elements, negative if the
+    element volume becomes negative or is left-handed w.r.t. the reference coordinates.
+    If the coordinate field or reference field is not suitable for calculating
+    the Jacobian determinant, the function returns None.
+    :param coordinates: Geometric coordinate field.
+    :param reference_coordinates: Reference geometric coordinate field.
+    :param name: String to set as the name of the field.
+    :return: Jacobian determinant field.
+    """
+    jacobian_determinant = None
+    if _is_3_component_real_valued_field(coordinates) and _is_3_component_real_valued_field(reference_coordinates):
+        fm = coordinates.getFieldmodule()
+        with ChangeManager(fm):
+            jacobian_determinant = fm.createFieldDeterminant(fm.createFieldGradient(coordinates, reference_coordinates))
+            if name is not None:
+                jacobian_determinant.setName(name)
+
+    return jacobian_determinant
+
+
+def create_xi_reference_jacobian_determinant_field(coordinates, name=None):
+    """
+    Calculate the Jacobian determinant w.r.t the 'xi' field.
+    Returns None if the coordinates field is not a valid field.
+    See also :func:create_jacobian_determinant_field.
+    :param coordinates: Geometric coordinate field.
+    :param name: String to set as the name of the field, optional defaults to 'Jacobian_determinant_wrt_xi'.
+    :return: Jacobian determinant field.
+    """
+    jacobian_determinant = None
+    if isinstance(coordinates, Field) and coordinates.isValid():
+        fm = coordinates.getFieldmodule()
+        jacobian_determinant = create_jacobian_determinant_field(coordinates, fm.findFieldByName("xi"), "Jacobian_determinant_wrt_xi" if name is None else name)
+
+    return jacobian_determinant
+
+
 # Create C++ style aliases for names of functions.
 createFieldsDisplacementGradients = create_fields_displacement_gradients
 createFieldsTransformations = create_fields_transformations

src/cmlibs/utils/zinc/finiteelement.py

@@ -1,6 +1,8 @@
 """
 Utilities for creating and working with Zinc Finite Elements.
 """
+import math
+
 from cmlibs.maths import vectorops
 from cmlibs.utils.zinc.general import ChangeManager
 from cmlibs.zinc.element import Element, Elementbasis, Elementfieldtemplate, Mesh
@@ -617,6 +619,47 @@ def interpolate_cubic_hermite_derivative(v1, d1, v2, d2, xi):
     return [(f1 * v1[i] + f2 * d1[i] + f3 * v2[i] + f4 * d2[i]) for i in range(len(v1))]
 
 
+def _is_real_scalar_field(field):
+    return field is not None and field.isValid() and (field.getValueType() == Field.VALUE_TYPE_REAL) and (field.getNumberOfComponents() == 1)
+
+
+def get_scalar_field_minimum_in_mesh(real_field, mesh=None):
+    """
+    Evaluate the 3D mesh associated with the given (real) field and
+    report back on the element with the lowest field value. If the optional
+    mesh (or mesh_group) is not given then the highest dimension mesh available will be used.
+    Returns (-1, inf) if the evaluation is invalid.
+    :param real_field: A real scalar value field to find the lowest value of.
+    :param mesh: A mesh or mesh_group (optional, default None).
+    :return: A tuple of element identifier and minimum value within that element.
+    """
+    minimum_element_id = -1
+    minimum_element_value = math.inf
+    if _is_real_scalar_field(real_field):
+        fm = real_field.getFieldmodule()
+        fc = fm.createFieldcache()
+        fr = fc.createFieldrange()
+
+        in_use_mesh = mesh if mesh else get_highest_dimension_mesh(fm)
+        if in_use_mesh is None or not hasattr(in_use_mesh, "createElementiterator"):
+            return minimum_element_id, minimum_element_value
+
+        element_iter = in_use_mesh.createElementiterator()
+        element = element_iter.next()
+        while element.isValid():
+            fc.setElement(element)
+            result = real_field.evaluateFieldrange(fc, fr)
+            if result == RESULT_OK:
+                result, min_value, max_value = fr.getRangeReal(1)
+                if result == RESULT_OK and min_value < minimum_element_value:
+                    minimum_element_id = element.getIdentifier()
+                    minimum_element_value = min_value
+
+            element = element_iter.next()
+
+    return minimum_element_id, minimum_element_value
+
+
 createCubeElement = create_cube_element
 createSquareElement = create_square_element
 createLineElement = create_line_element

src/cmlibs/utils/zinc/mesh.py

@@ -181,7 +181,6 @@ def _find_connected(element_nodes, seed_index=None, ignore_element_indices=None,
 
     update_indexes = {}
     if progress_callback is not None:
-
         update_interval = max(1, int(num_elements * 0.01))
         update_indexes = set([i for i in range(update_interval)] + [i for i in range(update_interval, num_elements, update_interval)])