Merge pull request #32 from rchristie/heart1

Add partial arterial root

scaffoldmaker/meshtypes/meshtype_3d_heart1.py

@@ -99,12 +99,13 @@ def generateBaseMesh(cls, region, options):
         annotationGroups += [ lFibrousRingGroup, rFibrousRingGroup ]
 
         # annotation points
-        #dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
+        dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
         dataLabel = getOrCreateLabelField(fm, 'data_label')
         dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
 
         datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         datapointTemplateInternal = datapoints.createNodetemplate()
+        datapointTemplateInternal.defineField(dataCoordinates)
         datapointTemplateInternal.defineField(dataLabel)
         datapointTemplateInternal.defineField(dataElementXi)
 
@@ -361,11 +362,14 @@ def generateBaseMesh(cls, region, options):
 
         # annotation points
         cruxElement = mesh.findElementByIdentifier(cruxElementId)
+        cruxXi = [ 0.0, 0.5, 1.0 ]
+        cache.setMeshLocation(cruxElement, cruxXi)
+        result, cruxCoordinates = coordinates.evaluateReal(cache, 3)
         datapoint = datapoints.createNode(-1, datapointTemplateInternal)
         cache.setNode(datapoint)
-        #dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, vApexInnerx)
+        dataCoordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cruxCoordinates)
         dataLabel.assignString(cache, 'crux')
-        dataElementXi.assignMeshLocation(cache, cruxElement, [ 0.0, 0.5, 1.0 ])
+        dataElementXi.assignMeshLocation(cache, cruxElement, cruxXi)
 
         fm.endChange()
         return annotationGroups

scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py

@@ -188,12 +188,13 @@ def generateBaseMesh(cls, region, options):
         annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ]
 
         # annotation points
-        #dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
+        dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
         dataLabel = getOrCreateLabelField(fm, 'data_label')
         dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
 
         datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         datapointTemplateInternal = datapoints.createNodetemplate()
+        datapointTemplateInternal.defineField(dataCoordinates)
         datapointTemplateInternal.defineField(dataLabel)
         datapointTemplateInternal.defineField(dataElementXi)
 
@@ -246,11 +247,11 @@ def generateBaseMesh(cls, region, options):
         elementSizeUpApex = lengthUpApex/elementsCountUpLVApex
         elementSizeUpRV = (lengthUpRV - 0.5*elementSizeUpApex)/(elementsCountUpRV - 0.5)
         elementSizeUpRVTransition = 0.5*(elementSizeUpApex + elementSizeUpRV)
-        # apex points, noting s1, s2 is x, -y to get out outward s3
-        vApexInnerx = [ 0.0, 0.0, -lvInnerHeight ]
-        vApexInnerd1 = [ elementSizeUpApex, 0.0, 0.0 ]
-        vApexInnerd2 = [ 0.0, -elementSizeUpApex, 0.0, 0.0 ]
-        vApexInnerd3 = [ 0.0, 0.0, -lvApexThickness ]
+        # LV apex points, noting s1, s2 is x, -y to get out outward s3
+        lvApexInnerx = [ 0.0, 0.0, -lvInnerHeight ]
+        lvApexInnerd1 = [ elementSizeUpApex, 0.0, 0.0 ]
+        lvApexInnerd2 = [ 0.0, -elementSizeUpApex, 0.0, 0.0 ]
+        lvApexInnerd3 = [ 0.0, 0.0, -lvApexThickness ]
 
         lvInnerRadiansUp = []
         radiansUp = 0.0
@@ -290,10 +291,10 @@ def generateBaseMesh(cls, region, options):
         elementSizeUpRV = (lengthUpRV - 0.5*elementSizeUpApex)/(elementsCountUpRV - 0.5)
         elementSizeUpRVTransition = 0.5*(elementSizeUpApex + elementSizeUpRV)
         # apex points, noting s1, s2 is x, -y to get out outward s3
