Add path smooth cross derivatives function

src/scaffoldmaker/meshtypes/meshtype_1d_path1.py

@@ -4,13 +4,13 @@
 
 from __future__ import division
 import math
-from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates
+from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup
 from opencmiss.utils.zinc.general import ChangeManager
 from opencmiss.zinc.element import Element, Elementbasis
 from opencmiss.zinc.field import Field
 from opencmiss.zinc.node import Node
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
-from scaffoldmaker.utils.interpolation import DerivativeScalingMode, smoothCubicHermiteDerivativesLine
+from scaffoldmaker.utils.interpolation import DerivativeScalingMode, smoothCubicHermiteDerivativesLine, smoothCubicHermiteCrossDerivativesLine
 from scaffoldmaker.utils import vector
 from opencmiss.zinc.result import RESULT_OK
 
@@ -131,49 +131,70 @@ def generateBaseMesh(cls, region, options):
         return []
 
     @classmethod
-    def smoothPath(cls, region, options, mode : DerivativeScalingMode):
+    def smoothPath(cls, region, options, editGroupName, mode : DerivativeScalingMode):
         x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1])
         d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode)
-        setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ])
+        setPathParameters(region, [ Node.VALUE_LABEL_D_DS1 ], [ d1 ], editGroupName)
+        return False, True  # settings not changed, nodes changed
 
     @classmethod
-    def makeD2Normal(cls, region, options):
+    def makeD2Normal(cls, region, options, editGroupName):
         if not options['D2 derivatives']:
             return
         d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2])
         for c in range(len(d1)):
             td2 = vector.vectorRejection(d2[c], d1[c])
             d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c]))
-        setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2])
+        setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2], editGroupName)
+        return False, True  # settings not changed, nodes changed
 
     @classmethod
-    def makeD3Normal(cls, region, options):
+    def makeD3Normal(cls, region, options, editGroupName):
         if not options['D3 derivatives']:
             return
         if options['D2 derivatives']:
             d1, d2, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2,
                                                                   Node.VALUE_LABEL_D_DS3])
             for c in range(len(d1)):
                 d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c]))
-            setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3])
         else:
             d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3])
             for c in range(len(d1)):
                 td3 = vector.vectorRejection(d3[c], d1[c])
                 d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c]))
-            setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3])
+        setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3], editGroupName)
+        return False, True  # settings not changed, nodes changed
 
