From a453891d5b6cac59c598366d64c21a62ddbf7843 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 9 Oct 2020 18:16:17 +1300 Subject: [PATCH] Add path smooth cross derivatives function Add edit group name option to interactive functions and enhance return values. --- .../meshtypes/meshtype_1d_path1.py | 58 ++++++++++---- src/scaffoldmaker/meshtypes/scaffold_base.py | 9 ++- src/scaffoldmaker/utils/interpolation.py | 77 ++++++++++++++++++- 3 files changed, 125 insertions(+), 19 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 35618f55..4575c744 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/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,23 +131,25 @@ 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']: @@ -155,14 +157,31 @@ def makeD3Normal(cls, region, options): 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): @@ -170,10 +189,12 @@ 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 diff --git a/src/scaffoldmaker/meshtypes/scaffold_base.py b/src/scaffoldmaker/meshtypes/scaffold_base.py index 978493d8..af4ebb0f 100644 --- a/src/scaffoldmaker/meshtypes/scaffold_base.py +++ b/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 [] diff --git a/src/scaffoldmaker/utils/interpolation.py b/src/scaffoldmaker/utils/interpolation.py index 36cb866d..67f5408d 100644 --- a/src/scaffoldmaker/utils/interpolation.py +++ b/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