Merge pull request #15 from rchristie/ventricles2

Add atria 2 scaffold

scaffoldmaker/annotation/annotationgroup.py

@@ -64,7 +64,7 @@ def addSubelements(self):
             elementGroup = self._group.getFieldElementGroup(mesh)
             if elementGroup.isValid():
                 meshGroup = elementGroup.getMeshGroup()
-                print('Mesh group:', self._name, ', size', meshGroup.getSize())
+                #print('Mesh group:', self._name, ', size', meshGroup.getSize())
                 meshGroup.addElementsConditional(elementGroup)  # use FieldElementGroup as conditional field
 
 

scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py

@@ -684,20 +684,18 @@ def generateMesh(region, options):
         for elementId in [ ivcElementId, svcElementId ]:
             element = mesh.findElementByIdentifier(elementId)
             vcLength = vcInnerDiameter*0.5
-            tricubichermite.replaceElementWithInlet4(element, elementIdentifier, nodetemplate, nodeIdentifier, vcLength, vcInnerDiameter, vcWallThickness)
+            tricubichermite.replaceElementWithInlet4(element, elementIdentifier, nodetemplate, nodeIdentifier, vcLength, vcInnerDiameter*0.5, vcWallThickness)
             elementIdentifier += 4
             nodeIdentifier += 8
-            mesh.destroyElement(element)
 
         # add left atria inlets (pulmonary veins)
 
         for elementId in [ lapvElementId, lppvElementId, rapvElementId, rppvElementId ]:
             element = mesh.findElementByIdentifier(elementId)
             pvLength = pvInnerDiameter*0.5
-            tricubichermite.replaceElementWithInlet4(element, elementIdentifier, nodetemplate, nodeIdentifier, pvLength, pvInnerDiameter, pvWallThickness)
+            tricubichermite.replaceElementWithInlet4(element, elementIdentifier, nodetemplate, nodeIdentifier, pvLength, pvInnerDiameter*0.5, pvWallThickness)
             elementIdentifier += 4
             nodeIdentifier += 8
-            mesh.destroyElement(element)
 
         fm.endChange()
 

scaffoldmaker/scaffoldmaker.py

@@ -10,6 +10,7 @@
 from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1
 from scaffoldmaker.meshtypes.meshtype_3d_boxhole1 import MeshType_3d_boxhole1
 from scaffoldmaker.meshtypes.meshtype_3d_heartatria1 import MeshType_3d_heartatria1
+from scaffoldmaker.meshtypes.meshtype_3d_heartatria2 import MeshType_3d_heartatria2
 from scaffoldmaker.meshtypes.meshtype_3d_heartventricles1 import MeshType_3d_heartventricles1
 from scaffoldmaker.meshtypes.meshtype_3d_heartventricles2 import MeshType_3d_heartventricles2
 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase1 import MeshType_3d_heartventriclesbase1
@@ -30,6 +31,7 @@ def __init__(self):
             MeshType_3d_box1,
             MeshType_3d_boxhole1,
             MeshType_3d_heartatria1,
+            MeshType_3d_heartatria2,
             MeshType_3d_heartventricles1,
             MeshType_3d_heartventricles2,
             MeshType_3d_heartventriclesbase1,

scaffoldmaker/utils/geometry.py

@@ -24,7 +24,7 @@ def getEllipseArcLength(a, b, angle1Radians, angle2Radians):
     :param b: Minor axis length.(On y, PI/2, 3PI/2).
     :param angle1Radians: First angle anticlockwise from major axis.
     :param angle2Radians: Second angle anticlockwise from major axis.
