diff --git a/src/scaffoldmaker/annotation/heart_terms.py b/src/scaffoldmaker/annotation/heart_terms.py index f96194dd..3b7d202e 100644 --- a/src/scaffoldmaker/annotation/heart_terms.py +++ b/src/scaffoldmaker/annotation/heart_terms.py @@ -4,7 +4,7 @@ # convention: preferred name, preferred id, followed by any other ids and alternative names heart_terms = [ - ( "heart", "FMA:7088", "UBERON:0000948" ), + ( "heart", "UBERON:0000948", "FMA:7088" ), # ventricles ( "left ventricle myocardium", "FMA:9558" ), ( "right ventricle myocardium", "FMA:9535" ), @@ -52,7 +52,10 @@ ( "anterior cusp of pulmonary valve", "FMA:7249" ), ( "left cusp of pulmonary valve", "FMA:7247" ), # future: fiducial markers - ( "apex of heart", "FMA:7164", "UBERON:0002098") + ( "apex of heart", "UBERON:0002098", "FMA:7164"), + ( "crux of heart", "ILX:0777104", "FMA:7220" ), + ( "left atrium epicardium venous midpoint", "None"), # point at centre of pulmonary vein ostia on left atrium epicardium, on anterior/ventral side for rodents + ( "right atrium epicardium venous midpoint", "None") # point at centre of inferior & superior vena cavae on right atrium epicardium ] def get_heart_term(name : str): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py index 2d6dce10..4b2d92e8 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heart1.py @@ -54,15 +54,19 @@ def getOrderedOptionNames(): optionNames = MeshType_3d_heartventriclesbase1.getOrderedOptionNames() optionNamesAtria = MeshType_3d_heartatria1.getOrderedOptionNames() # insert new numbers of elements from atria; others came with ventriclesbase1 - optionNames.insert(7, 'Number of elements along vena cava inlet') - optionNames.insert(8, 'Number of elements over atria') + optionNames.insert(7, 'Number of elements along atrial appendages') + optionNames.insert(8, 'Number of elements along vena cava inlet') + optionNames.insert(9, 'Number of elements over atria') + optionNames.insert(10, 'Number of elements radial pulmonary vein annuli') # remove number of elements, unit scale and dependent options from atria options for key in [ + 'Number of elements along atrial appendages', 'Number of elements along vena cava inlet', 'Number of elements around atrial septum', 'Number of elements around left atrium free wall', 'Number of elements around right atrium free wall', 'Number of elements over atria', + 'Number of elements radial pulmonary vein annuli', 'Unit scale', 'Aorta outer plus diameter']: optionNamesAtria.remove(key) @@ -133,6 +137,8 @@ def generateBaseMesh(cls, region, options): ventriclesAnnotationGroups = MeshType_3d_heartventriclesbase1.generateBaseMesh(region, options) atriaAnnotationGroups = MeshType_3d_heartatria1.generateBaseMesh(region, options) annotationGroups = mergeAnnotationGroups(ventriclesAnnotationGroups, atriaAnnotationGroups) + heartGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("heart")) + cruxGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("crux of heart")) lFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("left fibrous ring")) rFibrousRingGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_heart_term("right fibrous ring")) @@ -211,6 +217,7 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + heartMeshGroup = heartGroup.getMeshGroup(mesh) lFibrousRingMeshGroup = lFibrousRingGroup.getMeshGroup(mesh) rFibrousRingMeshGroup = rFibrousRingGroup.getMeshGroup(mesh) @@ -233,7 +240,7 @@ def generateBaseMesh(cls, region, options): lavNodeId[0][0][n1], lavNodeId[0][0][n1 + 1], lavNodeId[0][1][n1], lavNodeId[0][1][n1 + 1], lavNodeId[1][0][n1], lavNodeId[1][0][n1 + 1], lavNodeId[1][1][n1], lavNodeId[1][1][n1 + 1]] scalefactors = None - meshGroups = [ lFibrousRingMeshGroup ] + meshGroups = [ heartMeshGroup, lFibrousRingMeshGroup ] if e == -1: # interatrial groove straddles left and right atria, collapsed to 6 node wedge @@ -295,7 +302,7 @@ def generateBaseMesh(cls, region, options): ravNodeId[0][0][n1], ravNodeId[0][0][n1 + 1], ravNodeId[0][1][n1], ravNodeId[0][1][n1 + 1], ravNodeId[1][0][n1], ravNodeId[1][0][n1 + 1], ravNodeId[1][1][n1], ravNodeId[1][1][n1 + 1]] scalefactors = None - meshGroups = [ rFibrousRingMeshGroup ] + meshGroups = [ heartMeshGroup, rFibrousRingMeshGroup ] if e == -1: # interatrial groove straddles left and right atria, collapsed to 6 node wedge @@ -350,7 +357,7 @@ def generateBaseMesh(cls, region, options): meshGroup.addElement(element) # fibrous ring septum: - meshGroups = [ lFibrousRingMeshGroup, rFibrousRingMeshGroup ] + meshGroups = [ heartMeshGroup, lFibrousRingMeshGroup, rFibrousRingMeshGroup ] for e in range(elementsCountAroundAtrialSeptum): eft1 = eftFibrousRing nlm = e - elementsCountAroundAtrialSeptum @@ -391,13 +398,13 @@ def generateBaseMesh(cls, region, options): # annotation fiducial points cruxElement = mesh.findElementByIdentifier(cruxElementId) cruxXi = [ 0.5, 0.5, 1.0 ] - cache.setMeshLocation(cruxElement, cruxXi) - result, cruxCoordinates = coordinates.evaluateReal(cache, 3) markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerName.assignString(cache, "crux of heart") + markerName.assignString(cache, cruxGroup.getName()) markerLocation.assignMeshLocation(cache, cruxElement, cruxXi) + for group in [ heartGroup, cruxGroup ]: + group.getNodesetGroup(nodes).addNode(markerPoint) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py index 865c3b59..ba5479b0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartatria1.py @@ -92,11 +92,11 @@ class MeshType_3d_heartatria1(Scaffold_base): 'Ostium diameter' : 0.16, 'Ostium length' : 0.08, 'Ostium wall thickness' : 0.04, - 'Ostium inter-vessel distance' : 0.13, + 'Ostium inter-vessel distance' : 0.12, 'Ostium inter-vessel height' : 0.01, 'Use linear through ostium wall' : False, 'Vessel end length factor' : 1.5, - 'Vessel inner diameter' : 0.1, + 'Vessel inner diameter' : 0.09, 'Vessel wall thickness' : 0.01, 'Vessel angle 1 degrees' : 0.0, 'Vessel angle 1 spread degrees' : 0.0, @@ -116,11 +116,11 @@ class MeshType_3d_heartatria1(Scaffold_base): 'Ostium diameter' : 0.16, 'Ostium length' : 0.1, 'Ostium wall thickness' : 0.04, - 'Ostium inter-vessel distance' : 0.13, + 'Ostium inter-vessel distance' : 0.12, 'Ostium inter-vessel height' : 0.01, 'Use linear through ostium wall' : False, 'Vessel end length factor' : 1.5, - 'Vessel inner diameter' : 0.1, + 'Vessel inner diameter' : 0.09, 'Vessel wall thickness' : 0.01, 'Vessel angle 1 degrees' : 0.0, 'Vessel angle 1 spread degrees' : 0.0, @@ -223,11 +223,13 @@ def getDefaultOptions(cls, parameterSetName='Default'): ventriclesbaseOptions['LV outlet inner diameter'] = 0.28 ventriclesbaseOptions['LV outlet wall thickness'] = 0.022 options = {} + options['Number of elements along atrial appendages'] = 2 options['Number of elements along vena cava inlet'] = 2 options['Number of elements around atrial septum'] = 3 options['Number of elements around left atrium free wall'] = 8 options['Number of elements around right atrium free wall'] = 6 - options['Number of elements over atria'] = 6 + options['Number of elements over atria'] = 8 + options['Number of elements radial pulmonary vein annuli'] = 2 options['Unit scale'] = 1.0 options['Atria base inner major axis length'] = 0.47 options['Atria base inner minor axis length'] = 0.41 @@ -252,7 +254,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Atrial base side incline degrees'] = 20.0 options['Atria venous anterior over'] = 0.7 options['Atria venous midpoint over'] = 0.41 - options['Left atrium venous midpoint left'] = 0.55 + options['Left atrium venous midpoint left'] = 0.5 options['Right atrium venous right'] = 0.4 options['Left atrial appendage angle axial degrees'] = 0.0 options['Left atrial appendage angle left degrees'] = 20.0 @@ -420,11 +422,13 @@ def getDefaultOptions(cls, parameterSetName='Default'): @staticmethod def getOrderedOptionNames(): return [ + 'Number of elements along atrial appendages', 'Number of elements along vena cava inlet', 'Number of elements around atrial septum', 'Number of elements around left atrium free wall', 'Number of elements around right atrium free wall', 'Number of elements over atria', + 'Number of elements radial pulmonary vein annuli', 'Unit scale', 'Atria base inner major axis length', 'Atria base inner minor axis length', @@ -545,8 +549,13 @@ def checkOptions(cls, options): :return: True if dependent options changed, otherwise False. ''' dependentChanges = False - if options['Number of elements along vena cava inlet'] < 1: - options['Number of elements along vena cava inlet'] = 1 + for key in [ + 'Number of elements along atrial appendages', + 'Number of elements along vena cava inlet', + 'Number of elements radial pulmonary vein annuli' + ]: + if options[key] < 1: + options[key] = 1 if options['Number of elements around atrial septum'] < 2: options['Number of elements around atrial septum'] = 2 for key in [ @@ -692,12 +701,15 @@ def generateBaseMesh(cls, region, options): :return: list of AnnotationGroup """ cls.updateSubScaffoldOptions(options) + + elementsCountAlongAtrialAppendages = options['Number of elements along atrial appendages'] elementsCountAlongVCInlet = options['Number of elements along vena cava inlet'] elementsCountAroundAtrialSeptum = options['Number of elements around atrial septum'] elementsCountAroundLeftAtriumFreeWall = options['Number of elements around left atrium free wall'] elementsCountAroundLeftAtrium = elementsCountAroundLeftAtriumFreeWall + elementsCountAroundAtrialSeptum elementsCountAroundRightAtriumFreeWall = options['Number of elements around right atrium free wall'] elementsCountOverAtria = options['Number of elements over atria'] + elementsCountRadialPVAnnuli = options['Number of elements radial pulmonary vein annuli'] unitScale = options['Unit scale'] aBaseInnerMajorMag = unitScale*0.5*options['Atria base inner major axis length'] @@ -780,13 +792,14 @@ def generateBaseMesh(cls, region, options): coordinates = findOrCreateFieldCoordinates(fm) cache = fm.createFieldcache() + heartGroup = AnnotationGroup(region, get_heart_term("heart")) laGroup = AnnotationGroup(region, get_heart_term("left atrium myocardium")) raGroup = AnnotationGroup(region, get_heart_term("right atrium myocardium")) aSeptumGroup = AnnotationGroup(region, get_heart_term("interatrial septum")) fossaGroup = AnnotationGroup(region, get_heart_term("fossa ovalis")) laaGroup = AnnotationGroup(region, get_heart_term("left auricle")) raaGroup = AnnotationGroup(region, get_heart_term("right auricle")) - annotationGroups = [ laGroup, raGroup, aSeptumGroup, fossaGroup, laaGroup, raaGroup ] + annotationGroups = [ heartGroup, laGroup, raGroup, aSeptumGroup, fossaGroup, laaGroup, raaGroup ] lpvOstiumSettings = lpvOstium.getScaffoldSettings() lpvCount = lpvOstiumSettings['Number of vessels'] @@ -818,7 +831,9 @@ def generateBaseMesh(cls, region, options): svcInletGroup = AnnotationGroup(region, get_heart_term("superior vena cava inlet")) ivcGroup = AnnotationGroup(region, get_heart_term("inferior vena cava")) svcGroup = AnnotationGroup(region, get_heart_term("superior vena cava")) - annotationGroups += [ ivcInletGroup, svcInletGroup, ivcGroup, svcGroup ] + laeVenousMidpointGroup = AnnotationGroup(region, get_heart_term("left atrium epicardium venous midpoint")) + raeVenousMidpointGroup = AnnotationGroup(region, get_heart_term("right atrium epicardium venous midpoint")) + annotationGroups += [ laeVenousMidpointGroup, ivcGroup, ivcInletGroup, raeVenousMidpointGroup, svcGroup, svcInletGroup ] # av boundary nodes are put in left and right fibrous ring groups only so they can be found by heart1 lFibrousRingGroup = AnnotationGroup(region, get_heart_term("left fibrous ring")) rFibrousRingGroup = AnnotationGroup(region, get_heart_term("right fibrous ring")) @@ -854,6 +869,7 @@ def generateBaseMesh(cls, region, options): elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) + heartMeshGroup = heartGroup.getMeshGroup(mesh) laMeshGroup = laGroup.getMeshGroup(mesh) raMeshGroup = raGroup.getMeshGroup(mesh) aSeptumMeshGroup = aSeptumGroup.getMeshGroup(mesh) @@ -922,7 +938,7 @@ def generateBaseMesh(cls, region, options): lpvOstiumDirection = [ (cosAngle*-td2[c] + sinAngle*td1[c]) for c in range(3) ] nodeIdentifier, elementIdentifier, (lpvox, lpvod1, lpvod2, lpvod3, lpvoNodeId, lpvoPositions) = \ generateOstiumMesh(region, lpvOstiumSettings, laTrackSurface, lpvOstiumPosition, lpvOstiumDirection, nodeIdentifier, elementIdentifier, - vesselMeshGroups=[ [ meshGroup ] for meshGroup in lpvMeshGroups ], ostiumMeshGroups=[ laMeshGroup ]) + vesselMeshGroups=[ [ heartMeshGroup, meshGroup ] for meshGroup in lpvMeshGroups ], ostiumMeshGroups=[ heartMeshGroup, laMeshGroup ]) if not commonLeftRightPvOstium: # create right pulmonary vein ostium @@ -937,7 +953,7 @@ def generateBaseMesh(cls, region, options): rpvOstiumDirection = [ (cosAngle*-td2[c] + sinAngle*td1[c]) for c in range(3) ] nodeIdentifier, elementIdentifier, (rpvox, rpvod1, rpvod2, rpvod3, rpvoNodeId, rpvoPositions) = \ generateOstiumMesh(region, rpvOstiumSettings, laTrackSurface, rpvOstiumPosition, rpvOstiumDirection, nodeIdentifier, elementIdentifier, - vesselMeshGroups=[ [ meshGroup ] for meshGroup in rpvMeshGroups ], ostiumMeshGroups=[ laMeshGroup ]) + vesselMeshGroups=[ [ heartMeshGroup, meshGroup ] for meshGroup in rpvMeshGroups ], ostiumMeshGroups=[ heartMeshGroup, laMeshGroup ]) # get points over interatrial septum on exterior groove agn1Mid = elementsCountOverRightAtriumNonVenousAnterior + elementsCountOverRightAtriumVenous//2 @@ -1283,15 +1299,36 @@ def generateBaseMesh(cls, region, options): ravtd3[0].append(agd3[0]) # get points on left atrium over appendage, from aorta to laa end on vestibule top - agn1 = 2 if commonLeftRightPvOstium else 1 - asn1 = elementsCountAroundAtrialSeptum + 1 if commonLeftRightPvOstium else elementsCountAroundAtrialSeptum - startScaling = 2.0 if commonLeftRightPvOstium else 1.0 + endDerivative = [ -d for d in lavtd2[1][lan1Mid] ] + # for not commonLeftRightPvOstium need a reliable location for laml - sample 3 across and use point 1 laoax, laoad1, laoad2, laoad3, laoaProportions = laTrackSurface.createHermiteCurvePoints( - 1.0 - ragProportions[agn1], 0.0, + lavtProportions[1][0], lavtProportions[1][1], lavtProportions[lan1Mid][0], lavtProportions[lan1Mid][1], - elementsCount = elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV, - derivativeStart = [ startScaling*d for d in agd1[agn1] ], - derivativeEnd = [ -1.0*d for d in lavtd2[1][lan1Mid] ]) + elementsCount=(elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV) if commonLeftRightPvOstium else 3, + derivativeStart=lavtd2[1][1], + derivativeEnd=endDerivative) + if (not commonLeftRightPvOstium) and (elementsCountOverSideLeftAtriumLPV != 2): + # resample from the laml line and merge + _laoax, _laoad1, _laoad2, _laoad3, _laoaProportions = laTrackSurface.createHermiteCurvePoints( + laoaProportions[1][0], laoaProportions[1][1], + lavtProportions[lan1Mid][0], lavtProportions[lan1Mid][1], + elementsCount=elementsCountOverSideLeftAtriumLPV, + derivativeStart=laoad1[1], + derivativeEnd=endDerivative) + laoax = laoax [:1] + _laoax + laoad1 = laoad1[:1] + _laoad1 + laoad2 = laoad2[:1] + _laoad2 + laoad3 = laoad3[:1] + _laoad3 + laoaProportions = laoaProportions[:1] + _laoaProportions + # d2 magnitudes are the same as d1 by default; want these to be even instead + # transition between start and end side derivatives d1 + mag0 = vector.magnitude(lavtd1[1][1]) + magn = vector.magnitude(lavtd1[1][lan1Mid]) + # could blend accurately with arc lengths; try indexes: + laoaCount = len(laoaProportions) + for n1 in range(1, laoaCount - 1): + xi = n1 / laoaCount + laoad2[n1] = vector.setMagnitude(laoad2[n1], (1.0 - xi)*mag0 + xi*magn) # get inner points laoax = [ [ None ], laoax ] laoad1 = [ [ None ], laoad1 ] @@ -1305,37 +1342,79 @@ def generateBaseMesh(cls, region, options): laoad3[0].append(d3) laoad3[1][n] = d3 # substitute known start and end coordinates - laoax [0][ 0] = asx [0][asn1] - laoad1[0][ 0] = asd1[0][asn1] - laoad2[0][ 0] = asd2[0][asn1] - laoad3[0][ 0] = asd3[0][asn1] - laoax [1][ 0] = agx [agn1] - laoad1[1][ 0] = agd1[agn1] - laoad2[1][ 0] = agd2[agn1] - laoad3[1][ 0] = agd3[agn1] for n3 in range(2): + laoax [n3][ 0] = lavtx [n3][1] + laoad1[n3][ 0] = lavtd1[n3][1] + laoad2[n3][ 0] = lavtd2[n3][1] + laoad3[n3][ 0] = lavtd3[n3][1] laoax [n3][-1] = lavtx [n3][lan1Mid] laoad1[n3][-1] = lavtd1[n3][lan1Mid] laoad2[n3][-1] = lavtd2[n3][lan1Mid] laoad3[n3][-1] = lavtd3[n3][lan1Mid] - # smooth d2 to fit adjacent LPV derivative 2 - if commonLeftRightPvOstium: - n1Start = 1 - n1lpvStart = interp.getNearestPointIndex(lpvox[1], agx[agn1Mid]) - elementsCountOverLeftAtriumVenous//2 - 1 - else: - n1Start = elementsCountAroundLeftAtriumRPV + 1 - n1lpvStart = -2 - elementsCountOverLeftAtriumVenous//2 - for n1 in range(n1Start, len(laoax[1]) - 1): - n1lpv = n1lpvStart - (n1 - n1Start) - #print('n1',n1,'n1lpv',n1lpv) - #print('smooth laoa n1',n1,'lpv',n1lpv) + + if not commonLeftRightPvOstium: + # get left atrium venous mid line from 3rd point on laoa to point on lavt at venous midpoint + # find and pass through midpoint between left and right PVs + n1lpv = 0 + n1rpv = elementsCountOverLeftAtriumVenous + elementsCountAroundLeftAtriumRPV + n1End = elementsCountAroundLeftAtriumFreeWall - elementsCountAroundLeftAtriumRPV + rpvoProportion1, rpvoProportion2 = laTrackSurface.getProportion(rpvoPositions[n1rpv]) + lpvoProportion1, lpvoProportion2 = laTrackSurface.getProportion(lpvoPositions[n1lpv]) + mpx, mpd1, mpd2, _, mpProportions = laTrackSurface.createHermiteCurvePoints(rpvoProportion1, rpvoProportion2, lpvoProportion1, lpvoProportion2, + elementsCount = 2, derivativeStart = rpvod2[1][n1rpv], derivativeEnd = [ -d for d in lpvod2[1][n1lpv] ]) + # scale mid derivative 2 to be mean of d1 in LPV, RPV + d2mag = 0.5*vector.magnitude(lpvod1[1][n1lpv]) + 0.5*vector.magnitude(rpvod1[1][n1rpv]) + mpd2[1] = vector.setMagnitude(mpd2[1], d2mag) + lamlx, lamld2, lamld1, lamld3, lamlProportions = laTrackSurface.createHermiteCurvePoints( + laoaProportions[1][0], laoaProportions[1][1], + mpProportions[1][0], mpProportions[1][1], + elementsCount = elementsCountOverLeftAtriumVenous//2 + 1, + derivativeStart = laoad2[1][1], + derivativeEnd = mpd2[1]) + _lamlx, _lamld2, _lamld1, _lamld3, _lamlProportions = laTrackSurface.createHermiteCurvePoints( + mpProportions[1][0], mpProportions[1][1], + lavtProportions[n1End][0], lavtProportions[n1End][1], + elementsCount = elementsCountOverLeftAtriumVenous//2 + 1, + derivativeStart = mpd2[1], + derivativeEnd = [ -d for d in lavtd2[0][n1End]]) + lamlx += _lamlx [1:] + lamld1 += _lamld1[1:] + lamld2 += _lamld2 [1:] + lamld3 += _lamld3 [1:] + lamlProportions += _lamlProportions[1:] + # reverse d1 + lamld1 = [ [ -d for d in d1 ] for d1 in lamld1 ] + # smooth d2 + lamld2 = interp.smoothCubicHermiteDerivativesLine(lamlx, lamld2, fixAllDirections=True, fixStartDerivative=True, fixEndDerivative=True) + # give d1 the same magnitude as the smoothed d2 + for n in range(1, len(lamld1)): + lamld1[n] = vector.setMagnitude(lamld1[n], vector.magnitude(lamld2[n])) + # get inner points + lamlx = [ [ None ], lamlx ] + lamld1 = [ [ None ], lamld1 ] + lamld2 = [ [ None ], lamld2 ] + lamld3 = [ [ None ], lamld3 ] + lamlProportions += _lamlProportions[1:] + for n2 in range(1, elementsCountOverLeftAtriumVenous + 3): + x, d2, d1, d3 = interp.projectHermiteCurvesThroughWall(lamlx[1], lamld2[1], [ [ -d for d in d1 ] for d1 in lamld1[1] ], n2, -laVenousFreeWallThickness) + lamlx [0].append(x) + lamld1[0].append([ -d for d in d1 ]) + lamld2[0].append(d2) + lamld3[0].append(d3) + lamld3[1][n2] = d3 + # substitute known start and end coordinates for n3 in range(2): - nx = [ laoax [n3][n1], lpvox [n3][n1lpv] ] - nd1 = [ laoad2[n3][n1], [ -d for d in lpvod2[n3][n1lpv] ] ] - laoad2[n3][n1] = interp.smoothCubicHermiteDerivativesLine(nx, nd1, fixAllDirections = True, fixEndDerivative = True)[0] + lamlx [n3][ 0] = laoax [n3][1] + lamld1[n3][ 0] = laoad1[n3][1] + lamld2[n3][ 0] = laoad2[n3][1] + lamld3[n3][ 0] = laoad3[n3][1] + lamlx [n3][-1] = lavtx [n3][n1End] + lamld1[n3][-1] = lavtd1[n3][n1End] + lamld2[n3][-1] = lavtd2[n3][n1End] + lamld3[n3][-1] = lavtd3[n3][n1End] if not commonLeftRightPvOstium: - # get points on row above left atrium venous anterior, from interatrial septum to nearest point on LPV ostium + # get points on row above left atrium venous anterior, from interatrial septum to laml[1] # sample points up to venous midpoint, between RPV and laoa agn1va = elementsCountOverLeftAtriumNonVenousAnterior lavbx = [ agx [agn1va] ] @@ -1344,134 +1423,56 @@ def generateBaseMesh(cls, region, options): lavbd3 = [ agd3[agn1va] ] lavbd2inner = [ None ] lavbProportions = [ [ aVenousAnteriorOver, 0.0 ] ] - for n1 in range(1, elementsCountAroundLeftAtriumRPV + 1): - n1rpv = -elementsCountOverLeftAtriumVenous//2 - n1 + for n1 in range(elementsCountAroundLeftAtriumRPV - 1): + n1rpv = -elementsCountOverLeftAtriumVenous//2 - n1 - 1 rpvoProportion1, rpvoProportion2 = laTrackSurface.getProportion(rpvoPositions[n1rpv]) + startDerivative = [ (1.5*laoad2[1][n1][c] - 0.5*laoad1[1][n1][c]) for c in range(3) ] if (n1 == 0) else laoad2[1][n1] # GRC fudge factor _x, _d2, _d1, _d3, _proportions = laTrackSurface.createHermiteCurvePoints( laoaProportions[n1][0], laoaProportions[n1][1], rpvoProportion1, rpvoProportion2, - elementsCount = 2, - derivativeStart = laoad2[1][n1], # GRC automatic value of equal magnitude to d1, looks ok - derivativeEnd = [ -d for d in rpvod2[1][n1rpv] ]) + elementsCount=1 + elementsCountRadialPVAnnuli, + derivativeStart=startDerivative, + derivativeEnd=[ -d for d in rpvod2[1][n1rpv] ]) lavbx .append(_x [1]) lavbd1.append([ -d for d in _d1[1] ]) lavbd2.append(_d2[1]) lavbd3.append(_d3[1]) lavbProportions.append(_proportions[1]) # get precise inner d2 - _, _d2inner, _, _ = interp.projectHermiteCurvesThroughWall(_x, _d2, _d1, 1, -laVenousFreeWallThickness) + _d2inner = interp.projectHermiteCurvesThroughWall(_x, _d2, _d1, 1, -laVenousFreeWallThickness)[1] lavbd2inner.append(_d2inner) - # add nearest point on LPV ostium - n1lpv = -elementsCountOverLeftAtriumVenous//2 - lpvoProportion1, lpvoProportion2 = laTrackSurface.getProportion(lpvoPositions[n1lpv]) - lavbx .append(lpvox [1][n1lpv]) - lavbd1.append([ -d for d in lpvod2[1][n1lpv] ]) - lavbd2.append(lpvod1[1][n1lpv]) - lavbd3.append(lpvod3[1][n1lpv]) - lavbProportions.append([ lpvoProportion1, lpvoProportion2 ]) + # use end values from laml[1][1], with d1 at 30 degree angle + cos30 = math.cos(math.pi/6.0) + sin30 = math.sin(math.pi/6.0) + lavbx .append(lamlx [1][1]) + lavbd1.append([ (cos30*lamld1[1][1][c] + sin30*lamld2[1][1][c]) for c in range(3) ]) + lavbd2.append(lamld2[1][1]) + lavbd3.append(lamld3[1][1]) + lavbd2inner.append(lamld2[0][1]) + lavbProportions.append(lamlProportions[1]) # smooth d1: lavbd1 = interp.smoothCubicHermiteDerivativesLine(lavbx, lavbd1, fixAllDirections = True, - fixStartDerivative = True, fixEndDerivative = True, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + fixStartDerivative=True, fixEndDerivative=True, magnitudeScalingMode=interp.DerivativeScalingMode.ARITHMETIC_MEAN) # get inner points - lavbx = [ [ None ], lavbx ] - lavbd1 = [ [ None ], lavbd1 ] - lavbd2 = [ [ None ], lavbd2 ] - lavbd3 = [ [ None ], lavbd3 ] - for n1 in range(1, len(lavbx[1]) - 1): + asn1va = elementsCountAroundAtrialSeptum - 1 + elementsCountOverLeftAtriumNonVenousAnterior + lavbx = [ [ asx [0][asn1va] ], lavbx ] + lavbd1 = [ [ asd1[0][asn1va] ], lavbd1 ] + lavbd2 = [ [ asd2[0][asn1va] ], lavbd2 ] + lavbd3 = [ [ asd3[0][asn1va] ], lavbd3 ] + for n1 in range(1, len(lavbx[1])): x, d1, _, d3 = interp.projectHermiteCurvesThroughWall(lavbx[1], lavbd1[1], lavbd2[1], n1, -laVenousFreeWallThickness) lavbx [0].append(x) lavbd1[0].append(d1) lavbd2[0].append(lavbd2inner[n1]) lavbd3[0].append(d3) lavbd3[1][n1] = d3 - lavbx [0].append(None) - lavbd1[0].append(None) - lavbd2[0].append(None) - lavbd3[0].append(None) - # transfer/substitute known start and end coordinates - asn1va = elementsCountAroundAtrialSeptum - 1 + elementsCountOverLeftAtriumNonVenousAnterior - lavbx [0][0] = asx [0][asn1va] - lavbd1[0][0] = asd1[0][asn1va] - lavbd2[0][0] = asd2[0][asn1va] - lavbd3[0][0] = asd3[0][asn1va] - lavbx [1][0] = agx [0][agn1va] - lavbd1[1][0] = agd1[0][agn1va] - lavbd2[1][0] = agd2[0][agn1va] - lavbd3[1][0] = agd3[0][agn1va] - for n3 in range(2): - lavbx [1][-1] = lpvox [n3][n1lpv] - lavbd1[1][-1] = lpvod1[n3][n1lpv] - lavbd2[1][-1] = lpvod2[n3][n1lpv] - lavbd3[1][-1] = lpvod3[n3][n1lpv] + # copy end d1 values back to lamld1 + lamld1[0][1] = lavbd1[0][-1] + lamld1[1][1] = lavbd1[1][-1] agn1vp = elementsCountOverLeftAtriumNonVenousAnterior + elementsCountOverLeftAtriumVenous - if commonLeftRightPvOstium: - # get points on row above left atrium venous posterior, from interatrial septum to laoa[-2] - # sample points up to venous midpoint, between LPV and vestibule top - lavqx = [ agx [agn1vp] ] - lavqd1 = [ agd1[agn1vp] ] - lavqd2 = [ agd2[agn1vp] ] - lavqd3 = [ agd3[agn1vp] ] - lavqd2inner = [ None ] - lavqProportions = [ [ laVenousLimitPosterior, 0.0 ] ] - n1lpvStart = interp.getNearestPointIndex(lpvox[1], agx[agn1Mid]) + elementsCountOverLeftAtriumVenous//2 - elementsCountAroundLpvOstium - for n1 in range(1, elementsCountAroundLeftAtriumRPV + elementsCountAroundLeftAtriumLPV): - n1lpv = n1lpvStart + n1 - n1cs = elementsCountAroundLeftAtriumFreeWall - n1 - lpvoProportion1, lpvoProportion2 = laTrackSurface.getProportion(lpvoPositions[n1lpv]) - _x, _d2, _d1, _d3, _proportions = laTrackSurface.createHermiteCurvePoints( - lpvoProportion1, lpvoProportion2, - lavtProportions[n1cs][0], lavtProportions[n1cs][1], - elementsCount = 2, - derivativeStart = lpvod2[1][n1lpv], derivativeEnd = [ -d for d in lavtd2[1][n1cs] ]) - lavqx .append(_x [1]) - lavqd1.append([ -d for d in _d1[1] ]) - lavqd2.append(_d2[1]) - lavqd3.append(_d3[1]) - lavqProportions.append(_proportions[1]) - # get precise inner d2 - _, _d2inner, _, _ = interp.projectHermiteCurvesThroughWall(_x, _d2, _d1, 1, -laVenousFreeWallThickness) - lavqd2inner.append(_d2inner) - lavqd2inner.append(laoad1[0][-2]) - # add 2nd last point on laoa - lavqx .append(laoax [1][-2]) - lavqd1.append([ -d for d in laoad2[1][-2] ]) - lavqd2.append(laoad1[1][-2]) - lavqd3.append(laoad3[1][-2]) - lavqProportions.append(laoaProportions[-2]) - # smooth d1: - lavqd1 = interp.smoothCubicHermiteDerivativesLine(lavqx, lavqd1, fixAllDirections = True, - fixStartDerivative = True, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) - # get inner points - lavqx = [ [ None ], lavqx ] - lavqd1 = [ [ None ], lavqd1 ] - lavqd2 = [ [ None ], lavqd2 ] - lavqd3 = [ [ None ], lavqd3 ] - for n1 in range(1, len(lavqx[1])): - x, d1, _, d3 = interp.projectHermiteCurvesThroughWall(lavqx[1], lavqd1[1], lavqd2[1], n1, -laVenousFreeWallThickness) - lavqx [0].append(x) - lavqd1[0].append(d1) - lavqd2[0].append(lavqd2inner[n1]) - lavqd3[0].append(d3) - lavqd3[1][n1] = d3 - # transfer/substitute known start and end coordinates - asn1vp = -1 - lavqx [0][0] = asx [0][asn1vp] - lavqd1[0][0] = asd1[0][asn1vp] - lavqd2[0][0] = asd2[0][asn1vp] - lavqd3[0][0] = asd3[0][asn1vp] - lavqx [1][0] = asx [0][agn1vp] - lavqd1[1][0] = asd1[0][agn1vp] - lavqd2[1][0] = asd2[0][agn1vp] - lavqd3[1][0] = asd3[0][agn1vp] - for n3 in range(2): - laoad2[n3][-2] = [ -d for d in lavqd1[n3][-1] ] # use final d1 on laoa - lavqx [n3][-1] = laoax [n3][-2] - lavqd1[n3][-1] = laoad1[n3][-2] - lavqd2[n3][-1] = laoad2[n3][-2] - lavqd3[n3][-1] = laoad3[n3][-2] - else: - # get points on row above left atrium venous posterior, from interatrial septum to nearest point on LPV ostium + if not commonLeftRightPvOstium: + # get points on row above left atrium venous posterior, from interatrial septum to laml[-2] # sample points up to venous midpoint, between RPV and vestibule top lavqx = [ agx [agn1vp] ] lavqd1 = [ agd1[agn1vp] ] @@ -1479,134 +1480,51 @@ def generateBaseMesh(cls, region, options): lavqd3 = [ agd3[agn1vp] ] lavqd2inner = [ None ] lavqProportions = [ [ laVenousLimitPosterior, 0.0 ] ] - for n1 in range(1, elementsCountAroundLeftAtriumRPV + 1): + for n1 in range(1, elementsCountAroundLeftAtriumRPV): n1rpv = elementsCountOverLeftAtriumVenous//2 + n1 n1cs = elementsCountAroundLeftAtriumFreeWall - n1 rpvoProportion1, rpvoProportion2 = laTrackSurface.getProportion(rpvoPositions[n1rpv]) _x, _d2, _d1, _d3, _proportions = laTrackSurface.createHermiteCurvePoints( rpvoProportion1, rpvoProportion2, lavtProportions[n1cs][0], lavtProportions[n1cs][1], - elementsCount = 2, + elementsCount = 1 + elementsCountRadialPVAnnuli, derivativeStart = rpvod2[1][n1rpv], derivativeEnd = [ -d for d in lavtd2[1][n1cs] ]) - lavqx .append(_x [1]) - lavqd1.append([ -d for d in _d1[1] ]) - lavqd2.append(_d2[1]) - lavqd3.append(_d3[1]) - lavqProportions.append(_proportions[1]) + lavqx .append(_x [-2]) + lavqd1.append([ -d for d in _d1[-2] ]) + lavqd2.append(_d2[-2]) + lavqd3.append(_d3[-2]) + lavqProportions.append(_proportions[-2]) # get precise inner d2 - _, _d2inner, _, _ = interp.projectHermiteCurvesThroughWall(_x, _d2, _d1, 1, -laVenousFreeWallThickness) + _d2inner = interp.projectHermiteCurvesThroughWall(_x, _d2, _d1, len(_x) - 2, -laVenousFreeWallThickness)[1] lavqd2inner.append(_d2inner) - # add nearest point on LPV ostium - n1lpv = elementsCountOverLeftAtriumVenous//2 - lpvoProportion1, lpvoProportion2 = laTrackSurface.getProportion(lpvoPositions[n1lpv]) - lavqx .append(lpvox [1][n1lpv]) - lavqd1.append([ -d for d in lpvod2[1][n1lpv] ]) - lavqd2.append(lpvod1[1][n1lpv]) - lavqd3.append(lpvod3[1][n1lpv]) - lavqProportions.append([ lpvoProportion1, lpvoProportion2 ]) + # use end values from laml[1][-2], with d1 at 30 degree angle + cos30 = math.cos(math.pi/6.0) + sin30 = math.sin(math.pi/6.0) + lavqx .append(lamlx [1][-2]) + lavqd1.append([ (cos30*lamld1[1][-2][c] - sin30*lamld2[1][-2][c]) for c in range(3) ]) + lavqd2.append(lamld2[1][-2]) + lavqd3.append(lamld3[1][-2]) + lavqd2inner.append(lamld2[0][-1]) + lavqProportions.append(lamlProportions[-2]) # smooth d1: lavqd1 = interp.smoothCubicHermiteDerivativesLine(lavqx, lavqd1, fixAllDirections = True, - fixStartDerivative = True, fixEndDerivative = True, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + fixStartDerivative=True, fixEndDerivative=True, magnitudeScalingMode=interp.DerivativeScalingMode.ARITHMETIC_MEAN) # get inner points - lavqx = [ [ None ], lavqx ] - lavqd1 = [ [ None ], lavqd1 ] - lavqd2 = [ [ None ], lavqd2 ] - lavqd3 = [ [ None ], lavqd3 ] - for n1 in range(1, len(lavqx[1]) - 1): + asn1vp = -1 + lavqx = [ [ asx [0][asn1vp] ], lavqx ] + lavqd1 = [ [ asd1[0][asn1vp] ], lavqd1 ] + lavqd2 = [ [ asd2[0][asn1vp] ], lavqd2 ] + lavqd3 = [ [ asd3[0][asn1vp] ], lavqd3 ] + for n1 in range(1, len(lavqx[1])): x, d1, _, d3 = interp.projectHermiteCurvesThroughWall(lavqx[1], lavqd1[1], lavqd2[1], n1, -laVenousFreeWallThickness) lavqx [0].append(x) lavqd1[0].append(d1) lavqd2[0].append(lavqd2inner[n1]) lavqd3[0].append(d3) lavqd3[1][n1] = d3 - lavqx [0].append(None) - lavqd1[0].append(None) - lavqd2[0].append(None) - lavqd3[0].append(None) - # transfer/substitute known start and end coordinates - asn1vp = -1 - lavqx [0][0] = asx [0][asn1vp] - lavqd1[0][0] = asd1[0][asn1vp] - lavqd2[0][0] = asd2[0][asn1vp] - lavqd3[0][0] = asd3[0][asn1vp] - lavqx [1][0] = asx [0][agn1vp] - lavqd1[1][0] = asd1[0][agn1vp] - lavqd2[1][0] = asd2[0][agn1vp] - lavqd3[1][0] = asd3[0][agn1vp] - for n3 in range(2): - lavqx [n3][-1] = lpvox [n3][n1lpv] - lavqd1[n3][-1] = lpvod1[n3][n1lpv] - lavqd2[n3][-1] = lpvod2[n3][n1lpv] - lavqd3[n3][-1] = lpvod3[n3][n1lpv] - - if not commonLeftRightPvOstium: - # get left atrium venous mid line from 2nd last points on lavb to lavq - # find and pass through midpoint between left and right PVs - n1lpv = 0 - n1rpv = elementsCountOverLeftAtriumVenous + elementsCountAroundLeftAtriumRPV - rpvoProportion1, rpvoProportion2 = laTrackSurface.getProportion(rpvoPositions[n1rpv]) - lpvoProportion1, lpvoProportion2 = laTrackSurface.getProportion(lpvoPositions[n1lpv]) - mpx, mpd1, mpd2, _, mpProportions = laTrackSurface.createHermiteCurvePoints(rpvoProportion1, rpvoProportion2, lpvoProportion1, lpvoProportion2, - elementsCount = 2, derivativeStart = rpvod2[1][n1rpv], derivativeEnd = [ -d for d in lpvod2[1][n1lpv] ]) - # scale mid derivative 2 to be mean of d1 in LPV, RPV - d2mag = 0.5*vector.magnitude(lpvod1[1][n1lpv]) + 0.5*vector.magnitude(rpvod1[1][n1rpv]) - mpd2[1] = vector.setMagnitude(mpd2[1], d2mag) - lamlx, lamld2, lamld1, lamld3, lamlProportions = laTrackSurface.createHermiteCurvePoints( - lavbProportions[elementsCountAroundLeftAtriumRPV][0], lavbProportions[elementsCountAroundLeftAtriumRPV][1], - mpProportions[1][0], mpProportions[1][1], - elementsCount = elementsCountOverLeftAtriumVenous//2, - derivativeStart = [ (lavbd1[1][elementsCountAroundLeftAtriumRPV][c] + lavbd2[1][elementsCountAroundLeftAtriumRPV][c]) for c in range(3) ], - derivativeEnd = mpd2[1]) - _lamlx, _lamld2, _lamld1, _lamld3, _lamlProportions = laTrackSurface.createHermiteCurvePoints( - mpProportions[1][0], mpProportions[1][1], - lavqProportions[elementsCountAroundLeftAtriumRPV][0], lavqProportions[elementsCountAroundLeftAtriumRPV][1], - elementsCount = elementsCountOverLeftAtriumVenous//2, - derivativeStart = mpd2[1], - derivativeEnd = [ (-lavqd1[1][elementsCountAroundLeftAtriumRPV][c] + lavqd2[1][elementsCountAroundLeftAtriumRPV][c]) for c in range(3) ]) - lamlx += _lamlx [1:] - lamld1 += _lamld1[1:] - lamld2 += _lamld2 [1:] - lamld3 += _lamld3 [1:] - lamlProportions += _lamlProportions[1:] - # reverse d1 - lamld1 = [ [ -d for d in d1 ] for d1 in lamld1 ] - if elementsCountOverLeftAtriumVenous == 2: - # recalculate central d2 to match end derivatives - lamld2[1] = interp.smoothCubicHermiteDerivativesLine(lamlx, lamld2, fixAllDirections = True, - fixStartDerivative = True, fixEndDerivative = True, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN)[1] - # get inner points - lamlx = [ [ None ], lamlx ] - lamld1 = [ [ None ], lamld1 ] - lamld2 = [ [ None ], lamld2 ] - lamld3 = [ [ None ], lamld3 ] - lamlProportions += _lamlProportions[1:] - for n2 in range(1, elementsCountOverLeftAtriumVenous + 1): - x, d2, d1, d3 = interp.projectHermiteCurvesThroughWall(lamlx[1], lamld2[1], [ [ -d for d in d1 ] for d1 in lamld1[1] ], n2, -laVenousFreeWallThickness) - lamlx [0].append(x) - lamld1[0].append([ -d for d in d1 ]) - lamld2[0].append(d2) - lamld3[0].append(d3) - lamld3[1][n2] = d3 - # smooth d1 between RPV, LPV - for n2 in range(1, elementsCountOverLeftAtriumVenous): - n1lpv = n2 - elementsCountOverLeftAtriumVenous//2 - 1 - n1rpv = n2 - elementsCountOverLeftAtriumVenous - elementsCountAroundLeftAtriumRPV - 1 - for n3 in range(2): - nx = [ rpvox [n3][n1rpv], lamlx [n3][n2], lpvox [n3][n1lpv] ] - nd1 = [ rpvod1[n3][n1rpv], lamld1[n3][n2], [ -d for d in lpvod1[n3][n1lpv] ] ] - d1 = interp.smoothCubicHermiteDerivativesLine(nx, nd1, fixAllDirections = True, - fixStartDerivative = True, fixEndDerivative = True, magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN)[1] - lamld1[n3][n2] = d1 - # substitute known start and end coordinates - for n3 in range(2): - lamlx [n3][ 0] = lavbx [n3][elementsCountAroundLeftAtriumRPV] - lamld1[n3][ 0] = lavbd1[n3][elementsCountAroundLeftAtriumRPV] - lamld2[n3][ 0] = lavbd2[n3][elementsCountAroundLeftAtriumRPV] - lamld3[n3][ 0] = lavbd3[n3][elementsCountAroundLeftAtriumRPV] - lamlx [n3][-1] = lavqx [n3][elementsCountAroundLeftAtriumRPV] - lamld1[n3][-1] = lavqd1[n3][elementsCountAroundLeftAtriumRPV] - lamld2[n3][-1] = lavqd2[n3][elementsCountAroundLeftAtriumRPV] - lamld3[n3][-1] = lavqd3[n3][elementsCountAroundLeftAtriumRPV] + # copy end d1 values back to lamld1 + lamld1[0][-2] = lavqd1[0][-1] + lamld1[1][-2] = lavqd1[1][-1] # get points on right atrium along crista terminalis from aorta to posterior venous limit xi = (1.0 - aVenousMidpointOver - ravtProportions[ran1Aorta][0])/(ravtProportions[ran1Ctp][0] - ravtProportions[ran1Aorta][0]) @@ -1722,8 +1640,8 @@ def generateBaseMesh(cls, region, options): ractProportions[nc][0], ractProportions[nc][1], ravtProportions[ran1raap][0], ravtProportions[ran1raap][1], elementsCount = elementsCountOverSideRightAtriumPouch, - #derivativeStart = [ (ractd2[1][nc][c] - ractd1[1][nc][c]) for c in range(3) ], - derivativeStart = [ -d for d in ractd1[1][nc] ], + derivativeStart = [ (0.25*ractd2[1][nc][c] - 1.5*ractd1[1][nc][c]) for c in range(3) ], + #derivativeStart = [ -d for d in ractd1[1][nc] ], derivativeEnd = [ -d for d in ravtd2[1][ran1raap] ]) # get inner points raapx = [ [ None ], raapx ] @@ -1755,7 +1673,7 @@ def generateBaseMesh(cls, region, options): ravtProportions[ran1raaq][0], ravtProportions[ran1raaq][1], elementsCount = elementsCountOverAtria//2 + elementsCountOverCristaTerminalisAnterior - 3, #derivativeStart = [ (ractd2[1][nc][c] - 0.5*ractd1[1][nc][c]) for c in range(3) ], - derivativeStart = [ (ractd2[1][nc][c] -ractd1[1][nc][c]) for c in range(3) ], + derivativeStart = [ (0.5*ractd2[1][nc][c] - ractd1[1][nc][c]) for c in range(3) ], derivativeEnd = [ -d for d in ravtd2[1][ran1raaq] ]) # get inner points raaqx = [ [ None ], raaqx ] @@ -1910,7 +1828,7 @@ def generateBaseMesh(cls, region, options): lavtNodeId = [ [ asNodeId[0][elementsCountAroundAtrialSeptum] ], [ agNodeId[1] ] ] for n3 in range(2): for n1 in range(1, len(lavtx[n3]) - 1): - if (elementsCountAroundLeftAtriumAorta < n1 < lan1Mid) or ((not commonLeftRightPvOstium) and (n1 == elementsCountAroundLeftAtriumAorta)): + if elementsCountAroundLeftAtriumAorta < n1 < lan1Mid: # left atrial appendage lavtNodeId[n3].append(None) continue @@ -1945,9 +1863,7 @@ def generateBaseMesh(cls, region, options): ravtNodeId[1].append(agNodeId[1]) # create nodes on left atrium over appendage - agn1 = 2 if commonLeftRightPvOstium else 1 - asn1 = (elementsCountAroundAtrialSeptum + 1) if commonLeftRightPvOstium else elementsCountAroundAtrialSeptum - laoaNodeId = [ [ asNodeId[0][asn1] ], [ agNodeId[agn1] ] ] + laoaNodeId = [ [ lavtNodeId[0][1] ], [ lavtNodeId[1][1] ] ] for n3 in range(2): for n1 in range(1, len(laoax[n3]) - 1): node = nodes.createNode(nodeIdentifier, nodetemplate) @@ -1959,12 +1875,25 @@ def generateBaseMesh(cls, region, options): laoaNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 laoaNodeId[n3].append(lavtNodeId[n3][lan1Mid]) - if not commonLeftRightPvOstium: - # use second laoa node as aorta vestibule top node - lavtNodeId[n3][lan1Aorta] = laoaNodeId[n3][1] if not commonLeftRightPvOstium: - # create nodes on row above left atrium venous anterior to LPV ostium + # create left atrium venous midline nodes + lamlNodeId = [ [], [] ] + for n3 in range(2): + lamlNodeId[n3].append(laoaNodeId[n3][2]) + for n2 in range(1, len(lamlx[n3]) - 1): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lamlx [n3][n2]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lamld1[n3][n2]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lamld2[n3][n2]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lamld3[n3][n2]) + lamlNodeId[n3].append(nodeIdentifier) + nodeIdentifier += 1 + lamlNodeId[n3].append(lavtNodeId[n3][-elementsCountAroundLeftAtriumRPV - 1]) + + if not commonLeftRightPvOstium: + # create nodes on row above left atrium venous anterior to laml[1] lavbNodeId = [ [ asNodeId[0][asn1va] ], [ agNodeId[elementsCountOverLeftAtriumNonVenousAnterior] ] ] n1lpv = -elementsCountOverLeftAtriumVenous//2 for n3 in range(2): @@ -1977,41 +1906,23 @@ def generateBaseMesh(cls, region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lavbd3[n3][n1]) lavbNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 - lavbNodeId[n3].append(lpvoNodeId[n3][n1lpv]) - - # create nodes on row above left atrium venous posterior to LPV ostium - lavqNodeId = [ [ asNodeId[0][-1] ], [ agNodeId[elementsCountOverLeftAtriumNonVenousAnterior + elementsCountOverLeftAtriumVenous ] ] ] - n1lpv = elementsCountOverLeftAtriumVenous//2 - for n3 in range(2): - for n1 in range(1, len(lavqx[n3]) - 1): - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lavqx [n3][n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lavqd1[n3][n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lavqd2[n3][n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lavqd3[n3][n1]) - lavqNodeId[n3].append(nodeIdentifier) - nodeIdentifier += 1 - if commonLeftRightPvOstium: - lavqNodeId[n3].append(laoaNodeId[n3][-2]) - else: - lavqNodeId[n3].append(lpvoNodeId[n3][n1lpv]) + lavbNodeId[n3].append(lamlNodeId[n3][1]) + # create nodes on row above left atrium venous posterior to laml[-2] if not commonLeftRightPvOstium: - # create left atrium venous midline nodes - lamlNodeId = [ [], [] ] + lavqNodeId = [ [ asNodeId[0][-1] ], [ agNodeId[elementsCountOverLeftAtriumNonVenousAnterior + elementsCountOverLeftAtriumVenous ] ] ] + n1lpv = elementsCountOverLeftAtriumVenous//2 for n3 in range(2): - lamlNodeId[n3].append(lavbNodeId[n3][elementsCountAroundLeftAtriumRPV]) - for n2 in range(1, len(lamlx[n3]) - 1): + for n1 in range(1, len(lavqx[n3]) - 1): node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lamlx [n3][n2]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lamld1[n3][n2]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lamld2[n3][n2]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lamld3[n3][n2]) - lamlNodeId[n3].append(nodeIdentifier) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, lavqx [n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, lavqd1[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, lavqd2[n3][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, lavqd3[n3][n1]) + lavqNodeId[n3].append(nodeIdentifier) nodeIdentifier += 1 - lamlNodeId[n3].append(lavqNodeId[n3][elementsCountAroundLeftAtriumRPV]) + lavqNodeId[n3].append(lamlNodeId[n3][-2]) # create right atrium crista terminalis nodes ractNodeId = [ [], [] ] @@ -2157,7 +2068,7 @@ def generateBaseMesh(cls, region, options): if None in nids: continue # left atrial appendage scalefactors = None - meshGroups = [ laMeshGroup ] + meshGroups = [ heartMeshGroup, laMeshGroup ] if e1 == -1: # cfb/anterior interatrial groove straddles left and right atria, collapsed to 6 node wedge nids[0] = rabNodeId[0][-elementsCountAroundAtrialSeptum] @@ -2197,7 +2108,7 @@ def generateBaseMesh(cls, region, options): if None in nids: continue # right atrial appendage scalefactors = None - meshGroups = [ raMeshGroup ] + meshGroups = [ heartMeshGroup, raMeshGroup ] # Anderson definition of right atrial appendage starts at crista terminalis: #if (e1 >= elementsCountAroundRightAtriumPosteriorVenous) and (e1 < elementsCountAroundRightAtriumFreeWall - elementsCountAroundRightAtriumAorta - 1): # meshGroups += [ raaMeshGroup ] @@ -2238,8 +2149,8 @@ def generateBaseMesh(cls, region, options): meshGroup.addElement(element) if commonLeftRightPvOstium: - # left atrium extra - meshGroups = [ laMeshGroup ] + # left atrium row above vestibule beside aorta + meshGroups = [ heartMeshGroup, laMeshGroup ] for e1 in range(elementsCountAroundLeftAtriumAorta): eft1 = eft elementtemplate1 = elementtemplate @@ -2247,11 +2158,15 @@ def generateBaseMesh(cls, region, options): lavtNodeId[1][e1], lavtNodeId[1][e1 + 1], laoaNodeId[1][e1], laoaNodeId[1][e1 + 1] ] scalefactors = None if e1 == 0: + nids[2] = asNodeId[0][elementsCountAroundAtrialSeptum + 1] + nids[6] = agNodeId[2] # general linear map d3 adjacent to interatrial groove eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, []) ]) remapEftNodeValueLabel(eft1, [ 5, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, []) ]) + remapEftNodeValueLabel(eft1, [ 4, 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 4, 8 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ) ]) scalefactors = [ -1.0 ] elementtemplateX.defineField(coordinates, -1, eft1) elementtemplate1 = elementtemplateX @@ -2265,22 +2180,30 @@ def generateBaseMesh(cls, region, options): elementIdentifier += 1 for meshGroup in meshGroups: meshGroup.addElement(element) - else: # not commonLeftRightPvOstium: - # left atrium first row of elements above appendage, anterior - meshGroups = [ laMeshGroup ] + else: + # left atrium extra row between appendage and RPV, anterior + meshGroups = [ heartMeshGroup, laMeshGroup ] for e1 in range(elementsCountAroundLeftAtriumRPV): eft1 = eft elementtemplate1 = elementtemplate - nids = [ laoaNodeId[0][e1], laoaNodeId[0][e1 + 1], lavbNodeId[0][e1], lavbNodeId[0][e1 + 1], - laoaNodeId[1][e1], laoaNodeId[1][e1 + 1], lavbNodeId[1][e1], lavbNodeId[1][e1 + 1] ] + nids = [ lavtNodeId[0][0] if (e1 == 0) else laoaNodeId[0][e1 - 1], laoaNodeId[0][e1], lavbNodeId[0][e1], lavbNodeId[0][e1 + 1], + lavtNodeId[1][0] if (e1 == 0) else laoaNodeId[1][e1 - 1], laoaNodeId[1][e1], lavbNodeId[1][e1], lavbNodeId[1][e1 + 1] ] scalefactors = None if e1 == 0: # general linear map d3 adjacent to interatrial groove eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, []) ]) remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, []) ]) remapEftNodeValueLabel(eft1, [ 5, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, []) ]) scalefactors = [ -1.0 ] + elif e1 == 1: + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 1, 5 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, []) ]) + remapEftNodeValueLabel(eft1, [ 1, 5 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, []) ]) + scalefactors = [ -1.0 ] + if eft != eft1: elementtemplateX.defineField(coordinates, -1, eft1) elementtemplate1 = elementtemplateX element = mesh.createElement(elementIdentifier, elementtemplate1) @@ -2294,37 +2217,38 @@ def generateBaseMesh(cls, region, options): for meshGroup in meshGroups: meshGroup.addElement(element) - # left atrium first row of elements above vestibule, posterior - meshGroups = [ laMeshGroup ] - scalefactors = [ -1.0 ] - elementsCount = len(lavqx[1]) - 1 if commonLeftRightPvOstium else elementsCountAroundLeftAtriumRPV - for e1 in range(elementsCount): - nc = elementsCountAroundLeftAtriumFreeWall - e1 - nids = [ lavqNodeId[0][e1], lavqNodeId[0][e1 + 1], lavtNodeId[0][nc], lavtNodeId[0][nc - 1], - lavqNodeId[1][e1], lavqNodeId[1][e1 + 1], lavtNodeId[1][nc], lavtNodeId[1][nc - 1] ] - eft1 = tricubichermite.createEftNoCrossDerivatives() - setEftScaleFactorIds(eft1, [1], []) - scaleEftNodeValueLabels(eft1, [ 3, 4, 7, 8 ], [ Node.VALUE_LABEL_D_DS1 ], [ 1 ]) - scaleEftNodeValueLabels(eft1, [ 3, 4, 7, 8 ], [ Node.VALUE_LABEL_D_DS2 ], [ 1 ]) - if e1 == 0: - # general linear map d3 adjacent to collapsed inter-atrial groove - remapEftNodeValueLabel(eft1, [ 1, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, []) ]) - remapEftNodeValueLabel(eft1, [ 3, 5 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, []) ]) - elif commonLeftRightPvOstium and (e1 == (elementsCount - 1)): - remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1]) ]) - remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, []) ]) - elementtemplateX.defineField(coordinates, -1, eft1) - elementtemplate1 = elementtemplateX - element = mesh.createElement(elementIdentifier, elementtemplate1) - result2 = element.setNodesByIdentifier(eft1, nids) - result3 = element.setScaleFactors(eft1, scalefactors) - #print('create element lavq', element.isValid(), elementIdentifier, result2, result3, nids) - elementIdentifier += 1 - for meshGroup in meshGroups: - meshGroup.addElement(element) + if not commonLeftRightPvOstium: + # left atrium first row of elements above vestibule, posterior + meshGroups = [ heartMeshGroup, laMeshGroup ] + scalefactors = [ -1.0 ] + elementsCount = len(lavqx[1]) - 1 if commonLeftRightPvOstium else elementsCountAroundLeftAtriumRPV + for e1 in range(elementsCount): + nc = elementsCountAroundLeftAtriumFreeWall - e1 + nids = [ lavqNodeId[0][e1], lavqNodeId[0][e1 + 1], lavtNodeId[0][nc], lavtNodeId[0][nc - 1], + lavqNodeId[1][e1], lavqNodeId[1][e1 + 1], lavtNodeId[1][nc], lavtNodeId[1][nc - 1] ] + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + scaleEftNodeValueLabels(eft1, [ 3, 4, 7, 8 ], [ Node.VALUE_LABEL_D_DS1 ], [ 1 ]) + scaleEftNodeValueLabels(eft1, [ 3, 4, 7, 8 ], [ Node.VALUE_LABEL_D_DS2 ], [ 1 ]) + if e1 == 0: + # general linear map d3 adjacent to collapsed inter-atrial groove + remapEftNodeValueLabel(eft1, [ 1, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, []) ]) + remapEftNodeValueLabel(eft1, [ 3, 5 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, []) ]) + elif commonLeftRightPvOstium and (e1 == (elementsCount - 1)): + remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1]) ]) + remapEftNodeValueLabel(eft1, [ 2, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, []) ]) + elementtemplateX.defineField(coordinates, -1, eft1) + elementtemplate1 = elementtemplateX + element = mesh.createElement(elementIdentifier, elementtemplate1) + result2 = element.setNodesByIdentifier(eft1, nids) + result3 = element.setScaleFactors(eft1, scalefactors) + #print('create element lavq', element.isValid(), elementIdentifier, result2, result3, nids) + elementIdentifier += 1 + for meshGroup in meshGroups: + meshGroup.addElement(element) # create atrial septum base row of elements - meshGroups = [ laMeshGroup, raMeshGroup, aSeptumMeshGroup ] + meshGroups = [ heartMeshGroup, laMeshGroup, raMeshGroup, aSeptumMeshGroup ] for e1 in range(elementsCountAroundAtrialSeptum): n1l = elementsCountAroundLeftAtriumFreeWall + e1 - elementsCountAroundLeftAtrium n1r = -e1 @@ -2355,7 +2279,7 @@ def generateBaseMesh(cls, region, options): meshGroup.addElement(element) # create fossa ovalis elements - meshGroups = [ laMeshGroup, raMeshGroup, aSeptumMeshGroup, fossaMeshGroup ] + meshGroups = [ heartMeshGroup, laMeshGroup, raMeshGroup, aSeptumMeshGroup, fossaMeshGroup ] radiansAround0 = fossaRadiansAround[-1] if radiansAround0 > fossaRadiansAround[0]: radiansAround0 -= 2.0*math.pi @@ -2391,7 +2315,7 @@ def generateBaseMesh(cls, region, options): radiansAround0, radiansAround1, radiansAround2 = radiansAround1, radiansAround2, radiansAround3 # create atrial septum elements around fossa - meshGroups = [ laMeshGroup, raMeshGroup, aSeptumMeshGroup ] + meshGroups = [ heartMeshGroup, laMeshGroup, raMeshGroup, aSeptumMeshGroup ] for e1 in range(elementsCountAroundFossa): e1p = (e1 + 1)%elementsCountAroundFossa nids = [ asNodeId[0][e1], asNodeId[0][e1p], foNodeId[0][e1], foNodeId[0][e1p], \ @@ -2448,7 +2372,7 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eftAGroove, [ 1, 2, 3, 4, 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) ln_map = [ 1, 2, 3, 4, 1, 2, 5, 6 ] remapEftLocalNodes(eftAGroove, 6, ln_map) - meshGroups = [ laMeshGroup, raMeshGroup ] + meshGroups = [ heartMeshGroup, laMeshGroup, raMeshGroup ] elementsCountOverArch = elementsCountOverAtria - 2 elementtemplateX.defineField(coordinates, -1, eftAGroove) scalefactors = [ -1.0 ] @@ -2490,7 +2414,7 @@ def generateBaseMesh(cls, region, options): if ix > 0: ix -= elementsCountAroundLpvOstium # down interatrial groove from anterior venous limit, including both corners - for n1 in range(elementsCountOverLeftAtriumVenous + 1): + for n1 in range(elementsCountOverLeftAtriumVenous + 2): ns = n1 - elementsCountOverLeftAtriumVenous - 1 ng = elementsCountOverLeftAtriumNonVenousAnterior + n1 lpvax [0][ix] = asx [0][ns] @@ -2506,26 +2430,30 @@ def generateBaseMesh(cls, region, options): if n1 == 0: lpvaDerivativesMap[0][ix] = ( (-1, 0, 0), (-1, -1, 0), (1, 0, 1), (0, 1, 0 ) ) lpvaDerivativesMap[1][ix] = ( (-1, 0, 0), (-1, -1, 0), (-1, 0, 1), (0, 1, 0 ) ) - elif n1 == elementsCountOverLeftAtriumVenous: - lpvaDerivativesMap[0][ix] = ( (0, 1, 0), (-1, 1, 0), (1, 0, 1), (1, 0, 0 ) ) - lpvaDerivativesMap[1][ix] = ( (0, 1, 0), (-1, 1, 0), (-1, 0, 1), (1, 0, 0 ) ) + elif ((n1 == elementsCountOverLeftAtriumVenous) or + ((elementsCountOverAtria > 6) and (n1 == (elementsCountOverLeftAtriumVenous - 1)))): + lpvaDerivativesMap[0][ix] = ( (0, 1, 0), (-1, 1, 0), (1, 0, 1) ) + lpvaDerivativesMap[1][ix] = ( (0, 1, 0), (-1, 1, 0), (-1, 0, 1) ) + elif n1 == (elementsCountOverLeftAtriumVenous + 1): + lpvaDerivativesMap[0][ix] = ( (0, -1, 0), (1, -1, 0), (1, 0, 1), (-1, 0, 0 ) ) + lpvaDerivativesMap[1][ix] = ( (0, -1, 0), (1, -1, 0), (-1, 0, 1), (-1, 0, 0 ) ) else: lpvaDerivativesMap[0][ix] = ( (0, 1, 0), (-1, 0, 0), (1, 0, 1) ) lpvaDerivativesMap[1][ix] = ( (0, 1, 0), (-1, 0, 0), (-1, 0, 1) ) ix += 1 - # left around posterior venous limit lavq + # left around posterior vestibule top lavt for n1 in range(1, elementsCountAroundLeftAtriumRPV + elementsCountAroundLeftAtriumLPV): - nc = n1 + nc = elementsCountAroundLeftAtriumFreeWall - n1 for n3 in range(2): - lpvax [n3][ix] = lavqx [n3][nc] - lpvad1[n3][ix] = lavqd1[n3][nc] - lpvad2[n3][ix] = lavqd2[n3][nc] - lpvad3[n3][ix] = lavqd3[n3][nc] - lpvaNodeId[n3][ix] = lavqNodeId[n3][nc] - lpvaDerivativesMap[n3][ix] = ( None, None, None ) + lpvax [n3][ix] = lavtx [n3][nc] + lpvad1[n3][ix] = lavtd1[n3][nc] + lpvad2[n3][ix] = lavtd2[n3][nc] + lpvad3[n3][ix] = lavtd3[n3][nc] + lpvaNodeId[n3][ix] = lavtNodeId[n3][nc] + lpvaDerivativesMap[n3][ix] = ( (-1, 0, 0), (0, -1, 0), None ) ix += 1 # up left atrium over appendage laoa to interatrial groove - for n1 in range(1, elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV): + for n1 in range(laoaCount - 1): no = -1 - n1 for n3 in range(2): lpvax [n3][ix] = laoax [n3][no] @@ -2533,8 +2461,10 @@ def generateBaseMesh(cls, region, options): lpvad2[n3][ix] = laoad2[n3][no] lpvad3[n3][ix] = laoad3[n3][no] lpvaNodeId[n3][ix] = laoaNodeId[n3][no] - if n1 == 1: - lpvaDerivativesMap[n3][ix] = ( (0, -1, 0), (1, -1, 0), None, (-1, 0, 0 ) ) + if n1 == 0: + lpvaDerivativesMap[n3][ix] = ( (-1, 0, 0), (-1, -1, 0), None, (0, 1, 0 ) ) + elif n1 == (laoaCount - 2): + lpvaDerivativesMap[n3][ix] = ( (-1, 0, 0), (-1, -1, 0), None, (0, 1, 0 ) ) else: lpvaDerivativesMap[n3][ix] = ( (-1, 0, 0), (0, -1, 0), None ) ix += 1 @@ -2543,19 +2473,14 @@ def generateBaseMesh(cls, region, options): # insert at indexes such that 0 is one past the midpoint on venous midline ix = -elementsCountOverLeftAtriumVenous//2 # down left atrium venous midpoint line - for n1 in range(elementsCountOverLeftAtriumVenous + 1): + for n1 in range(1, elementsCountOverLeftAtriumVenous + 2): for n3 in range(2): lpvax [n3][ix] = lamlx [n3][n1] lpvad1[n3][ix] = lamld1[n3][n1] lpvad2[n3][ix] = lamld2[n3][n1] lpvad3[n3][ix] = lamld3[n3][n1] lpvaNodeId[n3][ix] = lamlNodeId[n3][n1] - if n1 == 0: - lpvaDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, 0, 0), None, (1, 1, 0 ) ) - elif n1 == elementsCountOverLeftAtriumVenous: - lpvaDerivativesMap[n3][ix] = ( (-1, 1, 0), (-1, 0, 0), None, (0, 1, 0 ) ) - else: - lpvaDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, 0, 0), None ) + lpvaDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, 0, 0), None ) ix += 1 # left around left cs, including 2 corners for n1 in range(elementsCountAroundLeftAtriumLPV + 1): @@ -2591,11 +2516,36 @@ def generateBaseMesh(cls, region, options): #print('lpvaNodeId[1]',lpvaNodeId[1]) #print('lpvaDerivativesMap[0]',lpvaDerivativesMap[0]) #print('lpvaDerivativesMap[1]',lpvaDerivativesMap[1]) + #if commonLeftRightPvOstium: + lpvoProportions = [ list(laTrackSurface.getProportion(lpvoPositions[i])) for i in range(elementsCountAroundLpvOstium) ] + lpvaProportions = [ list(laTrackSurface.getProportion(laTrackSurface.findNearestPosition(lpvax[1][i], startPosition=lpvoPositions[i]))) + for i in range(elementsCountAroundLpvOstium) ] nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, lpvox, lpvod1, lpvod2, lpvod3, lpvoNodeId, None, lpvax, lpvad1, lpvad2, lpvad3, lpvaNodeId, lpvaDerivativesMap, - elementsCountRadial = 1, meshGroups = [ laMeshGroup ]) + maxEndThickness=laVenousFreeWallThickness, + elementsCountRadial = elementsCountRadialPVAnnuli, meshGroups = [ heartMeshGroup, laMeshGroup ], + tracksurface=laTrackSurface, startProportions=lpvoProportions, endProportions=lpvaProportions, + rescaleStartDerivatives=True, rescaleEndDerivatives=True, sampleBlend=0.0) + + # left atrium epicardium venous midpoint marker point + if commonLeftRightPvOstium: + laevmElementId = elementIdentifier - elementsCountAroundLpvOstium*elementsCountRadialPVAnnuli + (elementsCountAroundLpvOstium // 2); + if lpvOstiumSettings['Number of vessels'] == 3: + laevmElementId += 1 + laevmXi = [ 0.0, 0.0, 1.0 ] + else: + laevmElementId = elementIdentifier - elementsCountAroundLpvOstium + laevmXi = [ 0.0, 1.0, 1.0 ] + laevmElement = mesh.findElementByIdentifier(laevmElementId) + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, laeVenousMidpointGroup.getName()) + markerLocation.assignMeshLocation(cache, laevmElement, laevmXi) + for group in [ heartGroup, laGroup, laeVenousMidpointGroup ]: + group.getNodesetGroup(nodes).addNode(markerPoint) if not commonLeftRightPvOstium: # create right pulmonary vein annulus @@ -2644,7 +2594,7 @@ def generateBaseMesh(cls, region, options): ix += 1 # up left atrium venous midline, including both corners for n1 in range(elementsCountOverLeftAtriumVenous + 1): - nm = elementsCountOverLeftAtriumVenous - n1 + nm = elementsCountOverLeftAtriumVenous + 1 - n1 for n3 in range(2): rpvax [n3][ix] = lamlx [n3][nm] rpvad1[n3][ix] = lamld1[n3][nm] @@ -2652,9 +2602,9 @@ def generateBaseMesh(cls, region, options): rpvad3[n3][ix] = lamld3[n3][nm] rpvaNodeId[n3][ix] = lamlNodeId[n3][nm] if n1 == 0: - rpvaDerivativesMap[n3][ix] = ( (1, 0, 0), (0, 1, 0), None, (1, -1, 0) ) + rpvaDerivativesMap[n3][ix] = ( (1, 0, 0), (1, 1, 0), None, (0, -1, 0) ) elif n1 == elementsCountOverLeftAtriumVenous: - rpvaDerivativesMap[n3][ix] = ( (-1, -1, 0), (0, -1, 0), None, (-1, 0, 0) ) + rpvaDerivativesMap[n3][ix] = ( (0, -1, 0), (1, -1, 0), None, (-1, 0, 0) ) else: rpvaDerivativesMap[n3][ix] = ( (0, -1, 0), (1, 0, 0), None ) ix += 1 @@ -2673,11 +2623,17 @@ def generateBaseMesh(cls, region, options): #print('rpvaNodeId[1]',rpvaNodeId[1]) #print('rpvaDerivativesMap[0]',rpvaDerivativesMap[0]) #print('rpvaDerivativesMap[1]',rpvaDerivativesMap[1]) + rpvoProportions = [ list(laTrackSurface.getProportion(rpvoPositions[i])) for i in range(elementsCountAroundRpvOstium) ] + rpvaProportions = [ list(laTrackSurface.getProportion(laTrackSurface.findNearestPosition(rpvax[1][i], startPosition=rpvoPositions[i]))) + for i in range(elementsCountAroundRpvOstium) ] nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, rpvox, rpvod1, rpvod2, rpvod3, rpvoNodeId, None, rpvax, rpvad1, rpvad2, rpvad3, rpvaNodeId, rpvaDerivativesMap, - elementsCountRadial = 1, meshGroups = [ laMeshGroup ]) + maxEndThickness=laVenousFreeWallThickness, + elementsCountRadial = elementsCountRadialPVAnnuli, meshGroups = [ heartMeshGroup, laMeshGroup ], + tracksurface=laTrackSurface, startProportions=rpvoProportions, endProportions=rpvaProportions, + rescaleStartDerivatives=True, rescaleEndDerivatives=True, sampleBlend=0.0) # create inferior and superior vena cavae inlets for v in range(2): @@ -2730,9 +2686,9 @@ def generateBaseMesh(cls, region, options): vcaNodeId = [ [ None ]*elementsCountAroundVC, [ None ]*elementsCountAroundVC ] vcaDerivativesMap = [ [ None ]*elementsCountAroundVC, [ None ]*elementsCountAroundVC ] if v == 0: # ivc - # set points clockwise from interatrial groove at venous midpoint + # set points clockwise starting above crux of heart ix = 0 - # over interatrial groove from cristvenous midpoint to anterior, including both corners + # up interatrial groove to venous midpoint, including both corners for n1 in range(elementsCountOverRightAtriumVenous//2, -1, -1): ns = (elementsCountAroundAtrialSeptum - 1 + elementsCountOverRightAtriumNonVenousAnterior + elementsCountOverRightAtriumVenous//2 + n1) % elementsCountAroundFossa ng = elementsCountOverRightAtriumNonVenousAnterior + elementsCountOverRightAtriumVenous//2 + n1 @@ -2867,9 +2823,17 @@ def generateBaseMesh(cls, region, options): vcMeshGroup = ivcMeshGroup if (v == 0) else svcMeshGroup vcInletMeshGroup = ivcInletMeshGroup if (v == 0) else svcInletMeshGroup if elementsCountAlongVCInlet == 1: - rowMeshGroups = [ [ vcMeshGroup, vcInletMeshGroup, raMeshGroup] ] + rowMeshGroups = [ [ heartMeshGroup, vcMeshGroup, vcInletMeshGroup, raMeshGroup] ] else: - rowMeshGroups = [ [ vcMeshGroup ] ]*max(0, (elementsCountAlongVCInlet - 2)) + [ [ vcMeshGroup, vcInletMeshGroup ], [ vcInletMeshGroup, raMeshGroup] ] + rowMeshGroups = [] + for i in range(elementsCountAlongVCInlet): + xi = (i + 1)/elementsCountAlongVCInlet + meshGroups = [ heartMeshGroup ] + if xi < 0.67: + meshGroups += [ vcMeshGroup ] + if xi > 0.51: + meshGroups += [ vcInletMeshGroup, raMeshGroup ] + rowMeshGroups.append(meshGroups) nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, vcvx, vcvd1, vcvd2, None, None, None, @@ -2878,6 +2842,19 @@ def generateBaseMesh(cls, region, options): elementsCountRadial = elementsCountAlongVCInlet, meshGroups = rowMeshGroups) + if v == 0: # ivc + # right atrium epicardium venous midpoint marker point + raevmElementId = elementIdentifier - elementsCountAroundVC + elementsCountOverRightAtriumVenous//2 + elementsCountOverSideRightAtriumVC//2 + raevmXi = [ 0.5 if (elementsCountOverSideRightAtriumVC % 2) else 0.0, 1.0, 1.0 ] + raevmElement = mesh.findElementByIdentifier(raevmElementId) + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, raeVenousMidpointGroup.getName()) + markerLocation.assignMeshLocation(cache, raevmElement, raevmXi) + for group in [ heartGroup, raGroup, raeVenousMidpointGroup ]: + group.getNodesetGroup(nodes).addNode(markerPoint) + # create left atrial appendage position = laTrackSurface.createPositionProportion(laaMidpointOver, laaMidpointLeft) laamx, d1, d2 = laTrackSurface.evaluateCoordinates(position, derivatives = True) @@ -2893,10 +2870,7 @@ def generateBaseMesh(cls, region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, laamd2) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, laamd3) nodeIdentifier += 1 - elementsCountAroundLaa = elementsCountAroundLeftAtrialAppendageBase + elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV + 1 - if commonLeftRightPvOstium: - elementsCountAroundLaa += 1 - #print('elementsCountAroundLaa', elementsCountAroundLaa) + elementsCountAroundLaa = elementsCountAroundLeftAtrialAppendageBase + len(laoaProportions) + 1 # get start points, nodes, derivative maps laasx = [ [ None ]*elementsCountAroundLaa, [ None ]*elementsCountAroundLaa ] laasd1 = [ [ None ]*elementsCountAroundLaa, [ None ]*elementsCountAroundLaa ] @@ -2935,8 +2909,8 @@ def generateBaseMesh(cls, region, options): else: laasDerivativesMap[n3][ix] = ( None, None, None ) ix += 1 - # right over appendage laoa - for n1 in range(elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV): + # left over appendage laoa + for n1 in range(laoaCount): no = -1 - n1 for n3 in range(2): laasx [n3][ix] = laoax [n3][no] @@ -2946,8 +2920,8 @@ def generateBaseMesh(cls, region, options): laasNodeId[n3][ix] = laoaNodeId[n3][no] if n1 == 0: laasDerivativesMap[n3][ix] = ( (0, 1, 0), (-1, 0, 0), None ) - elif n1 == (elementsCountAroundLeftAtriumRPV + elementsCountOverSideLeftAtriumLPV - 1): - laasDerivativesMap[n3][ix] = ( (-1, 0, 0), (1, -1, 0), None, (0, -1, 0 ) ) + elif n1 == (laoaCount - 1): + laasDerivativesMap[n3][ix] = ( (0, -1, 0), (1, 0, 0), None ) else: laasDerivativesMap[n3][ix] = ( (-1, 0, 0), (0, -1, 0), None ) ix += 1 @@ -2955,11 +2929,10 @@ def generateBaseMesh(cls, region, options): #print('laasNodeId[1]',laasNodeId[1]) #print('laasDerivativesMap[0]',laasDerivativesMap[0]) #print('laasDerivativesMap[1]',laasDerivativesMap[1]) - elementsCountLaaRadial = 2 # get end points, nodes, derivative maps, expanding from wedge laawx, laawd1, laawd2, laawd3, elementsCountAcrossLaaWedge, laawPointsMap, laaeDerivativesMap = \ getAtrialAppendageWedgePoints(laamx, laamd1, laamd2, laamd3, laaAngleLeftRadians, laaAngleUpradians, laaAngleAxialRadians, laaBaseLength, - elementsCountAroundLaa, elementsCountLaaRadial, laaArcLength, laaArcRadius, laaWallThickness, laaWedgeAngleRadians) + elementsCountAroundLaa, elementsCountAlongAtrialAppendages, laaArcLength, laaArcRadius, laaWallThickness, laaWedgeAngleRadians) # create laa wedge nodes: laawNodeId = [ [], [] ] for n3 in range(2): @@ -2995,13 +2968,13 @@ def generateBaseMesh(cls, region, options): laaex, laaed1, laaed2, laaed3, laaeNodeId, laaeDerivativesMap, forceMidLinearXi3 = True, forceEndLinearXi3 = True, maxStartThickness = laaWallThickness, - elementsCountRadial = elementsCountLaaRadial, - meshGroups = [ laMeshGroup, laaMeshGroup ]) + elementsCountRadial = elementsCountAlongAtrialAppendages, + meshGroups = [ heartMeshGroup, laMeshGroup, laaMeshGroup ]) # create right atrium plain elements # Anderson considers these part of the right atrial appendage: #meshGroups = [ raMeshGroup, raaMeshGroup ] - meshGroups = [ raMeshGroup ] + meshGroups = [ heartMeshGroup, raMeshGroup ] for e2 in range(2): for e1 in range(elementsCountOverSideRightAtriumPouch): eft1 = eft @@ -3149,11 +3122,10 @@ def generateBaseMesh(cls, region, options): #print('raasNodeId[1]',raasNodeId[1]) #print('raasDerivativesMap[0]',raasDerivativesMap[0]) #print('raasDerivativesMap[1]',raasDerivativesMap[1]) - elementsCountRaaRadial = 2 # get end points, nodes, derivative maps, expanding from wedge raawx, raawd1, raawd2, raawd3, elementsCountAcrossRaaWedge, raawPointsMap, raaeDerivativesMap = \ getAtrialAppendageWedgePoints(raamx, raamd1, raamd2, raamd3, raaAngleLeftRadians, raaAngleUpradians, raaAngleAxialRadians, raaBaseLength, - elementsCountAroundRaa, elementsCountRaaRadial, raaArcLength, raaArcRadius, raaWallThickness, raaWedgeAngleRadians) + elementsCountAroundRaa, elementsCountAlongAtrialAppendages, raaArcLength, raaArcRadius, raaWallThickness, raaWedgeAngleRadians) # create raa wedge nodes: raawNodeId = [ [], [] ] for n3 in range(2): @@ -3189,8 +3161,8 @@ def generateBaseMesh(cls, region, options): raaex, raaed1, raaed2, raaed3, raaeNodeId, raaeDerivativesMap, forceMidLinearXi3 = True, forceEndLinearXi3 = True, maxStartThickness = raaWallThickness, - elementsCountRadial = elementsCountRaaRadial, - meshGroups = [ raMeshGroup, raaMeshGroup ]) + elementsCountRadial = elementsCountAlongAtrialAppendages, + meshGroups = [ heartMeshGroup, raMeshGroup, raaMeshGroup ]) if drawLaTrackSurface: mesh2d = fm.findMeshByDimension(2) @@ -3351,7 +3323,7 @@ def getLeftAtriumPulmonaryVeinOstiaElementsCounts(elementsCountAroundLeftAtriumF = getLeftAtriumBaseFreewallElementsCounts(elementsCountAroundLeftAtriumFreeWall) elementsCountAroundRpvOstium = 2*(elementsCountOverLeftAtriumVenous + elementsCountAroundLeftAtriumRPV) if commonLeftRightPvOstium: - elementsCountAroundLpvOstium = elementsCountOverLeftAtriumVenous + 2*(elementsCountAroundLeftAtriumLPV + elementsCountAroundLeftAtriumRPV) + elementsCountAroundLpvOstium = elementsCountOverLeftAtriumVenous + 2*(elementsCountAroundLeftAtriumLPV + elementsCountAroundLeftAtriumRPV) + 2 else: elementsCountAroundLpvOstium = elementsCountOverLeftAtriumVenous + 2 + 2*elementsCountAroundLeftAtriumLPV return elementsCountAroundLpvOstium, elementsCountAroundRpvOstium diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py index 9c9344f3..ae4c687e 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles1.py @@ -246,10 +246,12 @@ def generateBaseMesh(cls, region, options): mesh = fm.findMeshByDimension(3) + heartGroup = AnnotationGroup(region, get_heart_term("heart")) + apexGroup = AnnotationGroup(region, get_heart_term("apex of heart")) lvGroup = AnnotationGroup(region, get_heart_term("left ventricle myocardium")) rvGroup = AnnotationGroup(region, get_heart_term("right ventricle myocardium")) vSeptumGroup = AnnotationGroup(region, get_heart_term("interventricular septum")) - annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ] + annotationGroups = [ heartGroup, apexGroup, lvGroup, rvGroup, vSeptumGroup ] # annotation fiducial points markerGroup = findOrCreateFieldGroup(fm, "marker") @@ -702,6 +704,7 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + heartMeshGroup = heartGroup.getMeshGroup(mesh) lvMeshGroup = lvGroup.getMeshGroup(mesh) rvMeshGroup = rvGroup.getMeshGroup(mesh) vSeptumMeshGroup = vSeptumGroup.getMeshGroup(mesh) @@ -725,7 +728,7 @@ def generateBaseMesh(cls, region, options): # LV apex - meshGroups = [ lvMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup ] for e1 in range(elementsCountAroundLV): eft1 = eft @@ -777,7 +780,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft scalefactors = None - meshGroups = [ lvMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup ] va = e1 if (e1 >= 0) else (elementsCountAroundLV - 1) vb = e1 + 1 @@ -829,7 +832,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft scalefactors = None - meshGroups = [ rvMeshGroup ] + meshGroups = [ heartMeshGroup, rvMeshGroup ] ua = e1 ub = e1 + 1 va = elementsCountAroundLVFreeWall + e1 @@ -945,7 +948,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft scalefactors = None - meshGroups = [ lvMeshGroup, rvMeshGroup, vSeptumMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup, rvMeshGroup, vSeptumMeshGroup ] ua = -e1 ub = ua - 1 @@ -1023,8 +1026,10 @@ def generateBaseMesh(cls, region, options): markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerName.assignString(cache, 'APEX') # interlex.org has 'apex of heart' + markerName.assignString(cache, apexGroup.getName()) markerLocation.assignMeshLocation(cache, element1, [ 0.0, 0.0, 1.0 ]) + for group in [ heartGroup, lvGroup, apexGroup ]: + group.getNodesetGroup(nodes).addNode(markerPoint) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py index 0a22b329..bcf83d9f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py @@ -337,10 +337,12 @@ def generateBaseMesh(cls, region, options): nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) mesh = fieldmodule.findMeshByDimension(3) + heartGroup = AnnotationGroup(region, get_heart_term("heart")) + apexGroup = AnnotationGroup(region, get_heart_term("apex of heart")) lvGroup = AnnotationGroup(region, get_heart_term("left ventricle myocardium")) rvGroup = AnnotationGroup(region, get_heart_term("right ventricle myocardium")) vSeptumGroup = AnnotationGroup(region, get_heart_term("interventricular septum")) - annotationGroups = [ lvGroup, rvGroup, vSeptumGroup ] + annotationGroups = [ heartGroup, apexGroup, lvGroup, rvGroup, vSeptumGroup ] # annotation fiducial points markerGroup = findOrCreateFieldGroup(fieldmodule, "marker") @@ -924,22 +926,23 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + heartMeshGroup = heartGroup.getMeshGroup(mesh) lvMeshGroup = lvGroup.getMeshGroup(mesh) rvMeshGroup = rvGroup.getMeshGroup(mesh) vSeptumMeshGroup = vSeptumGroup.getMeshGroup(mesh) elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) #print("LV elements") - elementIdentifier = lvShield.generateElements(fieldmodule, coordinates, elementIdentifier, [ lvMeshGroup ]) + elementIdentifier = lvShield.generateElements(fieldmodule, coordinates, elementIdentifier, [ heartMeshGroup, lvMeshGroup ]) #print("RV elements") - elementIdentifier = rvShield.generateElements(fieldmodule, coordinates, elementIdentifier, [ rvMeshGroup ]) + elementIdentifier = rvShield.generateElements(fieldmodule, coordinates, elementIdentifier, [ heartMeshGroup, rvMeshGroup ]) tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) elementtemplate1 = mesh.createElementtemplate() elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) # interventricular sulcus elements - meshGroups = [ lvMeshGroup, rvMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup, rvMeshGroup ] lvNodeId = lvShield.nodeId rvNodeId = rvShield.nodeId n1ln, n2ln = lvShield.convertRimIndex(0) @@ -1031,8 +1034,10 @@ def generateBaseMesh(cls, region, options): markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerName.assignString(cache, 'APEX') # interlex.org has 'apex of heart' + markerName.assignString(cache, apexGroup.getName()) markerLocation.assignMeshLocation(cache, apexElement, [ 0.0, 0.0, 1.0 ]) + for group in [ heartGroup, lvGroup, apexGroup ]: + group.getNodesetGroup(nodes).addNode(markerPoint) return annotationGroups diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py index 0bba85c6..5523a558 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventriclesbase1.py @@ -300,6 +300,7 @@ def generateBaseMesh(cls, region, options): annotationGroups = MeshType_3d_heartventricles1.generateBaseMesh(region, options) # find/add annotation groups + heartGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("heart")) lvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("left ventricle myocardium")) rvGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("right ventricle myocardium")) vSeptumGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("interventricular septum")) @@ -726,6 +727,7 @@ def generateBaseMesh(cls, region, options): # Create elements ################# + heartMeshGroup = heartGroup.getMeshGroup(mesh) lvMeshGroup = lvGroup.getMeshGroup(mesh) rvMeshGroup = rvGroup.getMeshGroup(mesh) vSeptumMeshGroup = vSeptumGroup.getMeshGroup(mesh) @@ -739,14 +741,14 @@ def generateBaseMesh(cls, region, options): elementtemplate1 = mesh.createElementtemplate() elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - scalefactors5hanging = [ -1.0, 0.5, 0.25, 0.125, 0.75 ] + scalefactors4hanging = [ -1.0, 0.5, 0.125, 0.75 ] # LV base elements row 1, starting at anterior interventricular sulcus for e in range(-1, elementsCountAroundLVFreeWall): eft1 = eft nids = None scalefactors = None - meshGroups = [ lvMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup ] if e == -1: # 4 node collapsed tetrahedral element on anterior interventricular sulcus @@ -810,7 +812,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft nids = None scalefactors = None - meshGroups = [ rvMeshGroup ] + meshGroups = [ heartMeshGroup, rvMeshGroup ] addMarker = None noa = e @@ -850,16 +852,16 @@ def generateBaseMesh(cls, region, options): remapEftNodeValueLabel(eft1, [ 5, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, []) ]) elif elementsCountRVHanging: eft1 = tricubichermite.createEftNoCrossDerivatives() - setEftScaleFactorIds(eft1, [1, 102, 104, 108, 304], []) - scalefactors = scalefactors5hanging + setEftScaleFactorIds(eft1, [1, 102, 108, 304], []) + scalefactors = scalefactors4hanging if e in [ 1, 3 ]: # 1st of pair of elements with hanging nodes at xi1=0.5 on xi2 == 0 plane - tricubichermite.setEftMidsideXi1HangingNode(eft1, 2, 1, 1, 2, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi1HangingNode(eft1, 6, 5, 5, 6, [1, 2, 3, 4, 5]) + tricubichermite.setEftMidsideXi1HangingNode(eft1, 2, 1, 1, 2, [1, 2, 3, 4]) + tricubichermite.setEftMidsideXi1HangingNode(eft1, 6, 5, 5, 6, [1, 2, 3, 4]) else: # e in [ 2, 4 ]: # 2nd of pair of elements with hanging nodes at xi1=0.5 on xi2 == 0 plane - tricubichermite.setEftMidsideXi1HangingNode(eft1, 1, 2, 1, 2, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi1HangingNode(eft1, 5, 6, 5, 6, [1, 2, 3, 4, 5]) + tricubichermite.setEftMidsideXi1HangingNode(eft1, 1, 2, 1, 2, [1, 2, 3, 4]) + tricubichermite.setEftMidsideXi1HangingNode(eft1, 5, 6, 5, 6, [1, 2, 3, 4]) if e == 4: remapEftNodeValueLabel(eft1, [ 4, 8 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, []) ]) elif e == elementsCountRVFreeWallRegular: @@ -970,7 +972,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft nids = None scalefactors = None - meshGroups = [ lvMeshGroup, rvMeshGroup, vSeptumMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup, rvMeshGroup, vSeptumMeshGroup ] lv1 = elementsCountAroundLVFreeWall + e lv2 = (lv1 + 1)%elementsCountAroundLV @@ -1073,7 +1075,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft nids = None scalefactors = None - meshGroups = [ lvMeshGroup ] + meshGroups = [ heartMeshGroup, lvMeshGroup ] addMarker = None eft1 = tricubichermite.createEftNoCrossDerivatives() @@ -1196,7 +1198,7 @@ def generateBaseMesh(cls, region, options): eft1 = eft nids = None scalefactors = None - meshGroups = [ rvMeshGroup ] + meshGroups = [ heartMeshGroup, rvMeshGroup ] if e == 0: # 6 node collapsed vs-ra shim element diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 36515fce..087337d3 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -1309,15 +1309,14 @@ def setEftMidsideXi1HangingNode(self, eft, hangingBasisNode, otherBasisNode, loc from localNode1 at xi1=0 and localNode2 at xi1=1 to midside xi1 = 0.5. Note! Cross derivatives are not handled and are currently unmodified. :param otherBasisNode: Other node along xi1 which needs its dxi1 derivative halved - :param scaleFactorIndexes: Local scale factor indexes for general values -1.0 0.5 0.25 0.125 0.75 + :param scaleFactorIndexes: Local scale factor indexes for general values -1.0 0.5 0.125 0.75 ''' n = hangingBasisNode - 1 o = otherBasisNode - 1 sfneg1 = scaleFactorIndexes[0] sf05 = scaleFactorIndexes[1] - sf025 = scaleFactorIndexes[2] - sf0125 = scaleFactorIndexes[3] - sf075 = scaleFactorIndexes[4] + sf0125 = scaleFactorIndexes[2] + sf075 = scaleFactorIndexes[3] # otherBasisNode d/dxi1 must be halved eft.setTermScaling(o*8 + 2, 1, [sf05]) # workaround for Zinc limitation where faces are not found due to only first @@ -1363,15 +1362,14 @@ def setEftMidsideXi3HangingNode(self, eft, hangingBasisNode, otherBasisNode, loc from localNode1 at xi3=0 and localNode2 at xi3=1 to midside xi3 = 0.5. Note! Cross derivatives are not handled and are currently unmodified. :param otherBasisNode: Other node along xi3 which needs its dxi3 derivative halved - :param scaleFactorIndexes: Local scale factor indexes for general values -1.0 0.5 0.25 0.125 0.75 + :param scaleFactorIndexes: Local scale factor indexes for general values -1.0 0.5 0.125 0.75 ''' n = hangingBasisNode - 1 o = otherBasisNode - 1 sfneg1 = scaleFactorIndexes[0] sf05 = scaleFactorIndexes[1] - sf025 = scaleFactorIndexes[2] - sf0125 = scaleFactorIndexes[3] - sf075 = scaleFactorIndexes[4] + sf0125 = scaleFactorIndexes[2] + sf075 = scaleFactorIndexes[3] # otherBasisNode d/dxi3 must be halved eft.setTermScaling(o*8 + 5, 1, [sf05]) # workaround for Zinc limitation where faces are not found due to only first diff --git a/src/scaffoldmaker/utils/meshrefinement.py b/src/scaffoldmaker/utils/meshrefinement.py index 08c36e22..38ba5871 100644 --- a/src/scaffoldmaker/utils/meshrefinement.py +++ b/src/scaffoldmaker/utils/meshrefinement.py @@ -46,6 +46,7 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): minimumsField = None maximumsField = None self._sourceMesh = self._sourceFm.findMeshByDimension(3) + self._sourceNodes = self._sourceFm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) self._sourceElementiterator = self._sourceMesh.createElementiterator() self._octree = Octree(minimums, maximums) @@ -72,12 +73,13 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): self._sourceAnnotationGroups = sourceAnnotationGroups self._annotationGroups = [] self._sourceAndTargetMeshGroups = [] + self._sourceAndTargetNodesetGroups = [] for sourceAnnotationGroup in sourceAnnotationGroups: - sourceMeshGroup = sourceAnnotationGroup.getMeshGroup(self._sourceMesh) targetAnnotationGroup = AnnotationGroup(self._targetRegion, sourceAnnotationGroup.getTerm()) - targetMeshGroup = targetAnnotationGroup.getMeshGroup(self._targetMesh) self._annotationGroups.append(targetAnnotationGroup) - self._sourceAndTargetMeshGroups.append( ( sourceMeshGroup, targetMeshGroup) ) + self._sourceAndTargetMeshGroups.append( ( sourceAnnotationGroup.getMeshGroup(self._sourceMesh), targetAnnotationGroup.getMeshGroup(self._targetMesh)) ) + self._sourceAndTargetNodesetGroups.append( ( sourceAnnotationGroup.getNodesetGroup(self._sourceNodes), targetAnnotationGroup.getNodesetGroup(self._targetNodes)) ) + # prepare element -> marker point list map self.elementMarkerMap = {} sourceMarkerGroup = findOrCreateFieldGroup(self._sourceFm, "marker") @@ -96,7 +98,7 @@ def __init__(self, sourceRegion, targetRegion, sourceAnnotationGroups = []): if not markerList: markerList = [] self.elementMarkerMap[elementIdentifier] = markerList - markerList.append( ( markerName, xi ) ) + markerList.append( ( markerName, xi, node.getIdentifier() ) ) node = nodeIter.next() if self.elementMarkerMap: self._targetMarkerGroup = findOrCreateFieldGroup(self._targetFm, "marker") @@ -203,7 +205,8 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n elementOffset = [ 1, numberInXi1, numberInXi1*numberInXi2 ] targetXi = [ 0.0 ]*3 for marker in markerList: - markerName, sourceXi = marker + markerName, sourceXi, sourceNodeIdentifier = marker + sourceNode = self._sourceNodes.findNodeByIdentifier(sourceNodeIdentifier) node = self._targetMarkerNodes.createNode(self._nodeIdentifier, self._targetMarkerTemplate) self._targetCache.setNode(node) self._targetMarkerName.assignString(self._targetCache, markerName) @@ -221,6 +224,10 @@ def refineElementCubeStandard3d(self, sourceElement, numberInXi1, numberInXi2, n targetElement = self._targetMesh.findElementByIdentifier(targetElementIdentifier) result = self._targetMarkerLocation.assignMeshLocation(self._targetCache, targetElement, targetXi) self._nodeIdentifier += 1 + # add new node to matching annotation groups the previous one was in + for sourceAndTargetNodesetGroup in self._sourceAndTargetNodesetGroups: + if sourceAndTargetNodesetGroup[0].containsNode(sourceNode): + sourceAndTargetNodesetGroup[1].addNode(node) return nids, nx diff --git a/tests/test_heart.py b/tests/test_heart.py index 3e6a4424..1e503010 100644 --- a/tests/test_heart.py +++ b/tests/test_heart.py @@ -22,16 +22,22 @@ def test_heart1(self): self.assertEqual(parameterSetNames, [ "Default", "Human 1", "Mouse 1", "Pig 1", "Rat 1", "Unit Human 1", "Unit Mouse 1", "Unit Pig 1", "Unit Rat 1" ]); options = scaffold.getDefaultOptions("Human 1") - self.assertEqual(117, len(options)) + self.assertEqual(119, len(options)) self.assertEqual(0.9, options.get("LV outer height")) self.assertEqual(80.0, options.get("Unit scale")) self.assertEqual(7, options.get("Number of elements around LV free wall")) self.assertEqual(7, options.get("Number of elements around RV free wall")) + # simplify atria + self.assertEqual(8, options.get("Number of elements over atria")) + options["Number of elements over atria"] = 6 + self.assertEqual(2, options.get("Number of elements radial pulmonary vein annuli")) + options["Number of elements radial pulmonary vein annuli"] = 1 + self.assertFalse(scaffold.checkOptions(options)) context = Context("Test") region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = scaffold.generateMesh(region, options) - self.assertEqual(27, len(annotationGroups)) + self.assertEqual(32, len(annotationGroups)) fieldmodule = region.getFieldmodule() mesh3d = fieldmodule.findMeshByDimension(3) self.assertEqual(289, mesh3d.getSize()) @@ -40,7 +46,7 @@ def test_heart1(self): mesh1d = fieldmodule.findMeshByDimension(1) self.assertEqual(1356, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(528, nodes.getSize()) + self.assertEqual(530, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -63,10 +69,10 @@ def test_heart1(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 36529.60472245443, delta=1.0E-3) + self.assertAlmostEqual(surfaceArea, 36508.985114306226, delta=1.0E-2) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 221442.64891709128, delta=1.0E-3) + self.assertAlmostEqual(volume, 221285.56778831664, delta=1.0E-2) # check some annotationGroups: expectedSizes3d = { @@ -93,15 +99,34 @@ def test_heart1(self): size = group.getMeshGroup(mesh2d).getSize() self.assertEqual(expectedSizes2d[name], size, name) + # test finding a marker in scaffold + markerGroup = fieldmodule.findFieldByName("marker").castGroup() + markerNodes = markerGroup.getFieldNodeGroup(nodes).getNodesetGroup() + self.assertEqual(7, markerNodes.getSize()) + markerName = fieldmodule.findFieldByName("marker_name") + self.assertTrue(markerName.isValid()) + markerLocation = fieldmodule.findFieldByName("marker_location") + self.assertTrue(markerLocation.isValid()) + # test apex marker point + cache = fieldmodule.createFieldcache() + node = findNodeWithName(markerNodes, markerName, "apex of heart") + self.assertTrue(node.isValid()) + cache.setNode(node) + element, xi = markerLocation.evaluateMeshLocation(cache, 3) + self.assertEqual(1, element.getIdentifier()) + assertAlmostEqualList(self, xi, [ 0.0, 0.0, 1.0 ], 1.0E-10) + apexGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("apex of heart")) + self.assertTrue(apexGroup.getNodesetGroup(nodes).containsNode(node)) + # refine 2x2x2 and check result - # first remove any surface annotation groups as they are re-added by defineFaceAnnotations + # first remove any face (but not point) annotation groups as they are re-added by defineFaceAnnotations removeAnnotationGroups = [] for annotationGroup in annotationGroups: - if annotationGroup.getMeshGroup(mesh3d).getSize() == 0: + if (not annotationGroup.hasMeshGroup(mesh3d)) and (annotationGroup.hasMeshGroup(mesh2d) or annotationGroup.hasMeshGroup(mesh1d)): removeAnnotationGroups.append(annotationGroup) for annotationGroup in removeAnnotationGroups: annotationGroups.remove(annotationGroup) - self.assertEqual(18, len(annotationGroups)) + self.assertEqual(23, len(annotationGroups)) refineRegion = region.createRegion() refineFieldmodule = refineRegion.getFieldmodule() @@ -120,7 +145,7 @@ def test_heart1(self): for annotation in annotationGroups: if annotation not in oldAnnotationGroups: annotationGroup.addSubelements() - self.assertEqual(27, len(annotationGroups)) + self.assertEqual(32, len(annotationGroups)) mesh3d = refineFieldmodule.findMeshByDimension(3) self.assertEqual(2236, mesh3d.getSize()) @@ -129,7 +154,7 @@ def test_heart1(self): mesh1d = refineFieldmodule.findMeshByDimension(1) self.assertEqual(8623, mesh1d.getSize()) nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(3183, nodes.getSize()) + self.assertEqual(3185, nodes.getSize()) datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -153,18 +178,21 @@ def test_heart1(self): markerGroup = refineFieldmodule.findFieldByName("marker").castGroup() refinedNodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) markerNodes = markerGroup.getFieldNodeGroup(refinedNodes).getNodesetGroup() - self.assertEqual(5, markerNodes.getSize()) + self.assertEqual(7, markerNodes.getSize()) markerName = refineFieldmodule.findFieldByName("marker_name") self.assertTrue(markerName.isValid()) markerLocation = refineFieldmodule.findFieldByName("marker_location") self.assertTrue(markerLocation.isValid()) + # test apex marker point cache = refineFieldmodule.createFieldcache() - node = findNodeWithName(markerNodes, markerName, "APEX") + node = findNodeWithName(markerNodes, markerName, "apex of heart") self.assertTrue(node.isValid()) cache.setNode(node) element, xi = markerLocation.evaluateMeshLocation(cache, 3) self.assertEqual(5, element.getIdentifier()) assertAlmostEqualList(self, xi, [ 0.0, 0.0, 1.0 ], 1.0E-10) + apexGroup = getAnnotationGroupForTerm(annotationGroups, get_heart_term("apex of heart")) + self.assertTrue(apexGroup.getNodesetGroup(nodes).containsNode(node)) if __name__ == "__main__": unittest.main()