+    @classmethod
+    def smoothCrossDX(cls, region, options, editGroupName, valueLabel):
+        if valueLabel == Node.VALUE_LABEL_D_DS2:
+            if not options['D2 derivatives']:
+                return
+            crossDerivativeLabel = Node.VALUE_LABEL_D2_DS1DS2
+        elif valueLabel == Node.VALUE_LABEL_D_DS3:
+            if not options['D3 derivatives']:
+                return
+            crossDerivativeLabel = Node.VALUE_LABEL_D2_DS1DS3
+        else:
+            assert False, 'Invalid value label'
+        valueLabels = [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, valueLabel, crossDerivativeLabel]
+        x, d1, d2, d12 = extractPathParametersFromRegion(region, valueLabels)
+        d12 = smoothCubicHermiteCrossDerivativesLine(x, d1, d2, d12)
+        setPathParameters(region, [ crossDerivativeLabel ], [ d12 ], editGroupName)
+        return False, True  # settings not changed, nodes changed
 
     @classmethod
     def getInteractiveFunctions(cls):
         """
         Supply client with functions for smoothing path parameters.
         """
         return [
-            ("Smooth D1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)),
-            ("Smooth D1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)),
-            ("Make D2 normal", lambda region, options: cls.makeD2Normal(region, options)),
-            ("Make D3 normal", lambda region, options: cls.makeD3Normal(region, options))
+            ("Smooth D1 arithmetic", lambda region, options, editGroupName: cls.smoothPath(region, options, editGroupName, DerivativeScalingMode.ARITHMETIC_MEAN)),
+            ("Smooth D1 harmonic", lambda region, options, editGroupName: cls.smoothPath(region, options, editGroupName, DerivativeScalingMode.HARMONIC_MEAN)),
+            ("Make D2 normal", lambda region, options, editGroupName: cls.makeD2Normal(region, options, editGroupName)),
+            ("Make D3 normal", lambda region, options, editGroupName: cls.makeD3Normal(region, options, editGroupName)),
+            ("Smooth D2", lambda region, options, editGroupName: cls.smoothCrossDX(region, options, editGroupName, Node.VALUE_LABEL_D_DS2)),
+            ("Smooth D3", lambda region, options, editGroupName: cls.smoothCrossDX(region, options, editGroupName, Node.VALUE_LABEL_D_DS3))
         ]
 
 
@@ -217,11 +238,12 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None):
     return returnValues
 
 
-def setPathParameters(region, nodeValueLabels, nodeValues):
+def setPathParameters(region, nodeValueLabels, nodeValues, editGroupName=None):
     '''
     Set node parameters for coordinates field in path from listed values.
     :param nodeValueLabels: List of nodeValueLabels to set e.g. [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ]
     :param nodeValues: List of values for each type e.g. [ xlist, d1list ]
+    :param editGroupName: Optional name of existing or new Zinc group to record modified nodes in.
     '''
     fieldmodule = region.getFieldmodule()
     coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
@@ -232,6 +254,12 @@ def setPathParameters(region, nodeValueLabels, nodeValues):
     nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
     assert nodesCount == nodes.getSize()
     with ChangeManager(fieldmodule):
+        if editGroupName:
+            editGroup = findOrCreateFieldGroup(fieldmodule, editGroupName, managed=False)
+            editNodeGroup = editGroup.getFieldNodeGroup(nodes)
+            if not editNodeGroup.isValid():
+                editNodeGroup = editGroup.createFieldNodeGroup(nodes)
+            editNodesetGroup = editNodeGroup.getNodesetGroup()
         cache = fieldmodule.createFieldcache()
         nodeIter = nodes.createNodeiterator()
         node = nodeIter.next()
@@ -240,5 +268,7 @@ def setPathParameters(region, nodeValueLabels, nodeValues):
             cache.setNode(node)
             for v in range(nodeValueLabelsCount):
                 coordinates.setNodeParameters(cache, -1, nodeValueLabels[v], 1, nodeValues[v][n])
+            if editGroupName:
+                editNodesetGroup.addNode(node)
             node = nodeIter.next()
             n += 1

src/scaffoldmaker/meshtypes/scaffold_base.py

@@ -155,6 +155,13 @@ def getInteractiveFunctions(cls):
         """
         Override to return list of named interactive functions that client
         can invoke to modify mesh parameters with a push button control.
-        :return: list(tuples), (name : str, callable(region, options)).
+        Functions must take 3 arguments: Zinc region, scaffold options and
+        editGroupName (can be None) for optional name of Zinc group to
+        create or modify so changed nodes etc. are put in it.
+        Functions return 2 boolean values: optionsChanged, nodesChanged.
+        These tell the client whether to redisplay the options or process
+        the effects of node edits (which will be recorded in edit group if
+        its name is supplied).
+        :return: list(tuples), (name : str, callable(region, options, editGroupName)).
         """
         return []

src/scaffoldmaker/utils/interpolation.py

@@ -719,12 +719,82 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
     cmax = 0.0
     for n in range(nodesCount):
         for c in componentRange:
-            if math.fabs(md1[n][c] - lastmd1[n][c]) > cmax:
-                cmax = math.fabs(md1[n][c] - lastmd1[n][c])
+            cmax = max(cmax, math.fabs(md1[n][c] - lastmd1[n][c]))
     closeness = cmax / dtol
     print('smoothCubicHermiteDerivativesLine max iters reached:', iter + 1, ', cmax = ', round(closeness,2), 'x tolerance')
     return md1
 
+def smoothCubicHermiteCrossDerivativesLine(nx, nd1, nd2, nd12,
+        fixStartDerivative = False, fixEndDerivative = False, instrument=False):
+    """
+    Smooth derivatives of cross directions of hermite curves.
+    Assumes initial nd12 derivatives are zero or reasonable.
+    Where directions are smoothed the weighted/harmonic mean is used.
+    :param nx: List of coordinates of nodes along curves.
+    :param nd1: List of derivatives of nodes along curves.
+    :param nd2: List of lateral direction vectors of nodes along curves.
+    :param nd12: List of derivatives of lateral directions along curves.
+    :param fixStartDerivative, fixEndDerivative: Set to True to fix derivative direction and magnitude at respective end.
+    :return: Modified nd12
+    """
+    nodesCount = len(nx)
+    elementsCount = nodesCount - 1
+    assert elementsCount > 0, 'smoothCubicHermiteCrossDerivativesLine.  Too few nodes/elements'
+    assert len(nd1) == nodesCount, 'smoothCubicHermiteCrossDerivativesLine.  Mismatched number of derivatives'
+    md12 = copy.copy(nd12)
+    componentsCount = len(nx[0])
+    componentRange = range(componentsCount)
+    # special case where equal derivatives at each end are sought
+    if (elementsCount == 1) and not (fixStartDerivative or fixEndDerivative):
+        delta = [ (nd2[1][c] - nd2[0][c]) for c in componentRange ]
+        return [ delta, copy.deepcopy(delta) ]
+    tol = 1.0E-6
+    arcLengths = [ getCubicHermiteArcLength(nx[e], nd1[e], nx[e + 1], nd1[e + 1]) for e in range(elementsCount) ]
+    dtol = tol*sum(vector.magnitude(d) for d in nd2)
+    if instrument:
+        print('iter 0', md12)
+    for iter in range(100):
+        lastmd12 = copy.copy(md12)
+        # start
+        if not fixStartDerivative:
+            md12[0] = interpolateLagrangeHermiteDerivative(nd2[0], nd2[1], lastmd12[1], 0.0)
+        # middle
+        for n in range(1, nodesCount - 1):
+            nm = n - 1
+            # get mean of directions from point n to points (n - 1) and (n + 1)
+            np = n + 1
+            dirm = [ (nd2[n ][c] - nd2[nm][c]) for c in componentRange ]
+            dirp = [ (nd2[np][c] - nd2[n ][c]) for c in componentRange ]
+            # mean weighted by fraction towards that end, equivalent to harmonic mean
+            arcLengthmp = arcLengths[nm] + arcLengths[n]
+            wm = arcLengths[n ]/arcLengthmp
+            wp = arcLengths[nm]/arcLengthmp
+            md12[n] = [ (wm*dirm[c] + wp*dirp[c]) for c in componentRange ]
+        # end
+        if not fixEndDerivative:
+            md12[-1] = interpolateHermiteLagrangeDerivative(nd2[-2], lastmd12[-2], nd2[-1], 1.0)
+        if instrument:
+            print('iter', iter + 1, md12)
+        for n in range(nodesCount):
+            for c in componentRange:
+                if math.fabs(md12[n][c] - lastmd12[n][c]) > dtol:
+                    break
+            else:
+                continue
+            break
+        else:
+            if instrument:
+                print('smoothCubicHermiteCrossDerivativesLine converged after iter:', iter + 1)
+            return md12
+
+    cmax = 0.0
+    for n in range(nodesCount):
+        for c in componentRange:
+            cmax = max(cmax, math.fabs(md12[n][c] - lastmd12[n][c]))
+    closeness = cmax / dtol
+    print('smoothCubicHermiteCrossDerivativesLine max iters reached:', iter + 1, ', cmax = ', round(closeness,2), 'x tolerance')
+    return md12
+
 def smoothCubicHermiteDerivativesLoop(nx, nd1,
         fixAllDirections = False,
         magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN, instrument=False):
@@ -790,8 +860,7 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1,
     cmax = 0.0
     for n in range(nodesCount):
         for c in componentRange:
-            if math.fabs(md1[n][c] - lastmd1[n][c]) > cmax:
-                cmax = math.fabs(md1[n][c] - lastmd1[n][c])
+            cmax = max(cmax, math.fabs(md1[n][c] - lastmd1[n][c]))
     closeness = cmax / dtol
     print('smoothCubicHermiteDerivativesLoop max iters reached:', iter + 1, ', cmax = ', round(closeness,2) , 'x tolerance')
     return md1