Skip to content

Commit

Permalink
Merge pull request #257 from rchristie/determine_eft
Browse files Browse the repository at this point in the history 
Add utility to determine Hermite element field template
  • Loading branch information
@rchristie
rchristie authored May 15, 2024
2 parents 5c960ab + 9b9faa3 commit 887d70b
Show file tree
Hide file tree
Showing 3 changed files with 256 additions and 4 deletions.
139 changes: 137 additions & 2 deletions src/scaffoldmaker/utils/eft_utils.py
View file
Original file line number Diff line number Diff line change
@@ -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):
Expand Down Expand Up @@ -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) / 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
4 changes: 2 additions & 2 deletions src/scaffoldmaker/utils/interpolation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
117 changes: 117 additions & 0 deletions tests/test_general.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, \
Expand Down Expand Up @@ -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()

0 comments on commit 887d70b

Please sign in to comment.