Skip to content

Commit

Permalink
Merge pull request #23 from mlin865/bicubicSphereShell
Browse files Browse the repository at this point in the history 
Linear basis through wall for sphere shell
  • Loading branch information
@rchristie
rchristie authored Oct 5, 2018
2 parents ca31652 + 02162c6 commit 77a0f33
Show file tree
Hide file tree
Showing 6 changed files with 154 additions and 46 deletions.
6 changes: 3 additions & 3 deletions scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -863,7 +863,7 @@ def generateBaseMesh(cls, region, options):
radiansAround0 = 0.5*math.pi + s*radiansPerElementAroundEnd
radiansAround1 = radiansAround0 + radiansPerElementAroundEnd
# scale factor identifiers follow convention of offsetting by 100 for each 'version'
eft1 = tricubichermite.createEftShellApexTop(s*100, (s + 1)*100)
eft1 = tricubichermite.createEftShellPoleTop(s*100, (s + 1)*100)
scalefactors = [
-1.0,
math.cos(radiansAround0), math.sin(radiansAround0), radiansPerElementAroundEnd,
Expand Down Expand Up @@ -960,7 +960,7 @@ def generateBaseMesh(cls, region, options):
radiansAround0 = -0.5*math.pi + s*radiansPerElementAroundEnd
radiansAround1 = radiansAround0 + radiansPerElementAroundEnd
# scale factor identifiers follow convention of offsetting by 100 for each 'version'
eft1 = tricubichermite.createEftShellApexTop(s*100, (s + 1)*100)
eft1 = tricubichermite.createEftShellPoleTop(s*100, (s + 1)*100)
scalefactors = [
-1.0,
math.cos(radiansAround0), math.sin(radiansAround0), radiansPerElementAroundEnd,
Expand Down Expand Up @@ -1131,7 +1131,7 @@ def generateBaseMesh(cls, region, options):
for e1 in range(elementsCountAroundFossa):
va = e1
vb = (e1 + 1)%elementsCountAroundFossa
eft1 = tricubichermite.createEftShellApexTop(va*100, vb*100)
eft1 = tricubichermite.createEftShellPoleTop(va*100, vb*100)
elementtemplateX.defineField(coordinates, -1, eft1)
element = mesh.createElement(elementIdentifier, elementtemplateX)
nids = [ fossaNodeId[0][va], fossaNodeId[0][vb], fossaCentreNodeId[0], fossaNodeId[1][va], fossaNodeId[1][vb], fossaCentreNodeId[1] ]
Expand Down
4 changes: 2 additions & 2 deletions scaffoldmaker/meshtypes/meshtype_3d_heartatria2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1054,7 +1054,7 @@ def generateBaseMesh(cls, region, options):
radiansAroundNext = radiansAroundApex + radiansPerElementAroundApex
va = e1
vb = e1 + 1
eft1 = tricubichermite.createEftShellApexTop(s*100, (s + 1)*100)
eft1 = tricubichermite.createEftShellPoleTop(s*100, (s + 1)*100)
elementtemplateX.defineField(coordinates, -1, eft1)
element = mesh.createElement(elementIdentifier, elementtemplateX)
nodeIdentifiers = [ aNodeId[0][n2][va], aNodeId[0][n2][vb], apexNodeId[0], aNodeId[1][n2][va], aNodeId[1][n2][vb], apexNodeId[1] ]
Expand Down Expand Up @@ -1193,7 +1193,7 @@ def generateBaseMesh(cls, region, options):
for e1 in range(elementsCountAroundFossa):
va = e1
vb = (e1 + 1)%elementsCountAroundFossa
eft1 = tricubichermite.createEftShellApexTop(va*100, vb*100)
eft1 = tricubichermite.createEftShellPoleTop(va*100, vb*100)
elementtemplateX.defineField(coordinates, -1, eft1)
element = mesh.createElement(elementIdentifier, elementtemplateX)
nids = [ fossaNodeId[0][va], fossaNodeId[0][vb], fossaCentreNodeId[0], fossaNodeId[1][va], fossaNodeId[1][vb], fossaCentreNodeId[1] ]
Expand Down
2 changes: 1 addition & 1 deletion scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -653,7 +653,7 @@ def generateBaseMesh(cls, region, options):
vb = (e1 + 1)%elementsCountAroundLV
nids = [ 1 , 2 + va , 2 + vb,
1 + nowl, 2 + va + nowl, 2 + vb + nowl ]
eft1 = tricubichermite.createEftShellApexBottom(va*100, vb*100)
eft1 = tricubichermite.createEftShellPoleBottom(va*100, vb*100)
# calculate general linear map coefficients
if e1 < elementsCountAroundVSeptum:
deltaRadians = dRadians1 = dRadians2 = radiansPerElementAroundSeptum
Expand Down
48 changes: 30 additions & 18 deletions scaffoldmaker/meshtypes/meshtype_3d_sphereshell1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

from __future__ import division
import math
from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
from scaffoldmaker.utils.zinc_utils import *
from scaffoldmaker.utils.meshrefinement import MeshRefinement
Expand Down Expand Up @@ -33,6 +34,7 @@ def getDefaultOptions():
'Length ratio' : 1.0,
'Element length ratio equator/apex' : 1.0,
'Use cross derivatives' : False,
'Use linear through wall' : False,
'Refine' : False,
'Refine number of elements around' : 1,
'Refine number of elements up' : 1,
Expand All @@ -52,6 +54,7 @@ def getOrderedOptionNames():
'Length ratio',
'Element length ratio equator/apex',
'Use cross derivatives',
'Use linear through wall',
'Refine',
'Refine number of elements around',
'Refine number of elements up',
Expand Down Expand Up @@ -89,7 +92,8 @@ def checkOptions(options):
@staticmethod
def generateBaseMesh(region, options):
"""
Generate the base tricubic Hermite mesh. See also generateMesh().
Generate the base tricubic Hermite or bicubic Hermite linear mesh.
See also generateMesh().
:param region: Zinc region to define model in. Must be empty.
:param options: Dict containing options. See getDefaultOptions().
:return: None
Expand All @@ -98,6 +102,7 @@ def generateBaseMesh(region, options):
elementsCountUp = options['Number of elements up']
elementsCountThroughWall = options['Number of elements through wall']
useCrossDerivatives = options['Use cross derivatives']
useCubicHermiteThroughWall = not(options['Use linear through wall'])
excludeBottomRows = options['Exclude bottom rows']
excludeTopRows = options['Exclude top rows']
wallThickness = options['Wall thickness']
Expand All @@ -115,27 +120,30 @@ def generateBaseMesh(region, options):
nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)
if useCubicHermiteThroughWall:
nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)
if useCrossDerivatives:
nodetemplate = nodes.createNodetemplate()
nodetemplate.defineField(coordinates)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1)
if useCubicHermiteThroughWall:
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1)
else:
nodetemplate = nodetemplateApex

mesh = fm.findMeshByDimension(3)

tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives)
eft = tricubichermite.createEftBasic()

tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE)
if useCubicHermiteThroughWall:
eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives)
else:
eftfactory = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives)
eft = eftfactory.createEftBasic()

elementtemplate = mesh.createElementtemplate()
elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE)
Expand Down Expand Up @@ -232,7 +240,8 @@ def generateBaseMesh(region, options):
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, position[1] ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, vector2[0], 0.0 ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ vector2[0], 0.0, 0.0 ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, vector3[1] ])
if useCubicHermiteThroughWall:
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, vector3[1] ])
nodeIdentifier = nodeIdentifier + 1

elif n2 < elementsCountUp:
Expand Down Expand Up @@ -266,12 +275,14 @@ def generateBaseMesh(region, options):
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3)
if useCubicHermiteThroughWall:
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3)
if useCrossDerivatives:
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero)
if useCubicHermiteThroughWall:
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero)
nodeIdentifier = nodeIdentifier + 1

else:
Expand All @@ -281,7 +292,8 @@ def generateBaseMesh(region, options):
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, position[1] ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, vector2[0], 0.0 ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ -vector2[0], 0.0, 0.0 ])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, vector3[1] ])
if useCubicHermiteThroughWall:
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, vector3[1] ])
nodeIdentifier = nodeIdentifier + 1

# create elements
Expand Down Expand Up @@ -311,7 +323,7 @@ def generateBaseMesh(region, options):
for e1 in range(elementsCountAround):
va = e1
vb = (e1 + 1)%elementsCountAround
eft1 = tricubichermite.createEftShellApexBottom(va*100, vb*100)
eft1 = eftfactory.createEftShellPoleBottom(va*100, vb*100)
elementtemplate1.defineField(coordinates, -1, eft1)
element = mesh.createElement(elementIdentifier, elementtemplate1)
bni1 = no + 1
Expand Down Expand Up @@ -357,7 +369,7 @@ def generateBaseMesh(region, options):
for e1 in range(elementsCountAround):
va = e1
vb = (e1 + 1)%elementsCountAround
eft1 = tricubichermite.createEftShellApexTop(va*100, vb*100)
eft1 = eftfactory.createEftShellPoleTop(va*100, vb*100)
elementtemplate1.defineField(coordinates, -1, eft1)
element = mesh.createElement(elementIdentifier, elementtemplate1)
bni3 = no + now
Expand Down
96 changes: 96 additions & 0 deletions scaffoldmaker/utils/eftfactory_bicubichermitelinear.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,102 @@ def createEftNoCrossDerivatives(self):
assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftNoCrossDerivatives: Failed to validate eft'
return eft

def createEftShellPoleBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1):
'''
Create a bicubic hermite linear element field for closing bottom pole of a shell.
Element is collapsed in xi1 on xi2 = 0.
Each collapsed node has 3 scale factors giving the cos, sin coefficients
of the radial line from global derivatives, plus the arc subtended by
the element in radians, so the pole can be rounded.
Need to create a new template for each sector around pole giving common
nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and
add 100 for each radial line around pole.
:param nodeScaleFactorOffset0: offset of node scale factors at pole on xi1=0
:param nodeScaleFactorOffset1: offset of node scale factors at pole on xi1=1
:return: Element field template
'''
# start with full bicubic hermite linear to remap D2_DS1DS2 at pole
eft = self._mesh.createElementfieldtemplate(self._basis)
if not self._useCrossDerivatives:
for n in [ 2, 3, 6, 7 ]:
eft.setFunctionNumberOfTerms(n*4 + 4, 0)

# GRC: allow scale factor identifier for global -1.0 to be prescribed
setEftScaleFactorIds(eft, [1], [
nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3,
nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ])
# remap parameters before collapsing nodes
remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [])
for layer in range(2):
so = layer*6 + 1
ln = layer*4 + 1
# 2 terms for d/dxi2 via general linear map:
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS2, [so + 2]) ])
# 2 terms for cross derivative 1 2 to correct circular pole: -sin(theta).phi, cos(theta).phi
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3]), (Node.VALUE_LABEL_D_DS2, [1, so + 1, so + 3]) ])

ln = layer*4 + 2
# 2 terms for d/dxi2 via general linear map:
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 4]), (Node.VALUE_LABEL_D_DS2, [so + 5]) ])
# 2 terms for cross derivative 1 2 to correct circular pole: -sin(theta).phi, cos(theta).phi
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 5, so + 6]), (Node.VALUE_LABEL_D_DS2, [1, so + 4, so + 6]) ])

ln_map = [ 1, 1, 2, 3, 4, 4, 5, 6 ]
remapEftLocalNodes(eft, 6, ln_map)

assert eft.validate(), 'eftfactory_tricubichermite.createEftShellPoleBottom: Failed to validate eft'
return eft

def createEftShellPoleTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1):
'''
Create a bicubic hermite linear element field for closing top pole of a shell.
Element is collapsed in xi1 on xi2 = 1.
Each collapsed node has 3 scale factors giving the cos, sin coefficients
of the radial line from global derivatives, plus the arc subtended by
the element in radians, so the pole can be rounded.
Need to create a new template for each sector around pole giving common
nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and
add 100 for each radial line around pole.
:param nodeScaleFactorOffset0: offset of node scale factors at pole on xi1=0
:param nodeScaleFactorOffset1: offset of node scale factors at pole on xi1=1
:return: Element field template
'''
# start with full bicubic hermite linear to remap D2_DS1DS2 at pole
eft = self._mesh.createElementfieldtemplate(self._basis)
if not self._useCrossDerivatives:
for n in [ 0, 1, 4, 5 ]:
eft.setFunctionNumberOfTerms(n*4 + 4, 0)

# GRC: allow scale factor identifier for global -1.0 to be prescribed
setEftScaleFactorIds(eft, [1], [
nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3,
nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3,
nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ])
# remap parameters before collapsing nodes
remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS1, [])
for layer in range(2):
so = layer*6 + 1
ln = layer*4 + 3
# 2 terms for d/dxi2 via general linear map:
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS2, [so + 2]) ])
# 2 terms for cross derivative 1 2 to correct circular pole: -sin(theta).phi, cos(theta).phi
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3]) ])

ln = layer*4 + 4
# 2 terms for d/dxi2 via general linear map:
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS2, so + 5) ])
# 2 terms for cross derivative 1 2 to correct circular pole: -sin(theta).phi, cos(theta).phi
remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS2, [so + 4, so + 6]) ])

ln_map = [ 1, 2, 3, 3, 4, 5, 6, 6 ]
remapEftLocalNodes(eft, 6, ln_map)

assert eft.validate(), 'eftfactory_tricubichermite.createEftShellPoleTop: Failed to validate eft'
return eft

def createEftSplitXi1RightStraight(self):
'''
Create an element field template suitable for the inner elements of the
Expand Down
Loading

0 comments on commit 77a0f33

Please sign in to comment.