Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add choice of arithmetic or harmonic mean derivatives #26

Merged
merged 2 commits into from
Oct 17, 2018
Merged

Add choice of arithmetic or harmonic mean derivatives #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
08d0984
Add choice of arithmetic or harmonic mean scaling
rchristie Oct 17, 2018
bf93f4c
Use harmonic mean derivatives radially on inlet annulus
rchristie Oct 17, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions scaffoldmaker/utils/eftfactory_tricubichermite.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1352,13 +1352,14 @@ def createAnnulusMesh3d(self,
# smooth derivative 1 around inner loop
pd1[0][n2] = smoothCubicHermiteDerivativesLoop(px[0][n2], pd1[0][n2])

for n3 in range(1, nodesCountWall):
for n3 in range(0, nodesCountWall):
# smooth derivative 2 radially/along annulus
for n1 in range(nodesCountAround):
sd2 = smoothCubicHermiteDerivativesLine(
[ px [n3][n2][n1] for n2 in range(elementsCountRadial + 1) ],
[ pd2[n3][n2][n1] for n2 in range(elementsCountRadial + 1) ],
fixStartDerivative = True, fixEndDerivative = True)
fixStartDerivative = True, fixEndDerivative = True,
magnitudeScalingMode = DerivativeScalingMode.HARMONIC_MEAN)
for n2 in range(elementsCountRadial + 1):
pd2[n3][n2][n1] = sd2[n2]

Expand Down
57 changes: 42 additions & 15 deletions scaffoldmaker/utils/interpolation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
from __future__ import division
import collections
import copy
from enum import Enum
import math
import scaffoldmaker.utils.vector as vector

Expand Down Expand Up @@ -388,26 +389,37 @@ def getCubicHermiteCurvesPointAtArcDistance(nx, nd, arcDistance):
length += arcLength
return nx[-1], nd[-1], elementsCount - 1, 1.0

class DerivativeScalingMode(Enum):
ARITHMETIC_MEAN = 1 # derivative is half of sum of arclengths on either side
HARMONIC_MEAN = 2 # derivative is reciprocal of arithmetic mean of reciprocals of arclengths = arc lengths weighted by proportion from other side