-        vApexOuterx = [ 0.0, 0.0, -lvOuterHeight ]
-        vApexOuterd1 = [ elementSizeUpApex, 0.0, 0.0 ]
-        vApexOuterd2 = [ 0.0, -elementSizeUpApex, 0.0, 0.0 ]
-        vApexOuterd3 = [ 0.0, 0.0, -lvApexThickness ]
+        lvApexOuterx = [ 0.0, 0.0, -lvOuterHeight ]
+        lvApexOuterd1 = [ elementSizeUpApex, 0.0, 0.0 ]
+        lvApexOuterd2 = [ 0.0, -elementSizeUpApex, 0.0, 0.0 ]
+        lvApexOuterd3 = [ 0.0, 0.0, -lvApexThickness ]
 
         lvOuterRadiansUp = []
         radiansUp = 0.0
@@ -356,8 +357,8 @@ def generateBaseMesh(cls, region, options):
         for n1 in range(pointsCountAroundOuter):
             cosRadiansAround = math.cos(apexRadiansAround[n1])
             sinRadiansAround = math.sin(apexRadiansAround[n1])
-            apexd2 = [ (vApexOuterd1[c]*cosRadiansAround + vApexOuterd2[c]*sinRadiansAround) for c in range(3) ]
-            nx = [ vApexOuterx ]
+            apexd2 = [ (lvApexOuterd1[c]*cosRadiansAround + lvApexOuterd2[c]*sinRadiansAround) for c in range(3) ]
+            nx = [ lvApexOuterx ]
             nd2 = [ apexd2 ]
             for n2 in range(elementsCountUpLV):
                 nx.append(lvInnerx[n2][n1])
@@ -375,8 +376,8 @@ def generateBaseMesh(cls, region, options):
         for n1 in range(pointsCountAroundOuter):
             cosRadiansAround = math.cos(apexRadiansAround[n1])
             sinRadiansAround = math.sin(apexRadiansAround[n1])
-            apexd2 = [ (vApexOuterd1[c]*cosRadiansAround + vApexOuterd2[c]*sinRadiansAround) for c in range(3) ]
-            nx = [ vApexOuterx ]
+            apexd2 = [ (lvApexOuterd1[c]*cosRadiansAround + lvApexOuterd2[c]*sinRadiansAround) for c in range(3) ]
+            nx = [ lvApexOuterx ]
             nd2 = [ apexd2 ]
             for n2 in range(elementsCountUpLV):
                 nx.append(vOuterx[n2][n1])
@@ -555,10 +556,10 @@ def generateBaseMesh(cls, region, options):
 
         node = nodes.createNode(nodeIdentifier, nodetemplateApex)
         cache.setNode(node)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, vApexInnerx)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, vApexInnerd1)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, vApexInnerd2)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, vApexInnerd3)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lvApexInnerx)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lvApexInnerd1)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lvApexInnerd2)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lvApexInnerd3)
         apexNodeId[1] = nodeIdentifier
         nodeIdentifier = nodeIdentifier + 1
 
@@ -610,10 +611,10 @@ def generateBaseMesh(cls, region, options):
 
         node = nodes.createNode(nodeIdentifier, nodetemplateApex)
         cache.setNode(node)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, vApexOuterx)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, vApexOuterd1)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, vApexOuterd2)
-        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, vApexOuterd3)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lvApexOuterx)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lvApexOuterd1)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lvApexOuterd2)
+        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lvApexOuterd3)
         apexNodeId[1] = nodeIdentifier
         nodeIdentifier = nodeIdentifier + 1
 
@@ -973,10 +974,12 @@ def generateBaseMesh(cls, region, options):
         element1 = mesh.findElementByIdentifier(1)
         datapoint = datapoints.createNode(-1, datapointTemplateInternal)
         cache.setNode(datapoint)
+        dataCoordinates.assignReal(cache, lvApexInnerx)
         dataLabel.assignString(cache, 'apex endo')
         dataElementXi.assignMeshLocation(cache, element1, [ 0.0, 0.0, 0.0 ])
         datapoint = datapoints.createNode(-1, datapointTemplateInternal)
         cache.setNode(datapoint)
