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

Linear basis through wall for sphere shell #23

Merged
merged 6 commits into from
Oct 5, 2018
Merged

Linear basis through wall for sphere shell #23

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b55d5c7
Add linear basis through wall
mlin865 Oct 2, 2018
f38f55f
Merge remote-tracking branch 'abi/master' into bicubicSphereShell
mlin865 Oct 2, 2018
0df8652
Rename useCubicHermite as useCubicHermiteThroughWall
mlin865 Oct 4, 2018
9945ffa
Apply duck typing on eftfactory
mlin865 Oct 4, 2018
71bd05a
Correct function name for creating shell top and bottom in other mesh…
mlin865 Oct 4, 2018
02162c6
Delete lines made redundant from duck typing
mlin865 Oct 5, 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
6 changes: 3 additions & 3 deletions scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
54 changes: 36 additions & 18 deletions scaffoldmaker/meshtypes/meshtype_3d_sphereshell1.py
View file
Open in desktop
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,10 @@ def generateBaseMesh(region, options):
for e1 in range(elementsCountAround):
va = e1
vb = (e1 + 1)%elementsCountAround
eft1 = tricubichermite.createEftShellApexBottom(va*100, vb*100)
if useCubicHermiteThroughWall:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These 4 lines should just be:

eft1 = eftfactory.createEftShellPoleBottom(va*100, vb*100)

as the same code is called in both cases.

eft1 = eftfactory.createEftShellPoleBottom(va*100, vb*100)
else:
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 +372,10 @@ def generateBaseMesh(region, options):
for e1 in range(elementsCountAround):
va = e1
vb = (e1 + 1)%elementsCountAround
eft1 = tricubichermite.createEftShellApexTop(va*100, vb*100)
if useCubicHermiteThroughWall:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These 4 lines should just be:

eft1 = eftfactory.createEftShellPoleTop(va*100, vb*100)

as the same code is called in both cases.

eft1 = eftfactory.createEftShellPoleTop(va*100, vb*100)
else:
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
Open in desktop
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