-    :return: Perimeter lenge, positive if anticlockwise, otherwise negative.
+    :return: Perimeter length, positive if anticlockwise, otherwise negative.
     '''
     angle1 = min(angle1Radians, angle2Radians)
     angle2 = max(angle1Radians, angle2Radians)
@@ -57,7 +57,7 @@ def updateEllipseAngleByArcLength(a, b, inAngleRadians, arcLength):
     angle = inAngleRadians
     lengthMoved = 0.0
     lengthTol = (a + b)*1.0E-4  # broader tolerance due to reliance on inexact getEllipseArcLength()
-    #print('inAngleRadians', inAngleRadians, ', arcLength', arcLength)
+    #print('updateEllipseAngleByArcLength', a, b, 'inAngleRadians', inAngleRadians, ', arcLength', arcLength)
     while math.fabs(arcLength - lengthMoved) > lengthTol:
         t = ( -a*math.sin(angle), b*math.cos(angle) )
         dlength_dangle = math.sqrt(t[0]*t[0] + t[1]*t[1])
@@ -66,3 +66,32 @@ def updateEllipseAngleByArcLength(a, b, inAngleRadians, arcLength):
         #print('lengthMoved', lengthMoved)
     #print('updateEllipseAngleByArcLength a', a, 'b', b, ', angle', inAngleRadians, ', arcLength', arcLength, ' -> ', angle)
     return angle
+
+def getEllipseRadiansToX(ax, bx, dx, initialTheta):
+    '''
+    Iteratively compute theta in ax*cos(theta) + bx*sin(theta) = dx
+    from initialTheta. Uses Newton's method.
+    :param ax: x component of major axis.
+    :param bx: x component of minor axis.
+    :param dx: x distance from centre of ellipse.
+    :param initialTheta: Initial value of theta needed since multi-valued.
+    :return: theta in radians
+    '''
+    theta = initialTheta
+    iters = 0
+    fTol = math.sqrt(ax*ax + bx*bx)*1.0E-10
+    while True:
+        cosAngle = math.cos(theta)
+        sinAngle = math.sin(theta)
+        f = ax*cosAngle + bx*sinAngle - dx
+        if math.fabs(f) < fTol:
+            break;
+        df = -ax*sinAngle + bx*cosAngle
+        #print(iters, '. theta', theta, 'f', f,'df',df,'-->',theta - f/df)
+        theta -= f/df
+        iters += 1
+        if iters == 100:
+            print('getEllipseRadiansToX: did not converge!')
+            break
+    return theta
+

scaffoldmaker/utils/zinc_utils.py

@@ -124,3 +124,31 @@ def interpolateNodesCubicHermite(cache, coordinates, xi, normal_scale, \
     dx_ds_normal = [ normal_scale*d for d in radialVector ]
 
     return x, dx_ds, dx_ds_cross, dx_ds_normal
+
+def computeNodeDerivativeHermiteLagrange(cache, coordinates, node1, derivative1, scale1, node2, scale2):
+    """
+    Computes the derivative at node2 from quadratic Hermite-Lagrange interpolation of
+    node1 value and derivative1 to node2 value.
+    :param cache: Field cache to evaluate in.
+    :param coordinates: Coordinates field.
+    :param node1, node2: Start and end nodes.
+    :param derivative1: Node value label for derivative at node1.
+    :param scale1, scale2: Scaling to apply to derivatives at nodes, e.g. -1.0 to reverse.
+    :return: dx_dxi at node2
+    """
+    cache.setNode(node1)
+    result, v1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3 )
+    result, d1 = coordinates.getNodeParameters(cache, -1, derivative1, 1, 3 )
+    d1 = [ d*scale1 for d in d1 ]
+    cache.setNode(node2)
+    result, v2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3 )
+    xi = 1.0
+    #phi1 = 1 - xi*xi
+    #phi2 = xi - xi*xi
+    #phi3 = xi*xi
+    dphi1 = -2.0*xi
+    dphi2 = 1 - 2.0*xi
+    dphi3 = 2.0*xi
+    d2 = [ (v1[c]*dphi1 + d1[c]*dphi2 + v2[c]*dphi3) for c in range(3) ]
+    d2 = [ d*scale2 for d in d2 ]
+    return d2