Bifurcation tree start

src/scaffoldmaker/meshtypes/meshtype_1d_bifurcationtree1.py

@@ -134,6 +134,14 @@ def addChild(self, childTreeNode, d1=None, r=None):
             self._r.append(r)
         self._children.append(childTreeNode)
 
+    def getChild(self, childIndex):
+        '''
+        :param childIndex: Index starting at 0 into child tree node list.
+        Return x1, d1, r1, x2, d2, r2
+        '''
+        assert 0 <= childIndex < len(self._children)
+        return self._children[childIndex]
+
     def getChildCurve(self, childIndex):
         '''
         :param childIndex: Index starting at 0 into child tree node list.
@@ -223,7 +231,7 @@ def generateZincModel(self, region, nextNodeIdentifier=1, nextElementIdentifier=
             self._radius = findOrCreateFieldFiniteElement(self._fieldmodule, "radius", components_count=1, managed=True)
             self._fieldcache = self._fieldmodule.createFieldcache()
             parentNode = None
-            self._generateZincModelTree(self._rootNode, parentNode, nextNodeIdentifier, nextElementIdentifier)
+            nextNodeIdentifier, nextElementIdentifier = self._generateZincModelTree(self._rootNode, parentNode, nextNodeIdentifier, nextElementIdentifier)
         return nextNodeIdentifier, nextElementIdentifier
 
     def _getZincNodetemplate(self, d1VersionsCount, rVersionsCount):

src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcationtree1.py

@@ -0,0 +1,149 @@
+"""
+Generates a 2-D tube bifurcation mesh.
+"""
+
+from __future__ import division
+import math
+from opencmiss.utils.maths.vectorops import cross, mult, normalize, sub
+from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates
+from opencmiss.zinc.field import Field
+from opencmiss.zinc.node import Node
+from scaffoldmaker.meshtypes.meshtype_1d_bifurcationtree1 import MeshType_1d_bifurcationtree1
+from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
+from scaffoldmaker.scaffoldpackage import ScaffoldPackage
+from scaffoldmaker.utils.bifurcation import get_curve_circle_points, \
+    make_tube_bifurcation_points, make_tube_bifurcation_elements_2d
+from scaffoldmaker.utils.geometry import createCirclePoints
+from scaffoldmaker.utils.interpolation import getCubicHermiteArcLength
+
+
+class MeshType_2d_tubebifurcationtree1(Scaffold_base):
+    '''
+    Generates a 2-D tube bifurcation tree mesh.
+    '''
+
+    @staticmethod
+    def getName():
+        return '2D Tube Bifurcation Tree 1'
+
+    @staticmethod
+    def getParameterSetNames():
+        return ['Default']
+
+    @staticmethod
+    def getDefaultOptions(parameterSetName='Default'):
+        options = {}
+        options['Bifurcation tree'] = ScaffoldPackage(MeshType_1d_bifurcationtree1)
+        options['Number of elements around root'] = 6
+        options['Maximum element length'] = 0.5
+        return options
+
+    @staticmethod
+    def getOrderedOptionNames():
+        return [
+            'Bifurcation tree',
+            'Number of elements around root',
+            'Maximum element length']
+
+    @classmethod
+    def getOptionValidScaffoldTypes(cls, optionName):
+        if optionName == 'Bifurcation tree':
+            return [ MeshType_1d_bifurcationtree1 ]
+        return []
+
+    @classmethod
+    def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType):
+        assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), cls.__name__ + '.getOptionScaffoldTypeParameterSetNames.  ' + \
+            'Invalid option \'' + optionName + '\' scaffold type ' + scaffoldType.getName()
+        return scaffoldType.getParameterSetNames()
+
+    @classmethod
+    def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None):
+        '''
+        :param parameterSetName:  Name of valid parameter set for option Scaffold, or None for default.
+        :return: ScaffoldPackage.
+        '''
+        if parameterSetName:
+            assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \
+                'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName()
+        if optionName == 'Bifurcation tree':
+            if not parameterSetName:
+                parameterSetName = MeshType_1d_bifurcationtree1.getParameterSetNames()[0]
+            return ScaffoldPackage(MeshType_1d_bifurcationtree1, defaultParameterSetName=parameterSetName)
+        assert False, cls.__name__ + '.getOptionScaffoldPackage:  Option ' + optionName + ' is not a scaffold'
+
+    @staticmethod
+    def checkOptions(options):
+        '''
+        :return:  True if dependent options changed, otherwise False. This
+        happens where two or more options must change together to be valid.
+        '''
+        dependentChanges = False
+        for key in [
+            'Number of elements around root']:
+            if options[key] < 4:
+                options[key] = 4
+        if options['Maximum element length'] <= 0.0:
+            options['Maximum element length'] = 1.0
+        return dependentChanges
+
+    @classmethod
+    def generateBaseMesh(cls, region, options):
+        """
+        Generate the base bicubic Hermite mesh.
+        :param region: Zinc region to define model in. Must be empty.
+        :param options: Dict containing options. See getDefaultOptions().
+        :return: list of AnnotationGroup
+        """
+        bifurcationTreeScaffold = options['Bifurcation tree']
+        maxElementLength = options['Maximum element length']
+        elementsCountAroundRoot = options['Number of elements around root']
+        useCrossDerivatives = False
+
+        fm = region.getFieldmodule()
+        coordinates = findOrCreateFieldCoordinates(fm)
+        cache = fm.createFieldcache()
+
+        ##############
+        # Create nodes
+        ##############
+
+        bifurcationTree = bifurcationTreeScaffold.getScaffoldType().generateBifurcationTree(bifurcationTreeScaffold.getScaffoldSettings())
+        nodeIdentifier = bifurcationTree.generateZincModel(region)[0]
+
+        nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
+        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)
+
+        rootNode = bifurcationTree.getRootNode()
+        child1 = rootNode.getChild(0)
+
+        x1, xd1, r1, x2, xd2, r2 = child1.getChildCurve(0)
+        rd1 = rd2 = (r2 - r1)
+
+        curveLength = getCubicHermiteArcLength(x1, xd1, x2, xd2)
+        elementsCount = max(2, math.ceil(curveLength/maxElementLength))
+        elementLength = curveLength/elementsCount
+        side = [ 1.0, 0.0, 0.0 ]
+        for i in range(elementsCount + 1):
+            xi = i/elementsCount
+            x, d1, d2 = get_curve_circle_points(x1, xd1, x2, xd2, r1, rd1, r2, rd2, xi, elementLength, side, elementsCountAroundRoot)
+            for n in range(elementsCountAroundRoot):
+                #print('node',nodeIdentifier,i,n,x[n],d1[n],d2[n])
+                node = nodes.createNode(nodeIdentifier, nodetemplate)
+                cache.setNode(node)
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x [n])
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1[n])
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2[n])
+                nodeIdentifier = nodeIdentifier + 1
+
+        #################
+        # Create elements
+        #################
+
+        elementIdentifier = 1
+
+        return []

src/scaffoldmaker/scaffolds.py

@@ -10,6 +10,7 @@
 from scaffoldmaker.meshtypes.meshtype_2d_sphere1 import MeshType_2d_sphere1
 from scaffoldmaker.meshtypes.meshtype_2d_tube1 import MeshType_2d_tube1
 from scaffoldmaker.meshtypes.meshtype_2d_tubebifurcation1 import MeshType_2d_tubebifurcation1
+from scaffoldmaker.meshtypes.meshtype_2d_tubebifurcationtree1 import MeshType_2d_tubebifurcationtree1
 from scaffoldmaker.meshtypes.meshtype_3d_bladder1 import MeshType_3d_bladder1
 from scaffoldmaker.meshtypes.meshtype_3d_bladderurethra1 import MeshType_3d_bladderurethra1
 from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1
@@ -54,6 +55,7 @@ def __init__(self):
             MeshType_2d_sphere1,
             MeshType_2d_tube1,
             MeshType_2d_tubebifurcation1,
+            #MeshType_2d_tubebifurcationtree1,
             MeshType_3d_bladder1,
             MeshType_3d_bladderurethra1,
             MeshType_3d_box1,

src/scaffoldmaker/utils/bifurcation.py

@@ -3,15 +3,54 @@
 """
 
 from __future__ import division
-from opencmiss.utils.maths.vectorops import add, cross, mult, normalize, sub
+from opencmiss.utils.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub
 from opencmiss.zinc.element import Element, Elementbasis
 from opencmiss.zinc.field import Field
 from opencmiss.zinc.node import Node
 from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds
+from scaffoldmaker.utils.geometry import createCirclePoints
 from scaffoldmaker.utils.interpolation import DerivativeScalingMode, interpolateCubicHermite, interpolateCubicHermiteDerivative, \
-    smoothCubicHermiteDerivativesLine, interpolateLagrangeHermiteDerivative
+    interpolateCubicHermiteSecondDerivative, smoothCubicHermiteDerivativesLine, interpolateLagrangeHermiteDerivative
 
 
+def get_curve_circle_points(x1, xd1, x2, xd2, r1, rd1, r2, rd2, xi, dmag, side, elementsCountAround):
+    '''
+    :param dmag: Magnitude of derivative on curve.
+    :param side: Vector in side direction of first node around.
+    Need not be unit or exactly normal to curve at xi.
+    :return: x[], d1[] around, d2[] along
+    '''
+    cx = interpolateCubicHermite(x1, xd1, x2, xd2, xi)
+    cxd = interpolateCubicHermiteDerivative(x1, xd1, x2, xd2, xi)
+    mag_cxd = magnitude(cxd)
+    cxd2 = interpolateCubicHermiteSecondDerivative(x1, xd1, x2, xd2, xi)
+    mag_cxd2 = magnitude(cxd2)
+    r = interpolateCubicHermite([ r1 ], [ rd1 ], [ r2 ], [ rd2 ], xi)[0]
+    rd = interpolateCubicHermiteDerivative([ r1 ], [ rd1 ], [ r2 ], [ rd2 ], xi)[0]
+    axis1 = normalize(cxd)
+    axis3 = normalize(cross(axis1, side))
+    axis2 = cross(axis3, axis1)
+    x, d1 = createCirclePoints(cx, mult(axis2, r), mult(axis3, r), elementsCountAround)
+    curvatureVector = mult(cxd2, 1.0/(mag_cxd*mag_cxd))
+    d2 = []
+    radialGrowth = rd/(mag_cxd*r)
+    for e in range(elementsCountAround):
+        radialVector = sub(x[e], cx)
+        dmagFinal = dmag*(1.0 - dot(radialVector, curvatureVector))
+        # add curvature and radius change components:
+        d2.append(add(mult(cxd, dmagFinal/mag_cxd), mult(radialVector, dmagFinal*radialGrowth)))
+    return x, d1, d2
+
+
+def track_curve_side_axis(x1, d1, x2, d2, sideStart, xiStart, xiEnd):
+    '''
+    Get side vector normal to curve at xiEnd for smoothest transition from
+    sideStart at xiStart.
+    :param xi, d1, x2, d2: 
+    :param sideStart: Unit vector normal to curve
+    '''
+    pass
+
 def get_bifurcation_triple_point(p1x, p1d, p2x, p2d, p3x, p3d):
     '''
     Get coordinates and derivatives of triple point between p1, p2 and p3 with derivatives.