Skip to content

Commit

Permalink
Add LV outlet elements
Browse files Browse the repository at this point in the history
  • Loading branch information
rchristie committed Mar 13, 2018
1 parent 4fcf5ff commit 6049e3b
Show file tree
Hide file tree
Showing 2 changed files with 151 additions and 5 deletions.
45 changes: 40 additions & 5 deletions scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
from scaffoldmaker.utils.interpolation import *
from scaffoldmaker.utils.zinc_utils import *
from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
from opencmiss.zinc.element import Element, Elementbasis
from opencmiss.zinc.field import Field
from opencmiss.zinc.node import Node
Expand Down Expand Up @@ -90,7 +91,7 @@ def generateMesh(region, options):
"""
elementsCountUp = options['Number of elements up']
elementsCountAround = options['Number of elements around']
print('elementsCountAround', elementsCountAround)
#print('elementsCountAround', elementsCountAround)
elementsCountAcrossSeptum = options['Number of elements across septum']
lvWallThickness = options['LV wall thickness']
lvWallThicknessRatioBase = options['LV wall thickness ratio base']
Expand Down Expand Up @@ -344,7 +345,7 @@ def generateMesh(region, options):
qx = ox*laUnitMajorX + oy*laUnitMajorY
qy = ox*laUnitMinorX + oy*laUnitMinorY
cruxLeftRadians = math.atan2(atriaOuterMajorMag*qy, atriaOuterMinorMag*qx)
print('qx', qx, ' qy', qy, ' cruxLeftDegrees', cruxLeftRadians*180.0/math.pi)
#print('qx', qx, ' qy', qy, ' cruxLeftDegrees', cruxLeftRadians*180.0/math.pi)

cosRadiansAround = math.cos(laSeptumRadians)
sinRadiansAround = math.sin(laSeptumRadians)
Expand Down Expand Up @@ -422,7 +423,7 @@ def generateMesh(region, options):
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3)
nodeIdentifier += 1

if True:
if False:
# show axes of left atrium
node = nodes.createNode(nodeIdentifier, nodetemplateFull)
cache.setNode(node)
Expand Down Expand Up @@ -466,7 +467,7 @@ def generateMesh(region, options):
cruxRightRadians = math.atan2(atriaOuterMajorMag*qy, atriaOuterMinorMag*qx)
if cruxRightRadians < math.pi:
cruxRightRadians += 2.0*math.pi
print('qx', qx, ' qy', qy, ' cruxRightDegrees', cruxRightRadians*180.0/math.pi)
#print('qx', qx, ' qy', qy, ' cruxRightDegrees', cruxRightRadians*180.0/math.pi)

raRadians = [0.0]*elementsCountAroundAtria
raRadians[0] = raSeptumRadians
Expand Down Expand Up @@ -538,7 +539,7 @@ def generateMesh(region, options):
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3)
nodeIdentifier += 1

if True:
if False:
# show axes of right atrium
node = nodes.createNode(nodeIdentifier, nodetemplateFull)
cache.setNode(node)
Expand Down Expand Up @@ -955,5 +956,39 @@ def generateMesh(region, options):
print('create element rv', elementIdentifier, result, result2, result3, nids[e])
elementIdentifier += 1

# LV outlet ring elements: linear in xi3

bicubichermitelinear3 = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives, linearAxis = 3)

# wedge shape, collapsed across xi3 at xi2 == 0
eft_lvring = bicubichermitelinear3.createEftNoCrossDerivatives()
ln_map = [ 1, 2, 3, 4, 1, 2, 5, 6 ]
remapEftLocalNodes(eft_lvring, 6, ln_map)
lvring_nids = [ nidl + 2, nidl + 3, nidl + 4, nidl + 5, nidl + 6, laNodeId[0][1] ]
for e in range(elementsCountAroundOutlet):
e2 = (e + 1) % elementsCountAroundOutlet
eft1 = eft_lvring
if e >= 4:
eft1 = bicubichermitelinear3.createEftNoCrossDerivatives()
setEftScaleFactorIds(eft1, [1], [])
if e == 4:
remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS1, [ (Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS2, []) ])
localNodes = [ 2, 6 ] if (e == 4) else [ 1, 5 ]
remapEftNodeValueLabel(eft1, localNodes, Node.VALUE_LABEL_D_DS1, [ (Node.VALUE_LABEL_D_DS1, [1]) ])
remapEftNodeValueLabel(eft1, localNodes[0:1], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ])
remapEftNodeValueLabel(eft1, localNodes[1:], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, []) ])
ln_map = [ 1, 2, 3, 4, 1, 2, 5, 6 ]
remapEftLocalNodes(eft1, 6, ln_map)

elementtemplate1 = mesh.createElementtemplate()
elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE)
result = elementtemplate1.defineField(coordinates, -1, eft1)

element = mesh.createElement(elementIdentifier, elementtemplate1)
nids = [ lvring_nids[e], lvring_nids[e2], lvOutletNodeId[0][e], lvOutletNodeId[0][e2], lvOutletNodeId[1][e], lvOutletNodeId[1][e2] ]
result2 = element.setNodesByIdentifier(eft1, nids)
print('create element lv outlet', elementIdentifier, result, result2, nids)
elementIdentifier += 1

fm.endChange()

111 changes: 111 additions & 0 deletions scaffoldmaker/utils/eftfactory_bicubichermitelinear.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
'''
Definitions of standard element field templates using bicubic Hermite x linear Lagrange basis.
@author: Richard Christie
'''
from scaffoldmaker.utils.eft_utils import *
from opencmiss.zinc.element import Elementbasis, Elementfieldtemplate
from opencmiss.zinc.node import Node
from opencmiss.zinc.status import OK as ZINC_OK

class eftfactory_bicubichermitelinear:
'''
Factory class for creating element field templates for a 3-D mesh using bicubic Hermite x linear Lagrange basis.
'''

def __init__(self, mesh, useCrossDerivatives, linearAxis = 3,
d_ds1 = Node.VALUE_LABEL_D_DS1, d_ds2 = Node.VALUE_LABEL_D_DS2):
'''
:param mesh: Zinc mesh to create element field templates in.
:param useCrossDerivatives: Set to True if you want cross derivative terms.
:param linearAxis: 1, 2, or 3.
:param d_ds1: Node derivative to use in Hermite axis 1: Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2.
:param d_ds2: Node derivative to use in Hermite axis 2, > d_ds1: Node.VALUE_LABEL_D_DS2 or Node.VALUE_LABEL_D_DS3.
'''
assert mesh.getDimension() == 3, 'eftfactory_bicubichermitelinear: not a 3-D Zinc mesh'
assert linearAxis in [ 1, 2, 3 ], 'eftfactory_bicubichermitelinear: linearAxis must be 1, 2 or 3'
assert d_ds1 in [ Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2 ], 'eftfactory_bicubichermitelinear: invalid d_ds1'
assert d_ds2 in [ Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ] and (d_ds2 > d_ds1), 'eftfactory_bicubichermitelinear: invalid d_ds2'
self._mesh = mesh
self._useCrossDerivatives = useCrossDerivatives
self._linearAxis = linearAxis
self._d_ds1 = d_ds1
self._d_ds2 = d_ds2
self._d2_ds1ds2 = Node.VALUE_LABEL_D2_DS2DS3 if (d_ds1 == Node.VALUE_LABEL_D_DS2) \
else Node.VALUE_LABEL_D2_DS1DS3 if (d_ds2 == Node.VALUE_LABEL_D_DS3) \
else Node.VALUE_LABEL_D2_DS1DS2
self._fieldmodule = mesh.getFieldmodule()
self._basis = self._fieldmodule.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE)
self._basis.setFunctionType(linearAxis, Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE)

def _remapDefaultNodeDerivatives(self):
'''
Remap the Hermite node derivatives to those chosen in __init__.
Use only on first create.
'''
# must do d_ds2 first!
if self._d_ds2 != Node.VALUE_LABEL_D_DS2:
remapEftNodeValueLabel(eft, range(1, 9), Node.VALUE_LABEL_D_DS2, [ (self._d_ds2, []) ])
if self._d_ds1 != Node.VALUE_LABEL_D_DS1:
remapEftNodeValueLabel(eft, range(1, 9), Node.VALUE_LABEL_D_DS1, [ (self._d_ds1, []) ])
if self._d2_ds1ds2 != Node.VALUE_LABEL_D2_DS1DS2:
remapEftNodeValueLabel(eft, range(1, 9), Node.VALUE_LABEL_D2_DS1DS2, [ (self._d2_ds1ds2, []) ])

def createEftBasic(self):
'''
Create the basic biicubic Hermite x linear Lagrange element template with 1:1 mappings to
node derivatives ds1 & ds2, with or without cross derivatives accordinate as initialised.
:return: Element field template
'''
if not self._useCrossDerivatives:
return self.createEftNoCrossDerivatives()
eft = self._mesh.createElementfieldtemplate(self._basis)
self._remapDefaultNodeDerivatives()
assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftBasic: Failed to validate eft'
return eft

def createEftNoCrossDerivatives(self):
'''
Create a basic tricubic hermite element template with 1:1 mappings to
node derivatives ds1 & ds2, without cross derivatives.
:return: Element field template
'''
eft = self._mesh.createElementfieldtemplate(self._basis)
for n in range(8):
eft.setFunctionNumberOfTerms(n*4 + 4, 0)
self._remapDefaultNodeDerivatives()
assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftNoCrossDerivatives: Failed to validate eft'
return eft

def createEftSplitXi1RightStraight(self):
'''
Create an element field template suitable for the inner elements of the
join between left and right chambers, with xi1 bifurcating to right.
Straight through version.
Only works with linearAxis 2.
:return: Element field template
'''
assert linearAxis == 2, 'eftfactory_bicubichermitelinear.createEftSplitXi1RightStraight: Not linearAxis 2'
eft = self.createEftNoCrossDerivatives()
setEftScaleFactorIds(eft, [1], [])
remapEftNodeValueLabel(eft, [ 5, 7 ], self._d_ds1, [ (self._d_ds1, []), (self._d_ds2, [1]) ])
remapEftNodeValueLabel(eft, [ 6, 8 ], self._d_ds1, [ (self._d_ds1, []), (self._d_ds2, []) ])
assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftSplitXi1RightStraight: Failed to validate eft'
return eft

def createEftSplitXi1RightOut(self):
'''
Create an element field template suitable for the outer elements of the
join between left and right chambers, with xi1 bifurcating to right.
Right out version i.e. xi1 heading to right. h-shape.
Only works with linearAxis 2.
:return: Element field template
'''
assert linearAxis == 2, 'eftfactory_bicubichermitelinear.createEftSplitXi1RightOut: Not linearAxis 2'
eft = self.createEftNoCrossDerivatives()
setEftScaleFactorIds(eft, [1], [])
remapEftNodeValueLabel(eft, [ 1, 3 ], self._d_ds1, [ (self._d_ds1, [1]) ])
remapEftNodeValueLabel(eft, [ 1, 3 ], self._d_ds2, [ (self._d_ds1, [1]), (self._d_ds2, [1]) ])
remapEftNodeValueLabel(eft, [ 5, 7 ], self._d_ds2, [ (self._d_ds1, [1]), (self._d_ds2, []) ])
assert eft.validate(), 'eftfactory_bicubichermitelinear.createEftSplitXi1RightOut: Failed to validate eft'
return eft

0 comments on commit 6049e3b

Please sign in to comment.