def smoothCubicHermiteDerivativesLine(nx, nd1,
fixAllDirections = False,
fixStartDerivative = False, fixEndDerivative = False,
fixStartDirection = False, fixEndDirection = False):
fixStartDirection = False, fixEndDirection = False,
magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN):
"""
Modifies derivatives nd1 to be smoothly varying and near arc length.
Values are treated as being in a line.
Assumes initial 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 fixAllDirections: Set to True to only smooth magnitudes, otherwise both direction and magnitude are adjusted.
:param fixStartDerivative, fixEndDerivative: Set to True to fix derivative direction and magnitude at respective end.
:param fixStartDirection, fixEndDirection: Set to True to fix direction at respective end.
Redundant if fixAllDirections or respective fixStart/EndDerivative is True.
:param magnitudeScalingMode: A value from enum DerivativeScalingMode specifying
expression used to get derivative magnitude from adjacent arc lengths.
:return: Modified nd1
"""
nodesCount = len(nx)
elementsCount = nodesCount - 1
assert elementsCount > 0, 'smoothCubicHermiteDerivativesLine. Too few nodes/elements'
assert len(nd1) == nodesCount, 'smoothCubicHermiteDerivativesLine. Mismatched number of derivatives'
arithmeticMeanMagnitude = magnitudeScalingMode is DerivativeScalingMode.ARITHMETIC_MEAN
assert arithmeticMeanMagnitude or (magnitudeScalingMode is DerivativeScalingMode.HARMONIC_MEAN), \
'smoothCubicHermiteDerivativesLine. Invalid magnitude scaling mode'
md1 = copy.copy(nd1)
componentsCount = len(nx[0])
componentRange = range(componentsCount)
Expand All @@ -431,17 +443,21 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
# middle
for n in range(1, nodesCount - 1):
nm = n - 1
arcLengthmp = arcLengths[nm] + arcLengths[n]
wm = arcLengths[n ]/arcLengthmp
wp = arcLengths[nm]/arcLengthmp
if not fixAllDirections:
# average direction, weighted by fraction towards that end
np = (n + 1)%nodesCount
# get mean of directions from point n to points (n - 1) and (n + 1)
np = n + 1
dirm = [ (nx[n ][c] - nx[nm][c]) for c in componentRange ]
dirp = [ (nx[np][c] - nx[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
md1[n] = [ (wm*dirm[c] + wp*dirp[c]) for c in componentRange ]
# harmonic mean of magnitude
md1[n] = vector.setMagnitude(md1[n], wm*arcLengths[nm] + wp*arcLengths[n])
if arithmeticMeanMagnitude:
mag = 0.5*(arcLengths[nm] + arcLengths[n])
else: # harmonicMeanMagnitude
mag = 2.0/(1.0/arcLengths[nm] + 1.0/arcLengths[n])
md1[n] = vector.setMagnitude(md1[n], mag)
# end
if not fixEndDerivative:
if fixAllDirections or fixEndDirection:
Expand All @@ -453,19 +469,26 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
return md1

def smoothCubicHermiteDerivativesLoop(nx, nd1,
fixAllDirections = False):
fixAllDirections = False,
magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN):
"""
Modifies derivatives nd1 to be smoothly varying and near arc length.
Values are treated as being in a loop, so the first point follows the last.
Assumes initial 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 fixAllDirections: Set to True to only smooth magnitudes, otherwise both direction and magnitude are adjusted.
:param magnitudeScalingMode: A value from enum DerivativeScalingMode specifying
expression used to get derivative magnitude from adjacent arc lengths.
:return: Modified nd1
"""
nodesCount = elementsCount = len(nx)
assert elementsCount > 1, 'smoothCubicHermiteDerivativesLoop. Too few nodes/elements'
assert len(nd1) == elementsCount, 'smoothCubicHermiteDerivativesLoop. Mismatched number of derivatives'
arithmeticMeanMagnitude = magnitudeScalingMode is DerivativeScalingMode.ARITHMETIC_MEAN
assert arithmeticMeanMagnitude or (magnitudeScalingMode is DerivativeScalingMode.HARMONIC_MEAN), \
'smoothCubicHermiteDerivativesLine. Invalid magnitude scaling mode'
md1 = copy.copy(nd1)
componentsCount = len(nx[0])
componentRange = range(componentsCount)
Expand All @@ -482,17 +505,21 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1,
return md1
for n in range(nodesCount):
nm = n - 1
arcLengthmp = arcLengths[nm] + arcLengths[n]
wm = arcLengths[n ]/arcLengthmp
wp = arcLengths[nm]/arcLengthmp
if not fixAllDirections:
# average direction, weighted by fraction towards that end
# get mean of directions from point n to points (n - 1) and (n + 1)
np = (n + 1)%nodesCount
dirm = [ (nx[n ][c] - nx[nm][c]) for c in componentRange ]
dirp = [ (nx[np][c] - nx[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
md1[n] = [ (wm*dirm[c] + wp*dirp[c]) for c in componentRange ]
# harmonic mean of magnitude
md1[n] = vector.setMagnitude(md1[n], wm*arcLengths[nm] + wp*arcLengths[n])
if arithmeticMeanMagnitude:
mag = 0.5*(arcLengths[nm] + arcLengths[n])
else: # harmonicMeanMagnitude
mag = 2.0/(1.0/arcLengths[nm] + 1.0/arcLengths[n])
md1[n] = vector.setMagnitude(md1[n], mag)
lastArcLengths = arcLengths
print('smoothCubicHermiteDerivativesLoop max iters reached:',iter)
return md1
Expand Down