From 5bfeff79aa9eef41b3d3931abb21c9f062dd6d13 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Thu, 22 Nov 2018 19:20:42 +1300 Subject: [PATCH 1/8] Add initial heart arterial root scaffold --- .../meshtype_3d_heartarterialroot1.py | 278 ++++++++++++++++++ .../meshtype_3d_heartventriclesbase1.py | 1 - scaffoldmaker/scaffolds.py | 2 + 3 files changed, 280 insertions(+), 1 deletion(-) create mode 100644 scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py new file mode 100644 index 00000000..1d408e29 --- /dev/null +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -0,0 +1,278 @@ +""" +Generates a 3-D heart arterial root scaffold with semilunar valve, +for attaching to a 6-element-around bicubic-linear orifice. +""" + +from __future__ import division +import math +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findAnnotationGroupByName +from scaffoldmaker.utils.eft_utils import * +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 + +class MeshType_3d_heartarterialroot1(object): + ''' + Generates a 3-D heart arterial root scaffold with semilunar valve, + for attaching to a 6-element-around bicubic-linear orifice. + ''' + + @staticmethod + def getName(): + return '3D Heart Arterial Root 1' + + @staticmethod + def getDefaultOptions(): + return { + 'Unit scale' : 1.0, + 'Outer height' : 0.5, + 'Inner depth' : 0.2, + 'Cusp height' : 0.5, + 'Inner diameter': 1.0, + 'Sinus radial displacement': 0.05, + 'Wall thickness': 0.1, + 'Cusp thickness' : 0.04, + 'Aortic not pulmonary' : True, + 'Refine' : False, + 'Refine number of elements surface' : 4, + 'Use cross derivatives' : False + } + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Unit scale', + 'Outer height', + 'Inner depth', + 'Cusp height', + 'Inner diameter', + 'Sinus radial displacement', + 'Wall thickness', + 'Cusp thickness', + 'Aortic not pulmonary', + 'Refine', + 'Refine number of elements surface' + ] + + @staticmethod + def checkOptions(options): + ''' + :return: True if dependent options changed, otherwise False. + ''' + dependentChanges = False + for key in [ + 'Unit scale', + 'Outer height', + 'Inner depth', + 'Cusp height', + 'Inner diameter', + 'Sinus radial displacement', + 'Wall thickness', + 'Cusp thickness']: + if options[key] < 0.0: + options[key] = 0.0 + for key in [ + 'Refine number of elements surface']: + if options[key] < 1: + options[key] = 1 + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSide1=[ 0.0, -1.0, 0.0 ], axisUp=[ 0.0, 0.0, 1.0]): + """ + Generate the base bicubic-linear Hermite mesh. See also generateMesh(). + Optional extra parameters allow centre and axes to be set. + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :param baseCentre: Centre of valve on ventriculo-arterial junction. + :param axisSide: Unit vector in first side direction where angle around starts. + :param axisUp: Unit vector in outflow direction of valve. + :return: list of AnnotationGroup + """ + unitScale = options['Unit scale'] + outerHeight = unitScale*options['Outer height'] + innerDepth = unitScale*options['Inner depth'] + cuspHeight = unitScale*options['Cusp height'] + innerRadius = unitScale*0.5*options['Inner diameter'] + sinusRadialDisplacement = unitScale*options['Sinus radial displacement'] + wallThickness = unitScale*options['Wall thickness'] + cuspThickness = unitScale*options['Cusp thickness'] + aorticNotPulmonary = options['Aortic not pulmonary'] + useCrossDerivatives = False + + fm = region.getFieldmodule() + fm.beginChange() + coordinates = getOrCreateCoordinateField(fm) + cache = fm.createFieldcache() + + if aorticNotPulmonary: + arterialRootGroup = AnnotationGroup(region, 'root of aorta', FMANumber = 3740, lyphID = 'Lyph ID unknown') + cuspGroups = [ + AnnotationGroup(region, 'right cusp of aortic valve', FMANumber = 7252, lyphID = 'Lyph ID unknown'), + AnnotationGroup(region, 'left cusp of aortic valve', FMANumber = 7251, lyphID = 'Lyph ID unknown'), + AnnotationGroup(region, 'posterior cusp of aortic valve', FMANumber = 7253, lyphID = 'Lyph ID unknown') ] + else: + arterialRootGroup = AnnotationGroup(region, 'root of pulmonary trunk', FMANumber = 8612, lyphID = 'Lyph ID unknown') + cuspGroups = [ + AnnotationGroup(region, 'right cusp of pulmonary valve', FMANumber = 7250, lyphID = 'Lyph ID unknown'), + AnnotationGroup(region, 'anterior cusp of pulmonary valve', FMANumber = 7249, lyphID = 'Lyph ID unknown'), + AnnotationGroup(region, 'left cusp of pulmonary valve', FMANumber = 7247, lyphID = 'Lyph ID unknown') ] + + allGroups = [ arterialRootGroup ] # groups that all elements in scaffold will go in + annotationGroups = allGroups + cuspGroups + + ################# + # Create nodes + ################# + + 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) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + # most nodes in this scaffold do not have a DS3 derivative + nodetemplateLinearS3 = nodes.createNodetemplate() + nodetemplateLinearS3.defineField(coordinates) + nodetemplateLinearS3.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplateLinearS3.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplateLinearS3.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + + nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) + + elementsCountAround = 6 + axisSide2 = vector.crossproduct3(axisUp, axisSide1) + startRadians = (math.pi/3.0) if aorticNotPulmonary else 0.0 + outerRadius = innerRadius + wallThickness + + # lower points + lowerx, lowerd1 = getArterialRootLowerPoints(baseCentre, axisSide1, axisSide2, axisUp, innerRadius, outerRadius, innerDepth, sinusRadialDisplacement, startRadians, elementsCountAround) + + # upper points + topCentre = [ (baseCentre[c] + axisUp[c]*outerHeight) for c in range(3) ] + ix, id1 = createCirclePoints(topCentre, + [ axisSide1[c]*innerRadius for c in range(3) ], [ axisSide2[c]*innerRadius for c in range(3) ], + elementsCountAround, startRadians) + ox, od1 = createCirclePoints(topCentre, + [ axisSide1[c]*outerRadius for c in range(3) ], [ axisSide2[c]*outerRadius for c in range(3) ], + elementsCountAround, startRadians) + upperx, upperd1 = [ ix, ox ], [ id1, od1 ] + + for n3 in range(2): + for n1 in range(elementsCountAround): + node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lowerx [n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lowerd1[n3][n1]) + #coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lowerd2[n3][n1]) + nodeIdentifier += 1 + + for n3 in range(2): + for n1 in range(elementsCountAround): + node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, upperx [n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, upperd1[n3][n1]) + #coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, upperd2[n3][n1]) + nodeIdentifier += 1 + + + ################# + # Create elements + ################# + + mesh = fm.findMeshByDimension(3) + + allMeshGroups = [ allGroup.getMeshGroup(mesh) for allGroup in allGroups ] + cuspMeshGroups = [ cuspGroup.getMeshGroup(mesh) for cuspGroup in cuspGroups ] + + bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) + eft = bicubichermitelinear.createEftNoCrossDerivatives() + + tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + + elementIdentifier = startElementIdentifier = getMaximumElementIdentifier(mesh) + 1 + + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + scalefactors5hanging = [ -1.0, 0.5, 0.25, 0.125, 0.75 ] + + fm.endChange() + return annotationGroups + + @classmethod + def refineMesh(cls, meshrefinement, options): + """ + Refine source mesh into separate region, with change of basis. + :param meshrefinement: MeshRefinement, which knows source and target region. + :param options: Dict containing options. See getDefaultOptions(). + """ + assert isinstance(meshrefinement, MeshRefinement) + refineElementsCountSurface = options['Refine number of elements surface'] + MeshType_3d_heartventricles1.refineMesh(meshrefinement, options) + numberInXi1 = refineElementsCountSurface + numberInXi2 = refineElementsCountSurface + numberInXi3 = 1 + for i in range(0): # eventually: (24): # fixed number of elements + element = meshrefinement._sourceElementiterator.next() + meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3) + + @classmethod + def generateMesh(cls, region, options): + """ + Generate base or refined mesh. + :param region: Zinc region to create mesh in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup for mesh. + """ + if not options['Refine']: + return cls.generateBaseMesh(region, options) + baseRegion = region.createRegion() + baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) + meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) + cls.refineMesh(meshrefinement, options) + return meshrefinement.getAnnotationGroups() + + +def getArterialRootLowerPoints(baseCentre, axisSide1, axisSide2, axisUp, innerRadius, outerRadius, innerDepth, sinusRadialDisplacement, + startRadians, elementsCountAround = 6): + ''' + :param startRadians: Angle from axisSide1 toward axisSide2 where first cusp begins. + :return: Coordinates, derivatives x[2][6], d1[2][6] for inner:outer and 6 around. + ''' + assert elementsCountAround == 6, 'getArterialRootLowerPoints. Only supports 6 elements around' + ix, id1 = createCirclePoints([ (baseCentre[c] - axisUp[c]*innerDepth) for c in range(3) ], + [ axisSide1[c]*innerRadius for c in range(3) ], [ axisSide2[c]*innerRadius for c in range(3) ], + elementsCountAround, startRadians) + + # every second outer points is displaced by sinusRadialDisplacement to make sinus dilatation + ox = [] + od1 = [] + outerRadiusPlus = outerRadius + sinusRadialDisplacement + radiansPerElementAround = 2.0*math.pi/elementsCountAround + radiansAround = startRadians + for n in range(elementsCountAround): + radius = outerRadiusPlus if (n % 2) else outerRadius + rcosRadiansAround = radius*math.cos(radiansAround) + rsinRadiansAround = radius*math.sin(radiansAround) + ox.append([ (baseCentre[c] + rcosRadiansAround*axisSide1[c] + rsinRadiansAround*axisSide2[c]) for c in range(3) ]) + od1.append([ radiansPerElementAround*(-rsinRadiansAround*axisSide1[c] + rcosRadiansAround*axisSide2[c]) for c in range(3) ]) + radiansAround += radiansPerElementAround + # smooth to get derivative in sinus + pd1 = smoothCubicHermiteDerivativesLine(ox[0:2], od1[0:2], fixStartDerivative=True, fixEndDirection=True) + od1[1] = pd1[1] + magSinus = vector.magnitude(pd1[1]) + for n in range(3, elementsCountAround, 2): + od1[n] = vector.setMagnitude(od1[n], magSinus) + + return [ ix, ox ], [ id1, od1 ] diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 89020c1f..4d184436 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -13,7 +13,6 @@ 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 diff --git a/scaffoldmaker/scaffolds.py b/scaffoldmaker/scaffolds.py index c110fae9..92c469be 100644 --- a/scaffoldmaker/scaffolds.py +++ b/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, From d353f1c309d888d4cf1ee465fe2e882536f02a1d Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 23 Nov 2018 09:09:46 +1300 Subject: [PATCH 2/8] Fix default name for element xi field --- scaffoldmaker/utils/zinc_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scaffoldmaker/utils/zinc_utils.py b/scaffoldmaker/utils/zinc_utils.py index 2c135b55..38eb0059 100644 --- a/scaffoldmaker/utils/zinc_utils.py +++ b/scaffoldmaker/utils/zinc_utils.py @@ -38,7 +38,7 @@ def getOrCreateCoordinateField(fieldmodule, name='coordinates', componentsCount= fieldmodule.endChange() return coordinates -def getOrCreateElementXiField(fieldmodule, name='label', mesh=None): +def getOrCreateElementXiField(fieldmodule, name='element_xi', mesh=None): ''' Finds or creates a stored mesh location field for storing locations in the supplied mesh e.g. for defining on annotation points with mesh locations. From 62ebea7ee961e2335e1a192d0b9e965909e6bee2 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sun, 25 Nov 2018 15:29:54 +1300 Subject: [PATCH 3/8] Fix convergence of derivative smoothing --- .../meshtype_3d_heartventriclesbase1.py | 2 +- scaffoldmaker/utils/interpolation.py | 72 ++++++++++++------- 2 files changed, 49 insertions(+), 25 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 4d184436..862a368d 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -547,7 +547,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] diff --git a/scaffoldmaker/utils/interpolation.py b/scaffoldmaker/utils/interpolation.py index 7c85c89e..411ff106 100644 --- a/scaffoldmaker/utils/interpolation.py +++ b/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 From 73e13c522fdaddca1e3aaa1e08f5fbbbc9a29db4 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sun, 25 Nov 2018 15:50:42 +1300 Subject: [PATCH 4/8] Improve performance of refinement Avoid searching for and adding interior nodes in octree. --- scaffoldmaker/utils/meshrefinement.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/scaffoldmaker/utils/meshrefinement.py b/scaffoldmaker/utils/meshrefinement.py index 1178699c..a13d84c4 100644 --- a/scaffoldmaker/utils/meshrefinement.py +++ b/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 From 723cd76acfb2e48a66a2e223f306f352d7069412 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sun, 25 Nov 2018 19:00:38 +1300 Subject: [PATCH 5/8] Add more arterial root nodes, derivatives --- .../meshtype_3d_heartarterialroot1.py | 155 ++++++++++++++---- 1 file changed, 121 insertions(+), 34 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index 1d408e29..ad6bc457 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -150,42 +150,134 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) elementsCountAround = 6 + radiansPerElementAround = 2.0*math.pi/elementsCountAround axisSide2 = vector.crossproduct3(axisUp, axisSide1) - startRadians = (math.pi/3.0) if aorticNotPulmonary else 0.0 outerRadius = innerRadius + wallThickness + cuspLowerLength2 = 0.5*getApproximateEllipsePerimeter(innerRadius, cuspHeight) + cuspLowerInnerArcLength = cuspLowerLength2*innerRadius/(innerRadius + cuspHeight) + cuspNoduleArcLength = cuspLowerLength2 - cuspLowerInnerArcLength + cuspUpperInnerd3 = interpolateHermiteLagrangeDerivative([ 0.0, 0.0 ], [ innerRadius, 0.0 ], [ innerRadius, outerHeight + innerDepth - cuspHeight ], 1.0) # lower points - lowerx, lowerd1 = getArterialRootLowerPoints(baseCentre, axisSide1, axisSide2, axisUp, innerRadius, outerRadius, innerDepth, sinusRadialDisplacement, startRadians, elementsCountAround) + ix, id1 = createCirclePoints([ (baseCentre[c] - axisUp[c]*innerDepth) for c in range(3) ], + [ axisSide1[c]*innerRadius for c in range(3) ], [ axisSide2[c]*innerRadius for c in range(3) ], + elementsCountAround) + ox, od1 = getSemilunarValveSinusPoints(baseCentre, axisSide1, axisSide2, outerRadius, sinusRadialDisplacement, + startMidCusp=aorticNotPulmonary) + lowerx, lowerd1 = [ ix, ox ], [ id1, od1 ] # upper points topCentre = [ (baseCentre[c] + axisUp[c]*outerHeight) for c in range(3) ] + # twice as many on inner: ix, id1 = createCirclePoints(topCentre, [ axisSide1[c]*innerRadius for c in range(3) ], [ axisSide2[c]*innerRadius for c in range(3) ], - elementsCountAround, startRadians) + elementsCountAround*2) ox, od1 = createCirclePoints(topCentre, [ axisSide1[c]*outerRadius for c in range(3) ], [ axisSide2[c]*outerRadius for c in range(3) ], - elementsCountAround, startRadians) + elementsCountAround) upperx, upperd1 = [ ix, ox ], [ id1, od1 ] + # get lower and upper derivative 2 + zero = [ 0.0, 0.0, 0.0 ] + nMidCusp = 0 if aorticNotPulmonary else 1 + upperd2factor = outerHeight*(2.0/3.0) + upd2 = [ d*upperd2factor for d in axisUp ] + lowerOuterd2 = smoothCubicHermiteDerivativesLine([ lowerx[1][nMidCusp], upperx[1][nMidCusp] ], [ upd2, upd2 ], + fixStartDirection=True, fixEndDerivative=True)[0] + lowerd2 = [ [ upd2 ]*elementsCountAround, [ lowerOuterd2 ]*elementsCountAround ] # lowerd2[0] to be fitted below + upperd2 = [ [ upd2 ]*(elementsCountAround*2), [ upd2 ]*elementsCountAround ] + + # get lower and upper derivative 3 for points which have them + lowerd3 = [ [ None ]*elementsCountAround, None ] + upperd3 = [ [ None ]*(elementsCountAround*2), None ] + for n1 in range(elementsCountAround): + radiansAround = n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + if (n1 % 2) == nMidCusp: + lowerd3[0][n1] = [ cuspLowerInnerArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ] + else: + upperd3[0][n1*2] = [ (cuspUpperInnerd3[0]*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) + cuspUpperInnerd3[1]*axisUp[c]) for c in range(3) ] + + # inner-wall mid sinus points + midDistance = 0.5*(outerHeight - innerDepth) + outerArcLength2 = getCubicHermiteArcLength(lowerx[1][nMidCusp], lowerd2[1][nMidCusp], upperx[1][nMidCusp], upperd2[1][nMidCusp]) + outerArcLength2Part = outerArcLength2*midDistance/outerHeight + mx, md2 = getCubicHermiteCurvesPointAtArcDistance([ lowerx[1][nMidCusp], upperx[1][nMidCusp] ], [ lowerd2[1][nMidCusp], upperd2[1][nMidCusp] ], outerArcLength2Part)[0:2] + mn = vector.setMagnitude(vector.crossproduct3(lowerd1[1][nMidCusp], md2), wallThickness) + mx = [ (mx[c] - mn[c]) for c in range(3) ] + id2 = smoothCubicHermiteDerivativesLine([ lowerx[0][nMidCusp], mx, upperx[0][nMidCusp] ], [ lowerd2[0][nMidCusp], md2, upperd2[0][nMidCusp] ], + fixAllDirections=True, fixEndDerivative=True) + lowerd2[0] = [ id2[0] ]*elementsCountAround + + startRadians = 0.0 if aorticNotPulmonary else math.pi/3.0 + cosStartRadians = math.cos(startRadians) + sinStartRadians = math.sin(startRadians) + axisRadial = [ (cosStartRadians*axisSide1[c] + sinStartRadians*axisSide2[c]) for c in range(3) ] + mr = vector.dotproduct([ (mx[c] - baseCentre[c]) for c in range(3) ], axisRadial) + mc = [ (mx[c] - mr*axisRadial[c]) for c in range(3) ] + md2a = vector.dotproduct(id2[1], axisUp) + md2r = vector.dotproduct(id2[1], axisRadial) + ix, id1 = getSemilunarValveSinusPoints(mc, axisSide1, axisSide2, innerRadius, mr - innerRadius, + startMidCusp=aorticNotPulmonary) + # resample to double number of points, halving derivative size: + jx = [] + jd1 = [] + xi = 0.5 + for n1 in range(elementsCountAround): + np = (n1 + 1)%elementsCountAround + jx .append(interpolateCubicHermite (ix[n1], id1[n1], ix[np], id1[np], xi)) + jd1.append(interpolateCubicHermiteDerivative(ix[n1], id1[n1], ix[np], id1[np], xi)) + sinusx = [] + sinusd1 = [] + for n1 in range(elementsCountAround): + sinusx .append(ix [n1]) + sinusx .append(jx [n1]) + sinusd1.append([ 0.5*d for d in id1[n1]]) + sinusd1.append([ 0.5*d for d in jd1[n1]]) + sinusd2 = [ None ]*(elementsCountAround*2) + for n1 in range(elementsCountAround*2): + radiansAround = 0.5*n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + n1Cusp = (n1 - nMidCusp*2)%4 + rf = 0.0 if (n1Cusp == 2) else (1.0 if (n1Cusp == 0) else 0.5) + sinusd2[n1] = [ (rf*md2r*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) + md2a*axisUp[c]) for c in range(3) ] + + # cusp nodule points + for n3 in range(2): for n1 in range(elementsCountAround): - node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) + node = nodes.createNode(nodeIdentifier, nodetemplate if (lowerd3[n3] and lowerd3[n3][n1]) else nodetemplateLinearS3) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lowerx [n3][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lowerd1[n3][n1]) - #coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lowerd2[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lowerd2[n3][n1]) + if lowerd3[n3] and lowerd3[n3][n1]: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lowerd3[n3][n1]) nodeIdentifier += 1 + for n1 in range(elementsCountAround*2): + if ((n1 - nMidCusp*2)%4) == 2: + continue + node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, sinusx [n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, sinusd1[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, sinusd2[n1]) + nodeIdentifier += 1 + for n3 in range(2): - for n1 in range(elementsCountAround): - node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) + for n1 in range(len(upperx[n3])): + node = nodes.createNode(nodeIdentifier, nodetemplate if (upperd3[n3] and upperd3[n3][n1]) else nodetemplateLinearS3) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, upperx [n3][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, upperd1[n3][n1]) - #coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, upperd2[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, upperd2[n3][n1]) + if upperd3[n3] and upperd3[n3][n1]: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, upperd3[n3][n1]) nodeIdentifier += 1 - ################# # Create elements ################# @@ -244,35 +336,30 @@ def generateMesh(cls, region, options): return meshrefinement.getAnnotationGroups() -def getArterialRootLowerPoints(baseCentre, axisSide1, axisSide2, axisUp, innerRadius, outerRadius, innerDepth, sinusRadialDisplacement, - startRadians, elementsCountAround = 6): +def getSemilunarValveSinusPoints(centre, axisSide1, axisSide2, radius, sinusRadialDisplacement, startMidCusp, elementsCountAround = 6): ''' - :param startRadians: Angle from axisSide1 toward axisSide2 where first cusp begins. - :return: Coordinates, derivatives x[2][6], d1[2][6] for inner:outer and 6 around. + Get points around a circle of radius with every second point displaced radially. + :return: x[], d1[] for 6 points around ''' assert elementsCountAround == 6, 'getArterialRootLowerPoints. Only supports 6 elements around' - ix, id1 = createCirclePoints([ (baseCentre[c] - axisUp[c]*innerDepth) for c in range(3) ], - [ axisSide1[c]*innerRadius for c in range(3) ], [ axisSide2[c]*innerRadius for c in range(3) ], - elementsCountAround, startRadians) - + px = [] + pd1 = [] # every second outer points is displaced by sinusRadialDisplacement to make sinus dilatation - ox = [] - od1 = [] - outerRadiusPlus = outerRadius + sinusRadialDisplacement + nMidCusp = 0 if startMidCusp else 1 + radiusPlus = radius + sinusRadialDisplacement radiansPerElementAround = 2.0*math.pi/elementsCountAround - radiansAround = startRadians for n in range(elementsCountAround): - radius = outerRadiusPlus if (n % 2) else outerRadius - rcosRadiansAround = radius*math.cos(radiansAround) - rsinRadiansAround = radius*math.sin(radiansAround) - ox.append([ (baseCentre[c] + rcosRadiansAround*axisSide1[c] + rsinRadiansAround*axisSide2[c]) for c in range(3) ]) - od1.append([ radiansPerElementAround*(-rsinRadiansAround*axisSide1[c] + rcosRadiansAround*axisSide2[c]) for c in range(3) ]) - radiansAround += radiansPerElementAround + radiansAround = n*radiansPerElementAround + midCusp = (n % 2) == nMidCusp + r = radiusPlus if midCusp else radius + rcosRadiansAround = r*math.cos(radiansAround) + rsinRadiansAround = r*math.sin(radiansAround) + px.append([ (centre[c] + rcosRadiansAround*axisSide1[c] + rsinRadiansAround*axisSide2[c]) for c in range(3) ]) + pd1.append([ radiansPerElementAround*(-rsinRadiansAround*axisSide1[c] + rcosRadiansAround*axisSide2[c]) for c in range(3) ]) # smooth to get derivative in sinus - pd1 = smoothCubicHermiteDerivativesLine(ox[0:2], od1[0:2], fixStartDerivative=True, fixEndDirection=True) - od1[1] = pd1[1] - magSinus = vector.magnitude(pd1[1]) - for n in range(3, elementsCountAround, 2): - od1[n] = vector.setMagnitude(od1[n], magSinus) + sd1 = smoothCubicHermiteDerivativesLine(px[1 - nMidCusp:3 - nMidCusp], pd1[1 - nMidCusp:3 - nMidCusp], fixStartDerivative=True, fixEndDirection=True) + magSinus = vector.magnitude(sd1[1]) + for n in range(nMidCusp, elementsCountAround, 2): + pd1[n] = vector.setMagnitude(pd1[n], magSinus) - return [ ix, ox ], [ id1, od1 ] + return px, pd1 From 0ca760ba26bf09f659d263efb1c264faf9de7cc4 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sun, 25 Nov 2018 20:21:16 +1300 Subject: [PATCH 6/8] Add cusp nodule nodes --- .../meshtype_3d_heartarterialroot1.py | 58 ++++++++++++++++--- 1 file changed, 50 insertions(+), 8 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index ad6bc457..6862687b 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -34,11 +34,11 @@ def getDefaultOptions(): 'Unit scale' : 1.0, 'Outer height' : 0.5, 'Inner depth' : 0.2, - 'Cusp height' : 0.5, + 'Cusp height' : 0.6, 'Inner diameter': 1.0, 'Sinus radial displacement': 0.05, 'Wall thickness': 0.1, - 'Cusp thickness' : 0.04, + 'Cusp thickness' : 0.02, 'Aortic not pulmonary' : True, 'Refine' : False, 'Refine number of elements surface' : 4, @@ -153,10 +153,14 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid radiansPerElementAround = 2.0*math.pi/elementsCountAround axisSide2 = vector.crossproduct3(axisUp, axisSide1) outerRadius = innerRadius + wallThickness - cuspLowerLength2 = 0.5*getApproximateEllipsePerimeter(innerRadius, cuspHeight) - cuspLowerInnerArcLength = cuspLowerLength2*innerRadius/(innerRadius + cuspHeight) - cuspNoduleArcLength = cuspLowerLength2 - cuspLowerInnerArcLength - cuspUpperInnerd3 = interpolateHermiteLagrangeDerivative([ 0.0, 0.0 ], [ innerRadius, 0.0 ], [ innerRadius, outerHeight + innerDepth - cuspHeight ], 1.0) + cuspOuterLength2 = 0.5*getApproximateEllipsePerimeter(innerRadius, cuspHeight) + cuspOuterWallArcLength = cuspOuterLength2*innerRadius/(innerRadius + cuspHeight) + noduleOuterAxialArcLength = cuspOuterLength2 - cuspOuterWallArcLength + noduleOuterRadialArcLength = innerRadius + cuspOuterWalld3 = interpolateHermiteLagrangeDerivative([ 0.0, 0.0 ], [ innerRadius, 0.0 ], [ innerRadius, outerHeight + innerDepth - cuspHeight ], 1.0) + cuspInnerLength2 = 0.5*getApproximateEllipsePerimeter(innerRadius - 2.0*cuspThickness, cuspHeight - cuspThickness) + noduleInnerAxialArcLength = cuspInnerLength2*(cuspHeight - cuspThickness)/(innerRadius + cuspHeight - 3.0*cuspThickness) + noduleInnerRadialArcLength = innerRadius - cuspThickness/math.tan(math.pi/3.0) # lower points ix, id1 = createCirclePoints([ (baseCentre[c] - axisUp[c]*innerDepth) for c in range(3) ], @@ -195,9 +199,9 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid cosRadiansAround = math.cos(radiansAround) sinRadiansAround = math.sin(radiansAround) if (n1 % 2) == nMidCusp: - lowerd3[0][n1] = [ cuspLowerInnerArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ] + lowerd3[0][n1] = [ cuspOuterWallArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ] else: - upperd3[0][n1*2] = [ (cuspUpperInnerd3[0]*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) + cuspUpperInnerd3[1]*axisUp[c]) for c in range(3) ] + upperd3[0][n1*2] = [ (cuspOuterWalld3[0]*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) + cuspOuterWalld3[1]*axisUp[c]) for c in range(3) ] # inner-wall mid sinus points midDistance = 0.5*(outerHeight - innerDepth) @@ -245,6 +249,34 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid sinusd2[n1] = [ (rf*md2r*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) + md2a*axisUp[c]) for c in range(3) ] # cusp nodule points + noduleCentre = [ (baseCentre[c] + axisUp[c]*(cuspHeight - innerDepth)) for c in range(3) ] + nodulex = [ [], [] ] + noduled1 = [ [], [] ] + noduled2 = [ [], [] ] + noduled3 = [ [], [] ] + for i in range(3): + nodulex[0].append(noduleCentre) + n1 = i*2 + nMidCusp + radiansAround = n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + nodulex[1].append([ (noduleCentre[c] + 2.0*cuspThickness*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c])) for c in range(3) ]) + n1 = i*2 - 1 + nMidCusp + radiansAround = n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + noduled1[0].append([ noduleOuterRadialArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ]) + noduled1[1].append(vector.setMagnitude(noduled1[0][i], noduleInnerRadialArcLength)) + n1 = i*2 + 1 + nMidCusp + radiansAround = n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + noduled2[0].append([ noduleOuterRadialArcLength*(cosRadiansAround*axisSide1[c] + sinRadiansAround*axisSide2[c]) for c in range(3) ]) + noduled2[1].append(vector.setMagnitude(noduled2[0][i], noduleInnerRadialArcLength)) + noduled3[0].append([ -noduleOuterAxialArcLength*axisUp[c] for c in range(3) ]) + noduled3[1].append([ -noduleInnerAxialArcLength*axisUp[c] for c in range(3) ]) + + # Create nodes for n3 in range(2): for n1 in range(elementsCountAround): @@ -267,6 +299,16 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, sinusd2[n1]) nodeIdentifier += 1 + for n3 in range(2): + for n1 in range(3): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodulex [n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, noduled1[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, noduled2[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, noduled3[n3][n1]) + nodeIdentifier += 1 + for n3 in range(2): for n1 in range(len(upperx[n3])): node = nodes.createNode(nodeIdentifier, nodetemplate if (upperd3[n3] and upperd3[n3][n1]) else nodetemplateLinearS3) From cbc7d5736a05758e485c1849693d9d0a76e22eeb Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Sun, 25 Nov 2018 23:37:29 +1300 Subject: [PATCH 7/8] Add arterial wall elements --- .../meshtype_3d_heartarterialroot1.py | 87 ++++++++++++++++--- 1 file changed, 76 insertions(+), 11 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py index 6862687b..20997407 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartarterialroot1.py @@ -36,7 +36,7 @@ def getDefaultOptions(): 'Inner depth' : 0.2, 'Cusp height' : 0.6, 'Inner diameter': 1.0, - 'Sinus radial displacement': 0.05, + 'Sinus radial displacement': 0.1, 'Wall thickness': 0.1, 'Cusp thickness' : 0.02, 'Aortic not pulmonary' : True, @@ -278,6 +278,7 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid # Create nodes + lowerNodeId = [ [], [] ] for n3 in range(2): for n1 in range(elementsCountAround): node = nodes.createNode(nodeIdentifier, nodetemplate if (lowerd3[n3] and lowerd3[n3][n1]) else nodetemplateLinearS3) @@ -287,18 +288,23 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lowerd2[n3][n1]) if lowerd3[n3] and lowerd3[n3][n1]: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lowerd3[n3][n1]) + lowerNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 + sinusNodeId = [] for n1 in range(elementsCountAround*2): if ((n1 - nMidCusp*2)%4) == 2: + sinusNodeId.append(None) continue node = nodes.createNode(nodeIdentifier, nodetemplateLinearS3) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, sinusx [n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, sinusd1[n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, sinusd2[n1]) + sinusNodeId.append(nodeIdentifier) nodeIdentifier += 1 + noduleNodeId = [ [], [] ] for n3 in range(2): for n1 in range(3): node = nodes.createNode(nodeIdentifier, nodetemplate) @@ -307,8 +313,10 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, noduled1[n3][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, noduled2[n3][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, noduled3[n3][n1]) + noduleNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 + upperNodeId = [ [], [] ] for n3 in range(2): for n1 in range(len(upperx[n3])): node = nodes.createNode(nodeIdentifier, nodetemplate if (upperd3[n3] and upperd3[n3][n1]) else nodetemplateLinearS3) @@ -318,8 +326,11 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, upperd2[n3][n1]) if upperd3[n3] and upperd3[n3][n1]: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, upperd3[n3][n1]) + upperNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 + + ################# # Create elements ################# @@ -332,14 +343,64 @@ def generateBaseMesh(cls, region, options, baseCentre=[ 0.0, 0.0, 0.0 ], axisSid bicubichermitelinear = eftfactory_bicubichermitelinear(mesh, useCrossDerivatives) eft = bicubichermitelinear.createEftNoCrossDerivatives() - tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + #tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) - elementIdentifier = startElementIdentifier = getMaximumElementIdentifier(mesh) + 1 + elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) elementtemplate1 = mesh.createElementtemplate() elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - scalefactors5hanging = [ -1.0, 0.5, 0.25, 0.125, 0.75 ] + # wall elements + for cusp in range(3): + n1 = cusp*2 - 1 + nMidCusp + n2 = n1*2 + for e in range(6): + eft1 = bicubichermitelinear.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + scalefactors = [ -1.0 ] + + if (e == 0) or (e == 5): + # 6 node collapsed wedge element expanding from zero width on outer wall of root, attaching to vertical part of cusp + if e == 0: + nids = [ lowerNodeId[0][n1], sinusNodeId[n2 + 1], upperNodeId[0][n2], upperNodeId[0][n2 + 1], lowerNodeId[1][n1], upperNodeId[1][n1] ] + remapEftNodeValueLabel(eft1, [ 1, 2 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, []) ]) + else: + nids = [ sinusNodeId[n2 + 3], lowerNodeId[0][n1 - 4], upperNodeId[0][n2 + 3], upperNodeId[0][n2 - 8], lowerNodeId[1][n1 - 4], upperNodeId[1][n1 - 4] ] + remapEftNodeValueLabel(eft1, [ 1, 2 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [1]) ]) + remapEftNodeValueLabel(eft1, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS1, []) + ln_map = [ 1, 2, 3, 4, 5, 5, 6, 6 ] + remapEftLocalNodes(eft1, 6, ln_map) + elif (e == 1) or (e == 4): + # 6 node collapsed wedge element on lower wall + if e == 1: + nids = [ lowerNodeId[0][n1], lowerNodeId[0][n1 + 1], sinusNodeId[n2 + 1], sinusNodeId[n2 + 2], lowerNodeId[1][n1], lowerNodeId[1][n1 + 1] ] + remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, []) ]) + else: + nids = [ lowerNodeId[0][n1 + 1], lowerNodeId[0][n1 - 4], sinusNodeId[n2 + 2], sinusNodeId[n2 + 3], lowerNodeId[1][n1 + 1], lowerNodeId[1][n1 - 4] ] + remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, []) ]) + remapEftNodeValueLabel(eft1, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) + ln_map = [ 1, 2, 3, 4, 5, 6, 5, 6 ] + remapEftLocalNodes(eft1, 6, ln_map) + else: + if e == 2: + nids = [ sinusNodeId[n2 + 1], sinusNodeId[n2 + 2], upperNodeId[0][n2 + 1], upperNodeId[0][n2 + 2], + lowerNodeId[1][n1], lowerNodeId[1][n1 + 1], upperNodeId[1][n1], upperNodeId[1][n1 + 1] ] + else: + nids = [ sinusNodeId[n2 + 2], sinusNodeId[n2 + 3], upperNodeId[0][n2 + 2], upperNodeId[0][n2 + 3], + lowerNodeId[1][n1 + 1], lowerNodeId[1][n1 - 4], upperNodeId[1][n1 + 1], upperNodeId[1][n1 - 4] ] + + result = elementtemplate1.defineField(coordinates, -1, eft1) + element = mesh.createElement(elementIdentifier, elementtemplate1) + result2 = element.setNodesByIdentifier(eft1, nids) + if scalefactors: + result3 = element.setScaleFactors(eft1, scalefactors) + else: + result3 = 7 + print('create arterial root wall', cusp, e, 'element',elementIdentifier, result, result2, result3, nids) + elementIdentifier += 1 + + for meshGroup in allMeshGroups: + meshGroup.addElement(element) fm.endChange() return annotationGroups @@ -353,13 +414,17 @@ def refineMesh(cls, meshrefinement, options): """ assert isinstance(meshrefinement, MeshRefinement) refineElementsCountSurface = options['Refine number of elements surface'] - MeshType_3d_heartventricles1.refineMesh(meshrefinement, options) - numberInXi1 = refineElementsCountSurface - numberInXi2 = refineElementsCountSurface - numberInXi3 = 1 - for i in range(0): # eventually: (24): # fixed number of elements - element = meshrefinement._sourceElementiterator.next() - meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3) + for cusp in range(3): + for e in range(6): + element = meshrefinement._sourceElementiterator.next() + numberInXi1 = refineElementsCountSurface + numberInXi2 = refineElementsCountSurface + numberInXi3 = 1 + if (e == 0) or (e == 5): + numberInXi1 = 1 + elif (e == 1) or (e == 4): + numberInXi2 = 1 + meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3) @classmethod def generateMesh(cls, region, options): From 8974654aaef167886cc0fd7d1e9e8238e71e0c31 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 27 Nov 2018 10:34:27 +1300 Subject: [PATCH 8/8] Define coordinates on all annotation points, transform --- scaffoldmaker/meshtypes/meshtype_3d_heart1.py | 10 ++-- .../meshtypes/meshtype_3d_heartventricles1.py | 47 ++++++++++--------- .../meshtype_3d_heartventriclesbase1.py | 19 +++++++- 3 files changed, 50 insertions(+), 26 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index 498330b1..2669e6b2 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/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 diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index 8f2e98af..cf65f91a 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/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 ]) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 862a368d..9724c29b 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -19,6 +19,7 @@ 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): ''' @@ -222,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() @@ -258,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