From a020fb9a11eb2719cde38412c3fc396843066d8f Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 2 May 2024 15:30:48 +1200 Subject: [PATCH 1/3] Add utility to determine Hermite element field template --- src/scaffoldmaker/utils/eft_utils.py | 139 ++++++++++++++++++++++++++- tests/test_general.py | 117 ++++++++++++++++++++++ 2 files changed, 254 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index 05856d5c..0e7e73c5 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -1,9 +1,11 @@ ''' Utility functions for element field templates shared by mesh generators. ''' -from cmlibs.utils.zinc.finiteelement import getElementNodeIdentifiers -from cmlibs.zinc.element import Elementfieldtemplate +from cmlibs.maths.vectorops import add, dot, magnitude, mult, sub +from cmlibs.zinc.element import Elementbasis, Elementfieldtemplate +from cmlibs.zinc.node import Node from cmlibs.zinc.result import RESULT_OK +import math def getEftTermScaling(eft, functionIndex, termIndex): @@ -321,3 +323,136 @@ def createEftElementSurfaceLayer(elementIn, eftIn, eftfactory, eftStd, removeNod assert eft.validate(), "Element " + str(elementIn.getIdentifier()) + " eft not validated" return eft, scalefactors + + +def determineTricubicHermiteEft(mesh, nodeParameters, nodeDerivativeFixedWeights=None, serendipity=False, + mapCrossDerivatives=False): + """ + Determine the tricubic Hermite normal/serendipity element field template for + interpolating node parameters at corners of a cube, by matching deltas + between corners with node derivatives. + Node lists use zinc ordering which varies nodes fastest over lower element coordinate. + :param mesh: A Zinc mesh of dimension 3. + :param nodeParameters: List over 8 local nodes in Zinc ordering of 4 parameter + vectors x, d1, d2, d3 each with 3 components. + :param nodeDerivativeFixedWeights: Optional list over 8 local nodes in Zinc ordering of + list of up to 3 derivatives each containing weights for global d1, d2, d3 to fix that + derivative, or None to use default search. + Example: [[], [None, [0.0, -1.0, 1.0]], [], [], [], [], [], []] + Local nodes 1, 3-8 use the standard search for all derivatives. + Local node 2 forces d/dxi2 = -d2 + d3, d/dxi1 and d/dxi3 use the standard search. + :param serendipity: Set to True for serendipity basis. + :param mapCrossDerivatives: For non-serendipity Hermite basis, map cross derivatives + as appropriate for first derivatives. For example, if d/dxi1 = d2 and d/dxi2 = d3, + d2/dxi1dxi2 should map to -d23. If any element derivative is a sum of node derivatives, + eliminate cross derivatives which use that direction. + :return: eft, scale factors list [-1.0] or None. Returned eft can be further modified. + """ + if nodeDerivativeFixedWeights is None: + nodeDerivativeFixedWeights = [[]] * 8 + assert mesh.getDimension() == 3 + assert len(nodeParameters) == 8 + assert len(nodeParameters[0]) == 4 + assert len(nodeDerivativeFixedWeights) == 8 + assert not (serendipity and mapCrossDerivatives) + delta12 = sub(nodeParameters[1][0], nodeParameters[0][0]) + delta34 = sub(nodeParameters[3][0], nodeParameters[2][0]) + delta56 = sub(nodeParameters[5][0], nodeParameters[4][0]) + delta78 = sub(nodeParameters[7][0], nodeParameters[6][0]) + delta13 = sub(nodeParameters[2][0], nodeParameters[0][0]) + delta24 = sub(nodeParameters[3][0], nodeParameters[1][0]) + delta57 = sub(nodeParameters[6][0], nodeParameters[4][0]) + delta68 = sub(nodeParameters[7][0], nodeParameters[5][0]) + delta15 = sub(nodeParameters[4][0], nodeParameters[0][0]) + delta26 = sub(nodeParameters[5][0], nodeParameters[1][0]) + delta37 = sub(nodeParameters[6][0], nodeParameters[2][0]) + delta48 = sub(nodeParameters[7][0], nodeParameters[3][0]) + deltas = [ + [delta12, delta13, delta15], + [delta12, delta24, delta26], + [delta34, delta13, delta37], + [delta34, delta24, delta48], + [delta56, delta57, delta15], + [delta56, delta68, delta26], + [delta78, delta57, delta37], + [delta78, delta68, delta48] + ] + fieldmodule = mesh.getFieldmodule() + tricubicHermiteBasis = fieldmodule.createElementbasis( + 3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE_SERENDIPITY if serendipity + else Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + eft = mesh.createElementfieldtemplate(tricubicHermiteBasis) + scalefactors = None + derivativeLabels = [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3] + crossDerivativeLabels = [Node.VALUE_LABEL_D2_DS1DS2, Node.VALUE_LABEL_D2_DS1DS3, + Node.VALUE_LABEL_D2_DS2DS3, Node.VALUE_LABEL_D3_DS1DS2DS3] + for n in range(8): + ln = n + 1 + derivativeExpressionTerms = [] + for ed in range(3): + if mapCrossDerivatives: + functionNumber = n * 8 + derivativeLabels[ed] + else: + functionNumber = n * 4 + ed + 2 + expressionTerms = [] + if nodeDerivativeFixedWeights[n] and (ed < len(nodeDerivativeFixedWeights[n])): + weights = nodeDerivativeFixedWeights[n][ed] + if weights is not None: + for nd in range(len(weights)): + weight = weights[nd] + if weight: + expressionTerms.append((derivativeLabels[nd], [1] if (weight < 0.0) else [])) + if not expressionTerms: + delta = deltas[n][ed] # delta side coordinates for nominal derivative direction + magDelta = magnitude(delta) + greatestSimilarity = 0.0 # can be negative in which can + derivativeLabel = None + for nd in range(3): + derivative = nodeParameters[n][nd + 1] + magDerivative = magnitude(derivative) + cosineSimilarity = dot(derivative, delta) / (magDerivative * magDelta) + magnitudeSimilarity = math.exp(-math.fabs(magDerivative - magDelta)) + similarity = cosineSimilarity * magnitudeSimilarity + if math.fabs(similarity) > math.fabs(greatestSimilarity): + greatestSimilarity = similarity + derivativeLabel = derivativeLabels[nd] + expressionTerms.append((derivativeLabel, [1] if (greatestSimilarity < 0.0) else [])) + termCount = len(expressionTerms) + if termCount != 1: + eft.setFunctionNumberOfTerms(functionNumber, termCount) + for t in range(termCount): + term = t + 1 + eft.setTermNodeParameter(functionNumber, term, ln, expressionTerms[t][0], 1) + scaling = expressionTerms[t][1] + if scaling: + if not scalefactors: + # assumes only + setEftScaleFactorIds(eft, [1], []) + scalefactors = [-1.0] + eft.setTermScaling(functionNumber, term, scaling) + if mapCrossDerivatives: + derivativeExpressionTerms.append(expressionTerms) + if mapCrossDerivatives: + for cd in range(4): + functionNumber = n * 8 + crossDerivativeLabels[cd] + derivativeIndexes = \ + [0, 1] if (cd == 0) else \ + [0, 2] if (cd == 1) else \ + [1, 2] if (cd == 2) else \ + [0, 1, 2] + sign = 1.0 + crossDerivativeLabel = Node.VALUE_LABEL_VALUE + for ed in derivativeIndexes: + if len(derivativeExpressionTerms[ed]) != 1: + crossDerivativeLabel = Node.VALUE_LABEL_VALUE + break + crossDerivativeLabel += derivativeExpressionTerms[ed][0][0] - Node.VALUE_LABEL_VALUE + if derivativeExpressionTerms[ed][0][1]: + sign *= -1.0 + if crossDerivativeLabel == Node.VALUE_LABEL_VALUE: + eft.setFunctionNumberOfTerms(functionNumber, 0) + else: + eft.setTermNodeParameter(functionNumber, 1, ln, crossDerivativeLabel, 1) + if sign < 0.0: + eft.setTermScaling(functionNumber, 1, [1]) + return eft, scalefactors \ No newline at end of file diff --git a/tests/test_general.py b/tests/test_general.py index 7c4ea984..c0831b2a 100644 --- a/tests/test_general.py +++ b/tests/test_general.py @@ -20,6 +20,7 @@ from scaffoldmaker.scaffoldpackage import ScaffoldPackage from scaffoldmaker.scaffolds import Scaffolds from scaffoldmaker.utils.bifurcation import SegmentTubeData +from scaffoldmaker.utils.eft_utils import determineTricubicHermiteEft from scaffoldmaker.utils.geometry import getEllipsoidPlaneA, getEllipsoidPolarCoordinatesFromPosition, \ getEllipsoidPolarCoordinatesTangents from scaffoldmaker.utils.interpolation import computeCubicHermiteSideCrossDerivatives, evaluateCoordinatesOnCurve, \ @@ -1857,6 +1858,122 @@ def test_smooth_side_cross_derivatives(self): self.assertAlmostEqual(targetLength, actualLength, delta=LENGTH_TOL) # print("xi", xi, "length", actualLength, "angle", actualAngle, targetAngle) + def test_determineHermiteSerendipityEft(self): + """ + Test algorithm for determining hermite serendipity eft from node derivative directions. + """ + context = Context("test_determineHermiteSerendipityEft") + region = context.getDefaultRegion() + fieldmodule = region.getFieldmodule() + mesh3d = fieldmodule.findMeshByDimension(3) + + nodeParameters = [ + [[-1.064, 0.152, 0.035], [0.776, -0.051, -0.204], [0.000, 3.000, 0.000], [0.000, 0.000, 1.000]], + [[-0.227, -0.113, 0.004], [0.069, -0.003, 0.880], [0.000, 3.000, 0.000], [-0.920, 0.188, 0.028]], + [[0.000, 3.000, 0.000], [1.000, 0.000, 0.000], [0.000, 3.000, 0.000], [0.000, 0.000, 1.000]], + [[0.700, 2.967, -0.194], [0.075, -0.196, -0.726], [0.000, 3.000, 0.000], [0.793, 0.064, 0.082]], + [[-1.033, 0.135, 1.042], [-0.196, 2.329, 0.099], [-0.330, -0.052, 0.698], [0.785, -0.320, -0.149]], + [[-0.174, 0.001, 0.883], [1.000, 0.000, 0.000], [0.000, 3.000, 0.000], [0.000, 0.000, 1.000]], + [[-0.904, 2.968, 1.362], [1.657, -0.955, -0.204], [-1.299, 2.712, 0.142], [0.000, 0.000, 1.000]], + [[0.337, 3.541, 1.461], [1.000, 0.000, 0.000], [-0.681, 2.850, -0.167], [0.000, 0.000, 1.000]] + ] + nodeDerivativeFixedWeights = [ + [], + [], + [], + [], + [], + [], + [None, [-1.0, 1.0]], + [] + ] + eft, scalefactors = determineTricubicHermiteEft(mesh3d, nodeParameters, nodeDerivativeFixedWeights, + serendipity=True) + regularTermExpressions =\ + [[(Node.VALUE_LABEL_D_DS1, [])], [(Node.VALUE_LABEL_D_DS2, [])], [(Node.VALUE_LABEL_D_DS3, [])]] + expectedNodeTermExpressions = [ + regularTermExpressions, + [[(Node.VALUE_LABEL_D_DS3, [1])], [(Node.VALUE_LABEL_D_DS2, [])], [(Node.VALUE_LABEL_D_DS1, [])]], + regularTermExpressions, + [[(Node.VALUE_LABEL_D_DS3, [])], [(Node.VALUE_LABEL_D_DS2, [])], [(Node.VALUE_LABEL_D_DS1, [1])]], + [[(Node.VALUE_LABEL_D_DS3, [])], [(Node.VALUE_LABEL_D_DS1, [])], [(Node.VALUE_LABEL_D_DS2, [])]], + regularTermExpressions, + [[(Node.VALUE_LABEL_D_DS1, [])], [(Node.VALUE_LABEL_D_DS1, [1]), (Node.VALUE_LABEL_D_DS2, [])], + [(Node.VALUE_LABEL_D_DS3, [])]], + regularTermExpressions + ] + self.assertEqual(scalefactors, [-1.0]) + self.assertEqual(eft.getNumberOfLocalScaleFactors(), 1) + self.assertEqual(eft.getScaleFactorType(1), eft.SCALE_FACTOR_TYPE_GLOBAL_GENERAL) + self.assertEqual(eft.getScaleFactorIdentifier(1), 1) + for n in range(8): + ln = n + 1 + for ed in range(3): + functionNumber = n * 4 + ed + 2 + expectedTermExpression = expectedNodeTermExpressions[n][ed] + expectedTermCount = len(expectedTermExpression) + self.assertEqual(eft.getFunctionNumberOfTerms(functionNumber), expectedTermCount) + for t in range(expectedTermCount): + term = t + 1 + self.assertEqual(eft.getTermLocalNodeIndex(functionNumber, term), ln) + self.assertEqual(eft.getTermNodeValueLabel(functionNumber, term), expectedTermExpression[t][0]) + self.assertEqual(eft.getTermNodeVersion(functionNumber, term), 1) + expectedScaleCount = len(expectedTermExpression[t][1]) + actualScaleCount, scalefactorIndexes = eft.getTermScaling(functionNumber, term, expectedScaleCount) + self.assertEqual(actualScaleCount, expectedScaleCount) + if expectedScaleCount == 1: + self.assertEqual(scalefactorIndexes, 1) + else: + for s in range(expectedScaleCount): + self.assertEqual(scalefactorIndexes[s], 1) + + # test mapping cross derivatives for full tricubic Hermite + eft, scalefactors = determineTricubicHermiteEft(mesh3d, nodeParameters, nodeDerivativeFixedWeights, + serendipity=False, mapCrossDerivatives=True) + + regularCrossTermExpressions =\ + [[(Node.VALUE_LABEL_D2_DS1DS2, [])], [(Node.VALUE_LABEL_D2_DS1DS3, [])], + [(Node.VALUE_LABEL_D2_DS2DS3, [])], [(Node.VALUE_LABEL_D3_DS1DS2DS3, [])]] + expectedNodeCrossTermExpressions = [ + regularCrossTermExpressions, + [[(Node.VALUE_LABEL_D2_DS2DS3, [1])], [(Node.VALUE_LABEL_D2_DS1DS3, [1])], + [(Node.VALUE_LABEL_D2_DS1DS2, [])], [(Node.VALUE_LABEL_D3_DS1DS2DS3, [1])]], + regularCrossTermExpressions, + [[(Node.VALUE_LABEL_D2_DS2DS3, [])], [(Node.VALUE_LABEL_D2_DS1DS3, [1])], + [(Node.VALUE_LABEL_D2_DS1DS2, [1])], [(Node.VALUE_LABEL_D3_DS1DS2DS3, [1])]], + [[(Node.VALUE_LABEL_D2_DS1DS3, [])], [(Node.VALUE_LABEL_D2_DS2DS3, [])], + [(Node.VALUE_LABEL_D2_DS1DS2, [])], [(Node.VALUE_LABEL_D3_DS1DS2DS3, [])]], + regularCrossTermExpressions, + [[], [(Node.VALUE_LABEL_D2_DS1DS3, [])], [], []], + regularCrossTermExpressions + ] + crossDerivativeLabels = [Node.VALUE_LABEL_D2_DS1DS2, Node.VALUE_LABEL_D2_DS1DS3, + Node.VALUE_LABEL_D2_DS2DS3, Node.VALUE_LABEL_D3_DS1DS2DS3] + self.assertEqual(scalefactors, [-1.0]) + self.assertEqual(eft.getNumberOfLocalScaleFactors(), 1) + self.assertEqual(eft.getScaleFactorType(1), eft.SCALE_FACTOR_TYPE_GLOBAL_GENERAL) + self.assertEqual(eft.getScaleFactorIdentifier(1), 1) + for n in range(8): + ln = n + 1 + for cd in range(4): + functionNumber = n * 8 + crossDerivativeLabels[cd] + expectedTermExpression = expectedNodeCrossTermExpressions[n][cd] + expectedTermCount = len(expectedTermExpression) + self.assertEqual(eft.getFunctionNumberOfTerms(functionNumber), expectedTermCount) + for t in range(expectedTermCount): + term = t + 1 + self.assertEqual(eft.getTermLocalNodeIndex(functionNumber, term), ln) + self.assertEqual(eft.getTermNodeValueLabel(functionNumber, term), expectedTermExpression[t][0]) + self.assertEqual(eft.getTermNodeVersion(functionNumber, term), 1) + expectedScaleCount = len(expectedTermExpression[t][1]) + actualScaleCount, scalefactorIndexes = eft.getTermScaling(functionNumber, term, expectedScaleCount) + self.assertEqual(actualScaleCount, expectedScaleCount) + if expectedScaleCount == 1: + self.assertEqual(scalefactorIndexes, 1) + else: + for s in range(expectedScaleCount): + self.assertEqual(scalefactorIndexes[s], 1) + if __name__ == "__main__": unittest.main() From b15a05126ecc344684accc154c289c832c2b1745 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 15 May 2024 12:56:40 +1200 Subject: [PATCH 2/3] Make magnitude similarity non-dimensional --- src/scaffoldmaker/utils/eft_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index 0e7e73c5..93f5c45a 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -411,7 +411,7 @@ def determineTricubicHermiteEft(mesh, nodeParameters, nodeDerivativeFixedWeights derivative = nodeParameters[n][nd + 1] magDerivative = magnitude(derivative) cosineSimilarity = dot(derivative, delta) / (magDerivative * magDelta) - magnitudeSimilarity = math.exp(-math.fabs(magDerivative - magDelta)) + magnitudeSimilarity = math.exp(-math.fabs((magDerivative - magDelta) / magDelta)) similarity = cosineSimilarity * magnitudeSimilarity if math.fabs(similarity) > math.fabs(greatestSimilarity): greatestSimilarity = similarity From 9b9faa378e39cae79632c9af1a96a30521ddf544 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 15 May 2024 12:57:43 +1200 Subject: [PATCH 3/3] Improve cubic hermite basis function documentation --- src/scaffoldmaker/utils/interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/utils/interpolation.py b/src/scaffoldmaker/utils/interpolation.py index 6ba3ee6b..42a68b86 100644 --- a/src/scaffoldmaker/utils/interpolation.py +++ b/src/scaffoldmaker/utils/interpolation.py @@ -26,7 +26,7 @@ def getCubicHermiteBasis(xi): """ - :return: 4 basis functions for x1, d1, x2, d2 + :return: 4 cubic Hermite basis function values for x1, d1, x2, d2 at xi. """ xi2 = xi*xi xi3 = xi2*xi @@ -54,7 +54,7 @@ def interpolateCubicHermite(v1, d1, v2, d2, xi): def getCubicHermiteBasisDerivatives(xi): """ - :return: 4 derivatives of basis functions for x1, d1, x2, d2 + :return: 4 cubic Hermite basis function first derivative values for x1, d1, x2, d2 at xi. """ xi2 = xi*xi df1 = -6.0*xi + 6.0*xi2