+        dataCoordinates.assignReal(cache, lvApexOuterx)
         dataLabel.assignString(cache, 'apex epi')
         dataElementXi.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ])
 

scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py

@@ -13,13 +13,13 @@
 from scaffoldmaker.utils.geometry import *
 from scaffoldmaker.utils.interpolation import *
 from scaffoldmaker.utils.zinc_utils import *
-from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
 from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
 from scaffoldmaker.utils.meshrefinement import MeshRefinement
 import scaffoldmaker.utils.vector as vector
 from opencmiss.zinc.element import Element, Elementbasis
 from opencmiss.zinc.field import Field
 from opencmiss.zinc.node import Node
+from opencmiss.zinc.result import RESULT_OK as ZINC_OK
 
 class MeshType_3d_heartventriclesbase1(object):
     '''
@@ -223,7 +223,7 @@ def generateBaseMesh(cls, region, options):
         # annotation points
         dataCoordinates = getOrCreateCoordinateField(fm, 'data_coordinates')
         dataLabel = getOrCreateLabelField(fm, 'data_label')
-        #dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
+        dataElementXi = getOrCreateElementXiField(fm, 'data_element_xi')
 
         datapoints = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         datapointTemplateExternal = datapoints.createNodetemplate()
@@ -259,10 +259,26 @@ def generateBaseMesh(cls, region, options):
         ventriclesOffset = fm.createFieldConstant([ vTranslationx, vTranslationy, -(fibrousRingThickness + baseHeight + baseThickness)])
         newCoordinates = fm.createFieldAdd(fm.createFieldMatrixMultiply(3, rotationMatrix, coordinates), ventriclesOffset)
         fieldassignment = coordinates.createFieldassignment(newCoordinates)
+        fieldassignment.setNodeset(nodes)
+        fieldassignment.assign()
+        # also transform data point coordinates and/or re-evaluate from embedded locations
+        newCoordinates = fm.createFieldAdd(fm.createFieldMatrixMultiply(3, rotationMatrix, dataCoordinates), ventriclesOffset)
+        fieldassignment = dataCoordinates.createFieldassignment(newCoordinates)
+        fieldassignment.setNodeset(datapoints)
         fieldassignment.assign()
         fieldassignment = None
         newCoordinates = None
         ventriclesOffset = None
+        dataHostCoordinates = fm.createFieldEmbedded(coordinates, dataElementXi)
+        iter = datapoints.createNodeiterator()
+        datapoint = iter.next()
+        while datapoint.isValid():
+            cache.setNode(datapoint)
+            result, datax = dataHostCoordinates.evaluateReal(cache, 3)
+            if result == ZINC_OK:
+                dataCoordinates.assignReal(cache, datax)
+            datapoint = iter.next()
+        dataHostCoordinates = None
 
         # discover ventricles top LV inner, RV inner, V Outer nodes, coordinates and derivatives
         startLVInnerNodeId = 2 + (elementsCountUpLV - 1)*elementsCountAroundLV
@@ -548,7 +564,7 @@ def generateBaseMesh(cls, region, options):
         fd1 = [ (rvOutletOuterd1[2][c] + rvOutletd2[c]) for c in range(3) ]
         pd1 = smoothCubicHermiteDerivativesLine([ ravOuterx[0][ravsvcn2], mx, rvOutletOuterx[1] ], [ sd1, zero, fd1 ],
             fixStartDerivative=True, fixEndDerivative = True)
-        pd2 = smoothCubicHermiteDerivativesLine([ mx, lvOutletOuterx[nf] ], [ md2, svcid2 ], fixStartDirection=True)
+        pd2 = smoothCubicHermiteDerivativesLine([ mx, lvOutletOuterx[nf] ], [ md2, svcid2 ], fixStartDirection=True)  # , instrument=True)
         svcox = [ mx[0], mx[1], mx[2] ]  # list components to avoid reference bug
         svcod1 = pd1[1]
         svcod2 = pd2[0]

scaffoldmaker/scaffolds.py

@@ -11,6 +11,7 @@
 from scaffoldmaker.meshtypes.meshtype_3d_centrallinetube1 import MeshType_3d_centrallinetube1
 from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1
 from scaffoldmaker.meshtypes.meshtype_3d_heart2 import MeshType_3d_heart2
+from scaffoldmaker.meshtypes.meshtype_3d_heartarterialroot1 import MeshType_3d_heartarterialroot1
 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
@@ -38,6 +39,7 @@ def __init__(self):
             MeshType_3d_centrallinetube1,
             MeshType_3d_heart1,
             MeshType_3d_heart2,
+            MeshType_3d_heartarterialroot1,
             MeshType_3d_heartatria1,
             MeshType_3d_heartatria2,
             MeshType_3d_heartventricles1,

scaffoldmaker/utils/interpolation.py

@@ -397,7 +397,7 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
         fixAllDirections = False,
         fixStartDerivative = False, fixEndDerivative = False,
         fixStartDirection = False, fixEndDirection = False,
-        magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN):
+        magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN, instrument=False):
     """
     Modifies derivatives nd1 to be smoothly varying and near arc length.
     Values are treated as being in a line.
@@ -424,22 +424,19 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
     componentsCount = len(nx[0])
     componentRange = range(componentsCount)
     tol = 1.0E-6
-    lastArcLengths = None
+    # special case where equal derivatives at each end are sought
+    equalDerivatives = (elementsCount == 1) and not (fixStartDerivative or fixEndDerivative)
+    if instrument:
+        print('iter 0', md1)
     for iter in range(100):
+        lastmd1 = copy.copy(md1)
         arcLengths = [ getCubicHermiteArcLength(nx[e], md1[e], nx[e + 1], md1[e + 1]) for e in range(elementsCount) ]
-        if lastArcLengths:
-            for n in range(elementsCount):
-                if (arcLengths[n] - lastArcLengths[n]) > (tol*arcLengths[n]):
-                    break
-            else:
-                #print('smoothCubicHermiteDerivativesLine converged after iter:',iter)
-                return md1
         # start
         if not fixStartDerivative:
             if fixAllDirections or fixStartDirection:
-                md1[0] = vector.setMagnitude(md1[0], 2.0*arcLengths[0] - vector.magnitude(md1[1]))
+                md1[0] = vector.setMagnitude(lastmd1[0], 2.0*arcLengths[0] - vector.magnitude(lastmd1[1]))
             else:
-                md1[0] = interpolateLagrangeHermiteDerivative(nx[0], nx[1], md1[1], 0.0)
+                md1[0] = interpolateLagrangeHermiteDerivative(nx[0], nx[1], lastmd1[1], 0.0)
         # middle
         for n in range(1, nodesCount - 1):
             nm = n - 1
@@ -461,16 +458,34 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
         # end
         if not fixEndDerivative:
             if fixAllDirections or fixEndDirection:
-                md1[-1] = vector.setMagnitude(md1[-1], 2.0*arcLengths[-1] - vector.magnitude(md1[-2]))
+                md1[-1] = vector.setMagnitude(lastmd1[-1], 2.0*arcLengths[-1] - vector.magnitude(lastmd1[-2]))
             else:
-                md1[-1] = interpolateHermiteLagrangeDerivative(nx[-2], md1[-2], nx[-1], 1.0)
-        lastArcLengths = arcLengths
+                md1[-1] = interpolateHermiteLagrangeDerivative(nx[-2], lastmd1[-2], nx[-1], 1.0)
+        if equalDerivatives:
+            mag = getCubicHermiteArcLength(nx[0], md1[0], nx[1], md1[1])
+            md1[0] = vector.setMagnitude(md1[0], mag)
+            md1[1] = vector.setMagnitude(md1[1], mag)
+        if instrument:
+            print('iter', iter + 1, md1)
+        dtol = tol*sum(arcLengths)/len(arcLengths)
+        for n in range(nodesCount):
+            for c in componentRange:
+                if math.fabs(md1[n][c] - lastmd1[n][c]) > dtol:
+                    break
+            else:
+                continue
+            break
+        else:
+            if instrument:
+                print('smoothCubicHermiteDerivativesLine converged after iter:',iter)
+            return md1
+
     print('smoothCubicHermiteDerivativesLine max iters reached:',iter)
     return md1
 
 def smoothCubicHermiteDerivativesLoop(nx, nd1,
         fixAllDirections = False,
-        magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN):
+        magnitudeScalingMode = DerivativeScalingMode.ARITHMETIC_MEAN, instrument=False):
     """
     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.
@@ -493,16 +508,11 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1,
     componentsCount = len(nx[0])
     componentRange = range(componentsCount)
     tol = 1.0E-6
-    lastArcLengths = None
+    if instrument:
+        print('iter 0', md1)
     for iter in range(100):
+        lastmd1 = copy.copy(md1)
         arcLengths = [ getCubicHermiteArcLength(nx[e], md1[e], nx[(e + 1)%elementsCount], md1[(e + 1)%elementsCount]) for e in range(elementsCount) ]
-        if lastArcLengths:
-            for n in range(nodesCount):
-                if (arcLengths[n] - lastArcLengths[n]) > (tol*arcLengths[n]):
-                    break
-            else:
-                #print('smoothCubicHermiteDerivativesLoop converged after iter:',iter)
-                return md1
         for n in range(nodesCount):
             nm = n - 1
             if not fixAllDirections:
@@ -520,7 +530,21 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1,
             else: # harmonicMeanMagnitude
                 mag = 2.0/(1.0/arcLengths[nm] + 1.0/arcLengths[n])
             md1[n] = vector.setMagnitude(md1[n], mag)
-        lastArcLengths = arcLengths
+        if instrument:
+            print('iter', iter + 1, md1)
+        dtol = tol*sum(arcLengths)/len(arcLengths)
+        for n in range(nodesCount):
+            for c in componentRange:
+                if math.fabs(md1[n][c] - lastmd1[n][c]) > dtol:
+                    break
+            else:
+                continue
+            break
+        else:
+            if instrument:
+                print('smoothCubicHermiteDerivativesLoop converged after iter:',iter)
+            return md1
+
     print('smoothCubicHermiteDerivativesLoop max iters reached:',iter)
     return md1
 

scaffoldmaker/utils/meshrefinement.py

@@ -96,20 +96,28 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n
         nids = []
         xi = [ 0.0, 0.0, 0.0 ]
         for k in range(numberInXi3 + 1):
+            kExterior = (k == 0) or (k == numberInXi3)
             xi[2] = k/numberInXi3
             for j in range(numberInXi2 + 1):
+                jExterior = kExterior or (j == 0) or (j == numberInXi2)
                 xi[1] = j/numberInXi2
                 for i in range(numberInXi1 + 1):
+                    iExterior = jExterior or (i == 0) or (i == numberInXi1)
                     xi[0] = i/numberInXi1
                     self._sourceCache.setMeshLocation(sourceElement, xi)
                     result, x = self._sourceCoordinates.evaluateReal(self._sourceCache, 3)
-                    nodeId = self._octree.findObjectByCoordinates(x)
+                    # only exterior points are ever common:
+                    if iExterior:
+                        nodeId = self._octree.findObjectByCoordinates(x)
+                    else:
+                        nodeId = None
                     if nodeId is None:
                         node = self._targetNodes.createNode(self._nodeIdentifier, self._nodetemplate)
                         self._targetCache.setNode(node)
                         result = self._targetCoordinates.setNodeParameters(self._targetCache, -1, Node.VALUE_LABEL_VALUE, 1, x)
                         nodeId = self._nodeIdentifier
-                        self._octree.addObjectAtCoordinates(x, nodeId)
+                        if iExterior:
+                            self._octree.addObjectAtCoordinates(x, nodeId)
                         self._nodeIdentifier += 1
                     nids.append(nodeId)
         # create elements