diff --git a/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py b/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py index eb52614b..598c458e 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_heartventricles2.py @@ -4,14 +4,19 @@ from __future__ import division import math +import scaffoldmaker.utils.vector as vector from scaffoldmaker.meshtypes.meshtype_3d_sphereshell1 import MeshType_3d_sphereshell1 from scaffoldmaker.utils.eft_utils import * -from scaffoldmaker.utils.zinc_utils import * +from scaffoldmaker.utils.geometry import * +from scaffoldmaker.utils.interpolation import * from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from scaffoldmaker.utils.zinc_utils import * from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node + class MeshType_3d_heartventricles2: ''' classdocs @@ -23,106 +28,125 @@ def getName(): @staticmethod def getDefaultOptions(): return { - 'Number of elements up' : 4, - 'Number of elements around' : 10, - 'Number of elements across septum' : 4, - 'Number of elements below septum' : 2, - 'Number of elements through LV wall' : 1, - 'LV wall thickness' : 0.15, - 'LV wall thickness ratio apex' : 0.5, - 'LV wall thickness ratio base': 0.5, - 'RV free wall thickness' : 0.05, - 'RV width' : 0.2, - 'Length ratio' : 2.0, - 'Element length ratio equator/apex' : 1.0, - 'Use cross derivatives' : False + 'Number of elements around LV free wall' : 5, + 'Number of elements around septum' : 6, + 'Number of elements up apex' : 1, + 'Number of elements up septum' : 3, + 'Total height' : 1.0, + 'LV outer radius' : 0.5, + 'LV free wall thickness' : 0.12, + 'LV apex thickness' : 0.06, + 'RV height' : 0.8, + 'RV arc around degrees' : 180.0, + 'RV free wall thickness' : 0.04, + 'RV width' : 0.4, + 'RV extra cross radius base' : 0.05, + 'Septum thickness' : 0.1, + 'Septum base radial displacement' : 0.15, + 'Use cross derivatives' : False, + 'Refine' : False, + 'Refine number of elements surface' : 1, + 'Refine number of elements through LV wall' : 1, + 'Refine number of elements through RV wall' : 1 } @staticmethod def getOrderedOptionNames(): return [ - 'Number of elements up', - 'Number of elements around', - 'Number of elements across septum', - 'Number of elements below septum', - 'Number of elements through LV wall', - 'LV wall thickness', - 'LV wall thickness ratio apex', - 'LV wall thickness ratio base', + 'Number of elements around LV free wall', + 'Number of elements around septum', + 'Number of elements up apex', + 'Number of elements up septum', + 'Total height', + 'LV outer radius', + 'LV free wall thickness', + 'LV apex thickness', + 'RV height', + 'RV arc around degrees', 'RV free wall thickness', 'RV width', - 'Length ratio', - 'Element length ratio equator/apex' + 'RV extra cross radius base', + 'Septum thickness', + 'Septum base radial displacement', + 'Refine', + 'Refine number of elements surface', + 'Refine number of elements through LV wall', + 'Refine number of elements through RV wall' ] @staticmethod def checkOptions(options): - if options['Number of elements up'] < 4: - options['Number of elements up'] = 4 - if options['Number of elements around'] < 4: - options['Number of elements around'] = 4 - if options['Number of elements across septum'] < 2: - options['Number of elements across septum'] = 2 - elif options['Number of elements across septum'] > (options['Number of elements around'] - 2): - options['Number of elements across septum'] = options['Number of elements around'] - 2 - if options['Number of elements below septum'] < 2: - options['Number of elements below septum'] = 2 - elif options['Number of elements below septum'] > (options['Number of elements up'] - 1): - options['Number of elements below septum'] = options['Number of elements up'] - 1 - if (options['Number of elements through LV wall'] < 1) : - options['Number of elements through LV wall'] = 1 - if options['LV wall thickness'] < 0.0: - options['LV wall thickness'] = 0.0 - elif options['LV wall thickness'] > 0.5: - options['LV wall thickness'] = 0.5 - if options['LV wall thickness ratio apex'] < 0.0: - options['LV wall thickness ratio apex'] = 0.0 - if options['LV wall thickness ratio base'] < 0.0: - options['LV wall thickness ratio base'] = 0.0 - if options['RV free wall thickness'] < 0.0: - options['RV free wall thickness'] = 0.0 - if options['RV width'] < 0.0: - options['RV width'] = 0.0 - if options['Length ratio'] < 1.0E-6: - options['Length ratio'] = 1.0E-6 - if options['Element length ratio equator/apex'] < 1.0E-6: - options['Element length ratio equator/apex'] = 1.0E-6 + for key in [ + 'Refine number of elements surface', + 'Refine number of elements through LV wall', + 'Refine number of elements through RV wall', + 'Number of elements up apex']: + if options[key] < 1: + options[key] = 1 + for key in [ + 'Number of elements around LV free wall', + 'Number of elements up septum']: + if options[key] < 2: + options[key] = 2 + for key in [ + 'Number of elements around septum']: + if options[key] < 3: + options[key] = 3 + for key in [ + 'Total height', + 'LV outer radius', + 'LV free wall thickness', + 'LV apex thickness', + 'RV height', + 'RV free wall thickness', + 'RV width', + 'RV extra cross radius base', + 'Septum thickness', + 'Septum base radial displacement']: + if options[key] < 0.0: + options[key] = 0.0 + if options['RV height'] > 0.99*options['Total height']: + options['RV height'] = 0.99*options['Total height'] + if options['RV arc around degrees'] < 45.0: + options['RV arc around degrees'] = 45.0 + elif options['RV arc around degrees'] > 270.0: + options['RV arc around degrees'] = 270.0 @staticmethod - def generateMesh(region, options): + def generateBaseMesh(region, options): """ + Generate the base tricubic Hermite mesh. See also generateMesh(). :param region: Zinc region to define model in. Must be empty. :param options: Dict containing options. See getDefaultOptions(). :return: None """ - elementsCountUp = options['Number of elements up'] - elementsCountAround = options['Number of elements around'] - elementsCountAcrossSeptum = options['Number of elements across septum'] - elementsCountBelowSeptum = options['Number of elements below septum'] - elementsCountThroughLVWall = options['Number of elements through LV wall'] - LVWallThickness = options['LV wall thickness'] - LVWallThicknessRatioApex = options['LV wall thickness ratio apex'] - LVWallThicknessRatioBase = options['LV wall thickness ratio base'] - RVWidthTop = options['RV width'] - RVFreeWallThickness = options['RV free wall thickness'] + elementsCountAroundLVFreeWall = options['Number of elements around LV free wall'] + elementsCountAroundSeptum = options['Number of elements around septum'] + elementsCountAroundLV = elementsCountAroundLVFreeWall + elementsCountAroundSeptum + elementsCountUpApex = options['Number of elements up apex'] + elementsCountUpSeptum = options['Number of elements up septum'] + elementsCountUp = elementsCountUpApex + elementsCountUpSeptum + totalHeight = options['Total height'] + lvOuterRadius = options['LV outer radius'] + lvFreeWallThickness = options['LV free wall thickness'] + lvApexThickness = options['LV apex thickness'] + rvHeight = options['RV height'] + rvArcAroundRadians = math.radians(options['RV arc around degrees']) + rvFreeWallThickness = options['RV free wall thickness'] + rvWidth = options['RV width'] + rvExtraCrossRadiusBase = options['RV extra cross radius base'] + septumThickness = options['Septum thickness'] + septumBaseRadialDisplacement = options['Septum base radial displacement'] useCrossDerivatives = options['Use cross derivatives'] - # generate a half sphere shell which will be edited - sphereShellOptions = MeshType_3d_sphereshell1.getDefaultOptions() - sphereShellOptions['Number of elements up'] = elementsCountUp*2 - sphereShellOptions['Number of elements around'] = elementsCountAround - sphereShellOptions['Number of elements through wall'] = elementsCountThroughLVWall - sphereShellOptions['Exclude top rows'] = elementsCountUp - sphereShellOptions['Wall thickness'] = LVWallThickness - sphereShellOptions['Wall thickness ratio apex'] = LVWallThicknessRatioApex - sphereShellOptions['Length ratio'] = options['Length ratio'] - sphereShellOptions['Element length ratio equator/apex'] = options['Element length ratio equator/apex'] - MeshType_3d_sphereshell1.generateMesh(region, sphereShellOptions) - fm = region.getFieldmodule() fm.beginChange() - # find the coordinates field created for the sphere shell - coordinates = fm.findFieldByName('coordinates').castFiniteElement() + coordinates = getOrCreateCoordinateField(fm) + cache = fm.createFieldcache() + + ################# + # Create nodes + ################# nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) nodetemplate = nodes.createNodetemplate() @@ -131,539 +155,872 @@ def generateMesh(region, options): 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) + nodetemplateApex = nodetemplate - mesh = fm.findMeshByDimension(3) + nodeIdentifier = 1 - cache = fm.createFieldcache() + # LV nodes - if LVWallThicknessRatioBase != 1.0: - # make LV walls thinner at base - # get inside node at middle of RV - now = 1 + elementsCountUp*elementsCountAround - midRVnid = now - elementsCountAround + 2 + (elementsCountAcrossSeptum // 2) - #print('midRVnid', midRVnid) - midRVnode = nodes.findNodeByIdentifier(midRVnid) - cp_coordinates = fm.createFieldCoordinateTransformation(coordinates) - cp_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_CYLINDRICAL_POLAR) - radius = fm.createFieldComponent(cp_coordinates, 1) - theta = fm.createFieldComponent(cp_coordinates, 2) - z = fm.createFieldComponent(cp_coordinates, 3) - cache.setNode(midRVnode) - #result, cp = cp_coordinates.evaluateReal(cache, 3) - #coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) - result, innerRadius = radius.evaluateReal(cache, 1) - #print('innerRadius', innerRadius) - ir = fm.createFieldConstant([innerRadius*0.9999]) - radius_gt_ir = fm.createFieldGreaterThan(radius, ir) - radius_minus_ir = fm.createFieldSubtract(radius, ir) - thickness_scale = fm.createFieldConstant([LVWallThicknessRatioBase]) - delta_radius = fm.createFieldMultiply(radius_minus_ir, thickness_scale) - new_radius = fm.createFieldAdd(ir, delta_radius) - new_cp_coordinates = fm.createFieldConcatenate([new_radius, theta, z]) - new_cp_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_CYLINDRICAL_POLAR) - new_coordinates = fm.createFieldCoordinateTransformation(new_cp_coordinates) - new_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN) - - baseNodeGroupField = fm.createFieldNodeGroup(nodes) - just_under_zero = fm.createFieldConstant([-0.0001]) - isBase = fm.createFieldGreaterThan(z, just_under_zero) - baseNodesetGroup = baseNodeGroupField.getNodesetGroup() - baseNodesetGroup.addNodesConditional(isBase) - #print('baseNodesetGroup.getSize()', baseNodesetGroup.getSize()) - fieldassignment = coordinates.createFieldassignment(new_coordinates) - result = fieldassignment.setNodeset(baseNodesetGroup) - #print('fieldassignment.setNodeset', result) - fieldassignment.assign() + rvOuterWidthBase = rvWidth - septumBaseRadialDisplacement + septumThickness - lvFreeWallThickness + rvFreeWallThickness - tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) - tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + # get element spacing up + # approximate correction to move up to next inner node, since bottom of RV is halfway up element + rvHeightCorrection = (elementsCountUpSeptum - 1)/(elementsCountUpSeptum - 0.5) + rvThetaUp = math.acos(rvHeight*rvHeightCorrection/totalHeight) + rvOuterHeight = 0.5*(rvHeight + totalHeight) + rvCrossHeight = rvHeight + radialDisplacementStartRadiansUp = math.acos(rvOuterHeight/totalHeight) + lvArcAroundRadians = 2.0*math.pi - rvArcAroundRadians + + norl = elementsCountAroundLV + nowl = elementsCountUp*elementsCountAroundLV + 1 + for n3 in range(2): + a = totalHeight + b = lvOuterRadius + bVS = lvOuterRadius - lvFreeWallThickness + septumThickness + if n3 == 0: + a -= lvApexThickness + b -= lvFreeWallThickness + bVS = b + + lengthUp = 0.25*getApproximateEllipsePerimeter(a, b) + lengthUpApex = getEllipseArcLength(a, b, 0.0, rvThetaUp) + lengthUpSeptum = lengthUp - lengthUpApex + elementSizeUpSeptum = lengthUpSeptum/(elementsCountUpSeptum - 1) + #if n3 == 1: + # elementSizeUpApex = (lengthUpApex - (2.0/3.0)*elementSizeUpSeptum)/(elementsCountUpApex + (2.0/3.0)) + # elementSizeUpSeptumTransition = (2.0/3.0)*(elementSizeUpApex + elementSizeUpSeptum) + #else: + elementSizeUpApex = (lengthUpApex - 0.5*elementSizeUpSeptum)/(elementsCountUpApex + 0.5) + elementSizeUpSeptumTransition = 0.5*(elementSizeUpApex + elementSizeUpSeptum) + # have reduced derivative around RV junction to fit extra row of elements: + d2FactorRvTransition = 0.5*elementSizeUpSeptumTransition/elementSizeUpApex # 2.0/3.0 - eft = tricubichermite.createEftBasic() - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate.defineField(coordinates, -1, eft) - - crossAngle = math.pi/8 - sinCrossAngle = math.sin(crossAngle) - cosCrossAngle = math.cos(crossAngle) - edgeRotateCrossAngle = math.pi/16 - sinEdgeRotateCrossAngle = math.sin(edgeRotateCrossAngle) - cosEdgeRotateCrossAngle = math.cos(edgeRotateCrossAngle) - - # create RV nodes and modify adjoining LV nodes - elementsCountUpRV = elementsCountUp - elementsCountBelowSeptum - now = 1 + elementsCountUp*elementsCountAround - nodeIdentifier = 1 + 2*elementsCountThroughLVWall*now - rv_nidsInner = [] - rv_nidsOuter = [] - baseExtraRVWidth = LVWallThickness*(1.0 - LVWallThicknessRatioBase) - for n3 in range(1, -1, -1): # only 1 element through RV free wall - - onInside = n3 == 0 - - for n2 in range(elementsCountUpRV + 1): - - #if n2 > (elementsCountUpRV - elementsCountUpBase): - # theta = math.pi*0.5 - #else: - # theta = (n2 / (elementsCountUpRV - elementsCountUpBase))*math.pi*0.5 - #theta = (n2 / elementsCountUpRV)*math.pi*0.5 - RVWidth = RVWidthTop #*math.sin(theta) - #RVWidth = RVWidthTop - - edgeOffsetInner = RVWidth*0.125 - if n2 == elementsCountUpRV: - RVWidth += baseExtraRVWidth - edgeOffsetInner = edgeOffsetInner + baseExtraRVWidth*0.5 - - onBottomEdge = (n2 == 0) - - for n1 in range(elementsCountAcrossSeptum + 1): - - onSideEdge1 = (n1 == 0) - onSideEdge2 = (n1 == elementsCountAcrossSeptum) - onSideEdge = onSideEdge1 or onSideEdge2 - - nid = 3 + elementsCountThroughLVWall*now + (elementsCountBelowSeptum + n2 - 1)*elementsCountAround + n1 - - baseNode = nodes.findNodeByIdentifier(nid) - cache.setNode(baseNode) - result, base_x = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) - result, base_dx_ds1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) - result, base_dx_ds2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) - result, base_dx_ds3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, 3) - #print('node ', nid, 'dx_ds3', result, base_dx_ds3) - mag = math.sqrt(base_dx_ds3[0]*base_dx_ds3[0] + base_dx_ds3[1]*base_dx_ds3[1] + base_dx_ds3[2]*base_dx_ds3[2]) - unitOutward = [ base_dx_ds3[0]/mag, base_dx_ds3[1]/mag, base_dx_ds3[2]/mag ] - baseRadiusAround = unitOutward[0]*base_x[0] + unitOutward[1]*base_x[1] + unitOutward[2]*base_x[2] - rotateEdge = False - if onInside: - if onBottomEdge and onSideEdge: - offset = 0.0 - elif onBottomEdge or onSideEdge: - offset = edgeOffsetInner + #print('length total apex septum', lengthUp, lengthUpApex, lengthUpSeptum, ' size apex septum', elementSizeUpApex, elementSizeUpSeptum) + + # apex node, noting s1, s2 is x, -y to get out outward s3 + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -a ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ elementSizeUpApex, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, -elementSizeUpApex, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, -lvApexThickness ]) + nodeIdentifier = nodeIdentifier + 1 + + radiansUp = 0.0 + for n2 in range(elementsCountUp): + + if n2 < elementsCountUpApex: + arcLengthUp = elementSizeUpApex + dArcLengthUp = elementSizeUpApex + else: + if n2 == elementsCountUpApex: + arcLengthUp = elementSizeUpSeptumTransition + else: + arcLengthUp = elementSizeUpSeptum + dArcLengthUp = elementSizeUpSeptum + radiansUp = updateEllipseAngleByArcLength(a, b, radiansUp, arcLengthUp) + if n2 == (elementsCountUp - 1): + radiansUp = 0.5*math.pi + #print(n3, n2, ' -> ', math.degrees(radiansUp), 'dArcLengthUp', dArcLengthUp) + cosRadiansUp = math.cos(radiansUp) + sinRadiansUp = math.sin(radiansUp) + z = -a*cosRadiansUp + radius = b*sinRadiansUp + bSeptum = b + if (n3 == 1) and (n2 >= elementsCountUpApex): + bSeptum -= (lvFreeWallThickness - septumThickness) + radiusSeptum = bSeptum*sinRadiansUp + + elementSizeCrossSeptum = 0.5*(radius*lvArcAroundRadians/elementsCountAroundLVFreeWall + septumThickness*sinRadiansUp) + elementSizeAroundLVFreeWall = (radius*lvArcAroundRadians - elementSizeCrossSeptum)/(elementsCountAroundLVFreeWall - 1) + elementSizeAroundLVFreeWallTransition = 0.5*(elementSizeAroundLVFreeWall + elementSizeCrossSeptum) + + # get radial displacement of centre of septum, a function of radiansUp + xiUp = max(0.0, (radiansUp - radialDisplacementStartRadiansUp)/(0.5*math.pi - radialDisplacementStartRadiansUp)) + radialDisplacement = interpolateCubicHermite([0.0], [0.0], [septumBaseRadialDisplacement], [0.0], xiUp)[0] + + radiansAround = -0.5*rvArcAroundRadians + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + if (n3 == 1) and (n2 < elementsCountUpApex): + # Follow RV on outside below septum + dEndMag = elementSizeCrossSeptum + xiUpZ = 1.0 + z/rvOuterHeight + xiUpCross = 1.0 + z/rvCrossHeight + rvAddWidthRadius, rvAddCrossRadius = getRVOuterSize(xiUpZ, xiUpCross, rvOuterWidthBase, rvExtraCrossRadiusBase) + sx, sd1 = getRvOuterPoints(rvArcAroundRadians, radius, rvAddWidthRadius, rvAddCrossRadius, elementsCountAroundSeptum, dEndMag, z, xiUpCross) + else: + sx, sd1 = getSeptumPoints(rvArcAroundRadians, radiusSeptum, radialDisplacement, elementsCountAroundLVFreeWall, elementsCountAroundSeptum, z, n3) + + # do finite difference to get d/dxi2 + dxi = 1.0E-3 + ds_dxi = dArcLengthUp + if (n3 == 1) and (n2 < elementsCountUpApex): + dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] + else: + dzr_dRadiansUp = [ a*sinRadiansUp, bSeptum*cosRadiansUp ] + ds_dRadiansUp = vector.magnitude(dzr_dRadiansUp) + dzr_ds = vector.normalise(dzr_dRadiansUp) + dz = dxi*ds_dxi*dzr_ds[0] + dr = dxi*ds_dxi*dzr_ds[1] + if (n3 == 1) and (n2 < elementsCountUpApex): + # Follow RV on outside below septum + radiansUpPlus = radiansUp + math.sqrt(dz*dz + dr*dr)/ds_dRadiansUp + dEndMag = 0.5*((radius + dr)*lvArcAroundRadians/elementsCountAroundLVFreeWall + septumThickness*math.sin(radiansUpPlus)) + xiUpZPlus = 1.0 + (z + dz)/rvOuterHeight + xiUpCross = 1.0 + (z + dz)/rvCrossHeight + rvAddWidthRadius, rvAddCrossRadius = getRVOuterSize(xiUpZPlus, xiUpCross, rvOuterWidthBase, rvExtraCrossRadiusBase) + tx, td1 = getRvOuterPoints(rvArcAroundRadians, radius + dr, rvAddWidthRadius, rvAddCrossRadius, elementsCountAroundSeptum, dEndMag, z + dz, xiUpCross) + else: + dxiUp_dxi = ds_dxi/(ds_dRadiansUp*(0.5*math.pi - radialDisplacementStartRadiansUp)) + dRadialDisplacement = dxi*dxiUp_dxi*interpolateCubicHermiteDerivative([0.0], [0.0], [septumBaseRadialDisplacement], [0.0], xiUp)[0] + tx, td1 = getSeptumPoints(rvArcAroundRadians, radiusSeptum + dr, radialDisplacement + dRadialDisplacement, elementsCountAroundLVFreeWall, elementsCountAroundSeptum, z + dz, n3) + # true values for LV freewall + dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] + dzr_ds = vector.normalise(dzr_dRadiansUp) + + for n1 in range(elementsCountAroundLV): + + if (n1 > 0) and (n1 < elementsCountAroundSeptum): + n1mid = (elementsCountAroundSeptum + 1)//2 + if n1 < n1mid: + nr = elementsCountAroundSeptum//2 - n1 + x = [ sx[nr][0], -sx[nr][1], sx[nr][2] ] + dx_ds1 = [ -sd1[nr][0], sd1[nr][1], sd1[nr][2] ] + xt = [ tx[nr][0], -tx[nr][1], tx[nr][2] ] else: - offset = RVWidth + nr = n1 - n1mid + x = sx[nr] + dx_ds1 = sd1[nr] + xt = tx[nr] + dx_ds2 = [ (xt[c] - x[c])/dxi for c in range(3) ] else: - if onBottomEdge and onSideEdge: - offset = RVFreeWallThickness - elif onBottomEdge or onSideEdge: - offset = edgeOffsetInner - rotateEdge = True + if n1 == 0: + radiansAround = -0.5*rvArcAroundRadians + dMagAround = elementSizeCrossSeptum + elif n1 == elementsCountAroundSeptum: + radiansAround = 0.5*rvArcAroundRadians + dMagAround = elementSizeCrossSeptum else: - offset = RVWidth + RVFreeWallThickness - x = [ - base_x[0] + unitOutward[0]*offset, - base_x[1] + unitOutward[1]*offset, - base_x[2] + unitOutward[2]*offset, - ] - scale1 = (baseRadiusAround + offset)/baseRadiusAround - dx_ds1 = [ - base_dx_ds1[0]*scale1, - base_dx_ds1[1]*scale1, - base_dx_ds1[2]*scale1 - ] - # GRC not sure if appropriate to use scale1 here: - dx_ds2 = [ - base_dx_ds2[0]*scale1, - base_dx_ds2[1]*scale1, - base_dx_ds2[2]*scale1 - ] - scale3 = RVFreeWallThickness - if onInside and (onBottomEdge or onSideEdge) and not (onBottomEdge and onSideEdge): - scale3 += 2.0*edgeOffsetInner - dx_ds3 = [ - unitOutward[0]*scale3, - unitOutward[1]*scale3, - unitOutward[2]*scale3 - ] + radiansAround = 0.5*rvArcAroundRadians + (elementSizeAroundLVFreeWallTransition + (n1 - elementsCountAroundSeptum - 1)*elementSizeAroundLVFreeWall)/radius + dMagAround = elementSizeAroundLVFreeWall + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + x = [ radius*cosRadiansAround, radius*sinRadiansAround, z ] + dx_ds1 = [ dMagAround*-sinRadiansAround, dMagAround*cosRadiansAround, 0.0 ] + dx_ds2 = [ cosRadiansAround*ds_dxi*dzr_ds[1], sinRadiansAround*ds_dxi*dzr_ds[1], ds_dxi*dzr_ds[0] ] - if rotateEdge: - if onSideEdge1: - rotatedOutward = [ - cosEdgeRotateCrossAngle*base_dx_ds3[0] - sinEdgeRotateCrossAngle*base_dx_ds1[0], - cosEdgeRotateCrossAngle*base_dx_ds3[1] - sinEdgeRotateCrossAngle*base_dx_ds1[1], - cosEdgeRotateCrossAngle*base_dx_ds3[2] - sinEdgeRotateCrossAngle*base_dx_ds1[2] - ] - elif onSideEdge2: - rotatedOutward = [ - cosEdgeRotateCrossAngle*base_dx_ds3[0] + sinEdgeRotateCrossAngle*base_dx_ds1[0], - cosEdgeRotateCrossAngle*base_dx_ds3[1] + sinEdgeRotateCrossAngle*base_dx_ds1[1], - cosEdgeRotateCrossAngle*base_dx_ds3[2] + sinEdgeRotateCrossAngle*base_dx_ds1[2] - ] - else: # onBottomEdge: - rotatedOutward = [ - cosEdgeRotateCrossAngle*base_dx_ds3[0] - sinEdgeRotateCrossAngle*base_dx_ds2[0], - cosEdgeRotateCrossAngle*base_dx_ds3[1] - sinEdgeRotateCrossAngle*base_dx_ds2[1], - cosEdgeRotateCrossAngle*base_dx_ds3[2] - sinEdgeRotateCrossAngle*base_dx_ds2[2] - ] - mag = math.sqrt(rotatedOutward[0]*rotatedOutward[0] + rotatedOutward[1]*rotatedOutward[1] + rotatedOutward[2]*rotatedOutward[2]) - unitRotatedOutward = [ rotatedOutward[0]/mag, rotatedOutward[1]/mag, rotatedOutward[2]/mag ] - scale3r = RVFreeWallThickness + edgeOffsetInner - x[0] += unitRotatedOutward[0]*scale3r - x[1] += unitRotatedOutward[1]*scale3r - x[2] += unitRotatedOutward[2]*scale3r - dx_ds3 = [ - unitRotatedOutward[0]*scale3r, - unitRotatedOutward[1]*scale3r, - unitRotatedOutward[2]*scale3r - ] - if onBottomEdge: - rscale2 = math.sqrt(dx_ds2[0]*dx_ds2[0] + dx_ds2[1]*dx_ds2[1] + dx_ds2[2]*dx_ds2[2]) - tang = [ - unitRotatedOutward[1]*dx_ds1[2] - unitRotatedOutward[2]*dx_ds1[1], - unitRotatedOutward[2]*dx_ds1[0] - unitRotatedOutward[0]*dx_ds1[2], - unitRotatedOutward[0]*dx_ds1[1] - unitRotatedOutward[1]*dx_ds1[0] - ] - rscale2 /= math.sqrt(tang[0]*tang[0] + tang[1]*tang[1] + tang[2]*tang[2]) - dx_ds2 = [ - tang[0]*rscale2, - tang[1]*rscale2, - tang[2]*rscale2 - ] - else: # onSideEdge - rscale1 = math.sqrt(dx_ds1[0]*dx_ds1[0] + dx_ds1[1]*dx_ds1[1] + dx_ds1[2]*dx_ds1[2]) - tang = [ - dx_ds2[1]*unitRotatedOutward[2] - dx_ds2[2]*unitRotatedOutward[1], - dx_ds2[2]*unitRotatedOutward[0] - dx_ds2[0]*unitRotatedOutward[2], - dx_ds2[0]*unitRotatedOutward[1] - dx_ds2[1]*unitRotatedOutward[0] - ] - rscale1 /= math.sqrt(tang[0]*tang[0] + tang[1]*tang[1] + tang[2]*tang[2]) - dx_ds1 = [ - tang[0]*rscale1, - tang[1]*rscale1, - tang[2]*rscale1 - ] - - if onInside and (onBottomEdge or onSideEdge): - node = baseNode - else: - node = nodes.createNode(nodeIdentifier, nodetemplate) - nodeIdentifier += 1 - cache.setNode(node) - if onInside: - rv_nidsInner.append(node.getIdentifier()) - else: - rv_nidsOuter.append(node.getIdentifier()) + if (n3 == 1) and (n2 == (elementsCountUpApex - 1)) and (n1 >= 0) and (n1 <= elementsCountAroundSeptum): + dx_ds2 = [ d*d2FactorRvTransition for d in dx_ds2 ] + + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) - - # create RV elements and modify adjoining LV element fields - elementIdentifier = elementsCountThroughLVWall*elementsCountUp*elementsCountAround + 1 - scalefactors5 = [ -1.0, sinCrossAngle, cosCrossAngle, sinCrossAngle, cosCrossAngle ] - scalefactors9 = [ -1.0, 0.5, 0.25, 0.125, 0.75, sinCrossAngle, cosCrossAngle, sinCrossAngle, cosCrossAngle ] - eow = elementsCountUp*elementsCountAround - RVSeptumElementIdBase = eow*(elementsCountThroughLVWall - 1) + elementsCountBelowSeptum*elementsCountAround + 2 - - # Refine elements around RV using hanging nodes - for n2 in range(-1, elementsCountUpRV): - - for n1 in range(-1, elementsCountAcrossSeptum + 1): - - if (n2 >= 0) and (n1 >= 0) and (n1 < elementsCountAcrossSeptum): - continue - - existingElementIdentifier = RVSeptumElementIdBase + n2*elementsCountAround + n1 - element = mesh.findElementByIdentifier(existingElementIdentifier) - eftInner = element.getElementfieldtemplate(coordinates, -1) - eftOuter = element.getElementfieldtemplate(coordinates, -1) - nodeIdentifiersInner = getElementNodeIdentifiers(element, eftInner) - nodeIdentifiersOuter = nodeIdentifiersInner[:] - # general scale factors 1 -> 1, 102 -> 1/2, 104 -> 1/4, 108 -> 1/8, 304 -> 3/4 - sfLimiter = 0 - if (n2 == -1) and ((n1 == -1) or (n1 == elementsCountAcrossSeptum)): - setEftScaleFactorIds(eftInner, [1, 102, 104, 108, 304], []) - setEftScaleFactorIds(eftOuter, [1, 102, 104, 108, 304], []) - elif (n2 == 0) or (n1 == 0) or (n1 == (elementsCountAcrossSeptum - 1)): - setEftScaleFactorIds(eftInner, [1, 102, 104, 108, 304], [1, 2]) - setEftScaleFactorIds(eftOuter, [1, 102, 104, 108, 304], [3, 4]) - sfLimiter = -2 # offset to make indexes 8/9 -> 6/7 - else: - setEftScaleFactorIds(eftInner, [1, 102, 104, 108, 304], [1, 2, 1, 2]) - setEftScaleFactorIds(eftOuter, [1, 102, 104, 108, 304], [3, 4, 3, 4]) - scaleFactorsCount = eftInner.getNumberOfLocalScaleFactors() - scaleFactorsSlice = scalefactors9[:scaleFactorsCount] - - if (n1 == -1) or (n2 == -1): - tricubichermite.setEftMidsideXi3HangingNode(eftInner, 5, 1, 1, 5, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi3HangingNode(eftOuter, 1, 5, 1, 5, [1, 2, 3, 4, 5]) - if n1 == -1: - tricubichermite.setEftMidsideXi3HangingNode(eftInner, 7, 3, 3, 7, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi3HangingNode(eftOuter, 3, 7, 3, 7, [1, 2, 3, 4, 5]) - if (n1 == elementsCountAcrossSeptum) or (n2 == -1): - tricubichermite.setEftMidsideXi3HangingNode(eftInner, 6, 2, 2, 6, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi3HangingNode(eftOuter, 2, 6, 2, 6, [1, 2, 3, 4, 5]) - if n1 == elementsCountAcrossSeptum: - tricubichermite.setEftMidsideXi3HangingNode(eftInner, 8, 4, 4, 8, [1, 2, 3, 4, 5]) - tricubichermite.setEftMidsideXi3HangingNode(eftOuter, 4, 8, 4, 8, [1, 2, 3, 4, 5]) - - if n2 == -1: # bottom - if (n1 > 0) and (n1 < elementsCountAcrossSeptum): - mapEftFunction1Node2Terms(eftInner, 6*8 + 5, 7, Node.VALUE_LABEL_D_DS2, 1, [6], Node.VALUE_LABEL_D_DS3, 1, [7]) - mapEftFunction1Node2Terms(eftOuter, 2*8 + 5, 3, Node.VALUE_LABEL_D_DS2, 1, [1, 6], Node.VALUE_LABEL_D_DS3, 1, [7]) - if n1 >= 0: - nodeIdentifiersOuter[2] = nodeIdentifiersInner[6] - nodeIdentifiersOuter[6] = rv_nidsOuter[n1] - if (n1 >= 0) and (n1 < (elementsCountAcrossSeptum - 1)): - mapEftFunction1Node2Terms(eftInner, 7*8 + 5, 8, Node.VALUE_LABEL_D_DS2, 1, [8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - mapEftFunction1Node2Terms(eftOuter, 3*8 + 5, 4, Node.VALUE_LABEL_D_DS2, 1, [1, 8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - if n1 < elementsCountAcrossSeptum: - nodeIdentifiersOuter[3] = nodeIdentifiersInner[7] - nodeIdentifiersOuter[7] = rv_nidsOuter[n1 + 1] - - if n1 == -1: # side 1 - nodeIdentifiersOuter[3] = nodeIdentifiersInner[7] - nodeIdentifiersOuter[7] = rv_nidsOuter[(n2 + 1)*(elementsCountAcrossSeptum + 1)] - if n2 > 0: - mapEftFunction1Node2Terms(eftInner, 5*8 + 5, 6, Node.VALUE_LABEL_D_DS1, 1, [6], Node.VALUE_LABEL_D_DS3, 1, [7]) - mapEftFunction1Node2Terms(eftOuter, 1*8 + 5, 2, Node.VALUE_LABEL_D_DS1, 1, [1, 6], Node.VALUE_LABEL_D_DS3, 1, [7]) - if n2 >= 0: - mapEftFunction1Node2Terms(eftInner, 7*8 + 5, 8, Node.VALUE_LABEL_D_DS1, 1, [8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - mapEftFunction1Node2Terms(eftOuter, 3*8 + 5, 4, Node.VALUE_LABEL_D_DS1, 1, [1, 8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - nodeIdentifiersOuter[1] = nodeIdentifiersInner[5] - nodeIdentifiersOuter[5] = rv_nidsOuter[n2*(elementsCountAcrossSeptum + 1)] - - if n1 == elementsCountAcrossSeptum: # side 2 - nodeIdentifiersOuter[2] = nodeIdentifiersInner[6] - nodeIdentifiersOuter[6] = rv_nidsOuter[(n2 + 1)*(elementsCountAcrossSeptum + 1) + n1] - if n2 > 0: - mapEftFunction1Node2Terms(eftInner, 4*8 + 5, 5, Node.VALUE_LABEL_D_DS1, 1, [1, 6], Node.VALUE_LABEL_D_DS3, 1, [7]) - mapEftFunction1Node2Terms(eftOuter, 0*8 + 5, 1, Node.VALUE_LABEL_D_DS1, 1, [6], Node.VALUE_LABEL_D_DS3, 1, [7]) - if n2 >= 0: - mapEftFunction1Node2Terms(eftInner, 6*8 + 5, 7, Node.VALUE_LABEL_D_DS1, 1, [1, 8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - mapEftFunction1Node2Terms(eftOuter, 2*8 + 5, 3, Node.VALUE_LABEL_D_DS1, 1, [8 + sfLimiter], Node.VALUE_LABEL_D_DS3, 1, [9 + sfLimiter]) - nodeIdentifiersOuter[0] = nodeIdentifiersInner[4] - nodeIdentifiersOuter[4] = rv_nidsOuter[n2*(elementsCountAcrossSeptum + 1) + n1] - - elementtemplateInner = mesh.createElementtemplate() - elementtemplateInner.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateInner.defineField(coordinates, -1, eftInner) - result1 = element.merge(elementtemplateInner) - result2 = element.setNodesByIdentifier(eftInner, nodeIdentifiersInner) - result3 = element.setScaleFactors(eftInner, scaleFactorsSlice) - #print('merge inner hanging element', existingElementIdentifier, result1, result2, result3, nodeIdentifiersInner, scaleFactorsSlice) - - elementtemplateOuter = mesh.createElementtemplate() - elementtemplateOuter.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateOuter.defineField(coordinates, -1, eftOuter) - element = mesh.createElement(elementIdentifier, elementtemplateOuter) - result2 = element.setNodesByIdentifier(eftOuter, nodeIdentifiersOuter) - result3 = element.setScaleFactors(eftOuter, scaleFactorsSlice) - #print('create outer hanging element', elementIdentifier, result2, result3, nodeIdentifiersOuter, scaleFactorsSlice) - elementIdentifier += 1 - - # Tweak RV septal wall - - # RV septal wall inner bottom elements - eftRVSeptalWallInnerBottom = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVSeptalWallInnerBottom, [1], [3, 4, 3, 4]) - for s in range(2): - n = 4 + s - ln = n + 1 - # d/dxi2 = -d/ds3, d/dxi3 = sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVSeptalWallInnerBottom, n*8 + 3, ln, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftRVSeptalWallInnerBottom, n*8 + 5, ln, Node.VALUE_LABEL_D_DS2, 1, [s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVSeptalWallInnerBottom = mesh.createElementtemplate() - elementtemplateRVSeptalWallInnerBottom.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVSeptalWallInnerBottom.defineField(coordinates, -1, eftRVSeptalWallInnerBottom) - #print('eftRVSeptalWallInnerBottom', result) - - # RV septal wall inner side 1 elements - eftRVSeptalWallInnerSide1 = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVSeptalWallInnerSide1, [1], [3, 4, 3, 4]) - for s in range(2): - n = 4 + s*2 - ln = n + 1 - # d/dxi1 = -d/ds3, d/dxi3 = sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVSeptalWallInnerSide1, n*8 + 2, ln, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftRVSeptalWallInnerSide1, n*8 + 5, ln, Node.VALUE_LABEL_D_DS1, 1, [s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVSeptalWallInnerSide1 = mesh.createElementtemplate() - elementtemplateRVSeptalWallInnerSide1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVSeptalWallInnerSide1.defineField(coordinates, -1, eftRVSeptalWallInnerSide1) - #print('eftRVSeptalWallInnerSide1', result) - - # RV septal wall inner side 2 elements - eftRVSeptalWallInnerSide2 = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVSeptalWallInnerSide2, [1], [3, 4, 3, 4]) - for s in range(2): - n = 5 + s*2 - ln = n + 1 - # d/dxi1 = d/ds3, d/dxi3 = -sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVSeptalWallInnerSide2, n*8 + 2, ln, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftRVSeptalWallInnerSide2, n*8 + 5, ln, Node.VALUE_LABEL_D_DS1, 1, [1, s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVSeptalWallInnerSide2 = mesh.createElementtemplate() - elementtemplateRVSeptalWallInnerSide2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVSeptalWallInnerSide2.defineField(coordinates, -1, eftRVSeptalWallInnerSide2) - #print('eftRVSeptalWallInnerSide2', result) + if n3 == 1: + # get derivative3 from difference in coordinates + innerNode = nodes.findNodeByIdentifier(nodeIdentifier - nowl) + cache.setNode(innerNode) + result, xInner = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) + dx_ds3 = [ (x[c] - xInner[c]) for c in range(3) ] + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) + nodeIdentifier = nodeIdentifier + 1 + + # RV nodes + + nidr = nodeIdentifier + elementsCountUpRV = elementsCountUpSeptum + 1 + elementsCountAroundRV = elementsCountAroundSeptum + 2 + + rxOuter = [] + rd1Outer = [] + rd2Outer = [] + rd3Outer = [] + elementSizeUpRvTransition1 = d2FactorRvTransition*elementSizeUpApex + elementSizeUpRv = (lengthUp - elementsCountUpApex*elementSizeUpApex - 1.5*elementSizeUpRvTransition1)/(elementsCountUpRV - 1.5) + elementSizeUpRvTransition2 = 0.5*(elementSizeUpRvTransition1 + elementSizeUpRv) + distUp = elementSizeUpApex*elementsCountUpApex + elementSizeUpRvTransition1 + radiansUp = updateEllipseAngleByArcLength(a, b, 0.0, distUp) + + #print('RV element size apex', elementSizeUpApex, ', transition1', elementSizeUpRvTransition2, ', ', transition2', elementSizeUpRvTransition1, ', rv', elementSizeUpRv, \ + # ', total', elementsCountUpApex*elementsCountUpApex + elementSizeUpRvTransition1 + elementSizeUpRvTransition2 + (elementsCountUpRV - 2)*elementSizeUpRv, ' vs ', lengthUp) for n2 in range(elementsCountUpRV): + if n2 == (elementsCountUpRV - 1): + radiansUp = 0.5*math.pi + cosRadiansUp = math.cos(radiansUp) + sinRadiansUp = math.sin(radiansUp) - for n1 in range(elementsCountAcrossSeptum): - existingElementIdentifier = RVSeptumElementIdBase + n2*elementsCountAround + n1 - - element = mesh.findElementByIdentifier(existingElementIdentifier) - eftTemp = element.getElementfieldtemplate(coordinates, -1) - nodeIdentifiers = getElementNodeIdentifiers(element, eftTemp) - #print('existing element', existingElementIdentifier, eftTemp.isValid(), nodeIdentifiers) - if n2 == 0: - if n1 == 0: - # RV free wall inner bottom side 1 elements: doubly curved but can't avoid a sharp corner at node 1 - eftTemp = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftTemp, [1], [3, 4, 3, 4]) - # node 6 d/dxi2 = -d/ds3, d/dxi3 = sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 5*8 + 3, 6, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftTemp, 5*8 + 5, 6, Node.VALUE_LABEL_D_DS2, 1, [2], Node.VALUE_LABEL_D_DS3, 1, [3]) - # node 7 d/dxi1 = -d/ds3, d/dxi3 = sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 6*8 + 2, 7, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftTemp, 6*8 + 5, 7, Node.VALUE_LABEL_D_DS1, 1, [1, 4], Node.VALUE_LABEL_D_DS3, 1, [5]) - #print('inner bottom side 1 eftTemp.isValid()',eftTemp.isValid()) - elementtemplateTemp = mesh.createElementtemplate() - elementtemplateTemp.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateTemp.defineField(coordinates, -1, eftTemp) - result = element.merge(elementtemplateTemp) - result = element.setNodesByIdentifier(eftTemp, nodeIdentifiers) - result = element.setScaleFactors(eftTemp, scalefactors5) - elif n1 == (elementsCountAcrossSeptum - 1): - # RV free wall inner bottom side 1 elements: doubly curved but can't avoid a sharp corner at node 1 - eftTemp = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftTemp, [1], [3, 4, 3, 4]) - # node 5 d/dxi2 = -d/ds3, d/dxi3 = sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 4*8 + 3, 5, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftTemp, 4*8 + 5, 5, Node.VALUE_LABEL_D_DS2, 1, [2], Node.VALUE_LABEL_D_DS3, 1, [3]) - # node 8 d/dxi1 = d/ds3, d/dxi3 = s2.d/ds1 + s3.d/ds3 - mapEftFunction1Node1Term(eftTemp, 7*8 + 2, 8, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftTemp, 7*8 + 5, 8, Node.VALUE_LABEL_D_DS1, 1, [1, 4], Node.VALUE_LABEL_D_DS3, 1, [5]) - #print('inner bottom side 2 eftTemp.isValid()',eftTemp.isValid()) - elementtemplateTemp = mesh.createElementtemplate() - elementtemplateTemp.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateTemp.defineField(coordinates, -1, eftTemp) - result = element.merge(elementtemplateTemp) - result = element.setNodesByIdentifier(eftTemp, nodeIdentifiers) - result = element.setScaleFactors(eftTemp, scalefactors5) - else: - result = element.merge(elementtemplateRVSeptalWallInnerBottom) - result = element.setNodesByIdentifier(eftRVSeptalWallInnerBottom, nodeIdentifiers) - result = element.setScaleFactors(eftRVSeptalWallInnerBottom, scalefactors5) + radiansAround = -0.5*rvArcAroundRadians + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + + lvOuterRadius = b*math.sin(radiansUp) + ax = [ lvOuterRadius*cosRadiansAround, lvOuterRadius*sinRadiansAround, -a*cosRadiansUp ] + z = ax[2] + elementSizeCrossSeptum = 0.5*(lvOuterRadius*lvArcAroundRadians/elementsCountAroundLVFreeWall + septumThickness*sinRadiansUp) + dEndMag = elementSizeCrossSeptum + ad2 = [ b*cosRadiansUp*cosRadiansAround, b*cosRadiansUp*sinRadiansAround, a*sinRadiansUp ] + scale = elementSizeUpRv/vector.magnitude(ad2) + if n2 == 0: + scale *= elementSizeUpRvTransition1/elementSizeUpRv + ad2 = [ d*scale for d in ad2 ] + + xiUpZ =1.0 + z/rvOuterHeight + xiUpCross = 1.0 + z/rvCrossHeight + + rvAddWidthRadius, rvAddCrossRadius = getRVOuterSize(xiUpZ, xiUpCross, rvOuterWidthBase, rvExtraCrossRadiusBase) + #print('z',z,'xiUpZ', xiUpZ, 'xiUpFast', xiUpFast, 'xiUpSlow', xiUpSlow) + #print('rvAddWidthRadius =', rvAddWidthRadius, ', rvAddCrossRadius=', rvAddCrossRadius) + rx, rd1 = getRvOuterPoints(rvArcAroundRadians, lvOuterRadius, rvAddWidthRadius, rvAddCrossRadius, elementsCountAroundRV, dEndMag, z, xiUpCross, True) + + if n2 < (elementsCountUpRV - 1): + # get points at node position plus increment of xi2, to calculate derivative 2 + dxi = 1.0E-3 + xPlus = [ ax[i] + ad2[i]*dxi for i in range(3) ] + lvOuterRadiusPlus = math.sqrt(xPlus[0]*xPlus[0] + xPlus[1]*xPlus[1]) + dArcLengthUp = elementSizeUpRv + ds_dxi = dArcLengthUp + dzr_dRadiansUp = [ a*sinRadiansUp, b*cosRadiansUp ] + ds_dRadiansUp = vector.magnitude(dzr_dRadiansUp) + dzr_ds = vector.normalise(dzr_dRadiansUp) + dz = dxi*ds_dxi*dzr_ds[0] + dr = dxi*ds_dxi*dzr_ds[1] + radiansUpPlus = radiansUp + math.sqrt(dz*dz + dr*dr)/ds_dRadiansUp + dEndMagPlus = 0.5*(lvOuterRadiusPlus*lvArcAroundRadians/elementsCountAroundLVFreeWall + septumThickness*math.sin(radiansUpPlus)) + xiUpZPlus = 1.0 + xPlus[2]/rvOuterHeight + xiUpCross = 1.0 + xPlus[2]/rvCrossHeight + rvAddWidthRadiusPlus, rvAddCrossRadiusPlus = getRVOuterSize(xiUpZPlus, xiUpCross, rvOuterWidthBase, rvExtraCrossRadiusBase) + px, pd1 = getRvOuterPoints(rvArcAroundRadians, lvOuterRadiusPlus, rvAddWidthRadiusPlus, rvAddCrossRadiusPlus, elementsCountAroundRV, dEndMagPlus, xPlus[2], xiUpCross, True) + + rd2 = [ [(px[n][c] - rx[n][c])/dxi for c in range(3)] for n in range(len(rx)) ] + else: + rd2 = [ad2]*len(rx) + + #if n2 == 0: + # for n in range(len(rd2)): + # rd2[n] = [ d*elementSizeUpRvTransition1/elementSizeUpRv for d in rd2[n] ] + + rd3 = [] + for n in range(len(rx)): + normal = vector.crossproduct3(rd1[n], rd2[n]) + mag = rvFreeWallThickness/vector.magnitude(normal) + rd3.append([ c*mag for c in normal ]) + + rxOuter.append(rx) + rd1Outer.append(rd1) + rd2Outer.append(rd2) + rd3Outer.append(rd3) + radiansUp = updateEllipseAngleByArcLength(a, b, radiansUp, elementSizeUpRv if (n2 > 0) else elementSizeUpRvTransition2) + + # get inner RV nodes from outer + + nidl = 2 + (elementsCountUpApex - 1)*norl + + rxInner = [] + rd1Inner = [] + rd2Inner = [] + rd3Inner = [] + # fix bottom row below + for n2 in range(1, elementsCountUpRV): + ix = [] + id1 = [] + id2 = [] + id3 = [] + rx = rxOuter[n2] + rd1 = rd1Outer[n2] + rd2 = rd2Outer[n2] + rd3 = rd3Outer[n2] + + for n1 in range(1, elementsCountAroundRV - 2): + ix.append([ (rx[n1][c] - rd3[n1][c]) for c in range(3) ] ) + unitRadial = vector.normalise(rd3[n1]) + + # calculate inner d1 from curvature around + curvature = 0.0 + count = 0 + if n1 > 0: + curvature -= getCubicHermiteCurvature(rx[n1 - 1], rd1[n1 - 1], rx[n1], rd1[n1], unitRadial, 1.0) + count += 1 + if n1 < (elementsCountAroundRV - 2): + curvature -= getCubicHermiteCurvature(rx[n1], rd1[n1], rx[n1 + 1], rd1[n1 + 1], unitRadial, 0.0) + count += 1 + curvature /= count + factor = 1.0 - curvature*rvFreeWallThickness + id1.append([ factor*c for c in rd1[n1] ]) + + # calculate inner d2 from curvature up + curvature = 0.0 + count = 0.0 + if n2 > 0: + curvature -= getCubicHermiteCurvature(rxOuter[n2 - 1][n1], rd2Outer[n2 - 1][n1], rx[n1], rd2[n1], unitRadial, 1.0) + count += 1 + if n2 < (elementsCountUpRV - 1): + curvature -= getCubicHermiteCurvature(rx[n1], rd2[n1], rxOuter[n2 + 1][n1], rd2Outer[n2 + 1][n1], unitRadial, 0.0) + count += 1 + curvature /= count + factor = 1.0 - curvature*rvFreeWallThickness + id2.append([ factor*c for c in rd2Outer[n2][n1] ]) + + id3.append(rd3[n1]) + + # RV inlet/outlet radius + + for n1 in [ 0, elementsCountAroundRV - 2 ]: + # interpolate node inside RV side edge + n1r = 0 if (n1 == 0) else (elementsCountAroundRV - 3) + nid = nidl + nowl + norl*n2 + (1 if (n1 == 0) else (elementsCountAroundSeptum - 1)) + node = nodes.findNodeByIdentifier(nid) + cache.setNode(node) + result, ax = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) + result, ad1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) + result, ad2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) + ad1 = [ -d for d in ad1 ] + bx = ix[n1r] + bd1 = id1[n1r] + bd2 = id2[n1r] + ad1 = [ 2.0*d for d in ad1 ] + bd1 = [ 2.0*d for d in bd1 ] + if n1 > 0: + ax, ad1, ad2, bx, bd1, bd2 = bx, bd1, bd2, ax, ad1, ad2 + x = list(interpolateCubicHermite(ax, ad1, bx, bd1, 0.5)) + dx_ds1 = interpolateCubicHermiteDerivative(ax, ad1, bx, bd1, 0.5) + dx_ds1 = [ 0.5*d for d in dx_ds1 ] + dx_ds2 = [ 0.5*(ad2[c] + bd2[c]) for c in range(3) ] + ox = rxOuter[n2][n1] + dx_ds3 = [ (ox[c] - x[c]) for c in range(3) ] + ix.insert(n1, x) + id1.insert(n1, dx_ds1) + id2.insert(n1, dx_ds2) + id3.insert(n1, dx_ds3) + rd3Outer[n2][n1] = dx_ds3 + + rxInner.append(ix) + rd1Inner.append(id1) + rd2Inner.append(id2) + rd3Inner.append(id3) + + # Bottom of RV + ix = [] + id1 = [] + id2 = [] + id3 = [] + for n1 in range(1, elementsCountAroundRV - 2): + # interpolate node inside RV bottom edge + nid = nidl + nowl + norl + n1 + #print(n2, 'n1', n1, 'n1r', n1r, 'LV nid', nid, ) + node = nodes.findNodeByIdentifier(nid) + cache.setNode(node) + result, ax = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3) + result, ad1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, 3) + result, ad2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3) + ad2 = [ -d for d in ad2 ] + bx = rxInner[0][n1] + bd1 = rd1Inner[0][n1] + bd2 = rd2Inner[0][n1] + ad2 = [ 1.5*d for d in ad2 ] + bd2 = [ 1.5*d for d in bd2 ] + x = list(interpolateCubicHermite(ax, ad2, bx, bd2, 0.5)) + dx_ds1 = [ 0.5*(ad1[c] + bd1[c]) for c in range(3) ] + dx_ds2 = interpolateCubicHermiteDerivative(ax, ad2, bx, bd2, 0.5) + dx_ds2 = [ (1.0/3.0)*d for d in dx_ds2 ] + ox = rxOuter[0][n1] + dx_ds3 = [ (ox[c] - x[c]) for c in range(3) ] + ix.append(x) + id1.append(dx_ds1) + id2.append(dx_ds2) + id3.append(dx_ds3) + rd3Outer[0][n1] = dx_ds3 + + # Bottom corners of RV + for n1 in [ 0, elementsCountAroundRV - 2 ]: + # interpolate node inside RV side edge + n1r = 0 if (n1 == 0) else (elementsCountAroundRV - 2) + ax = rxInner[0][n1r] + ad1 = rd2Inner[0][n1r] + n1r = 0 if (n1 == 0) else (elementsCountAroundRV - 3) + bx = ix[n1r] + bd1 = id1[n1r] + if n1 == 0: + ad1 = [ -d for d in ad1 ] + else: + ax, ad1, bx, bd1 = bx, bd1, ax, ad1 + ad1 = [ 1.5*d for d in ad1 ] + bd1 = [ 1.5*d for d in bd1 ] + x = list(interpolateCubicHermite(ax, ad1, bx, bd1, 0.5)) + dx_ds1 = interpolateCubicHermiteDerivative(ax, ad1, bx, bd1, 0.5) + dx_ds1 = [ (1.0/3.0)*d for d in dx_ds1 ] + if n1 == 0: + dx_ds2 = [ -d for d in dx_ds1 ] + else: + dx_ds2 = dx_ds1 + ox = rxOuter[0][n1] + dx_ds3 = [ (ox[c] - x[c]) for c in range(3) ] + ix.insert(n1, x) + id1.insert(n1, dx_ds1) + id2.insert(n1, dx_ds2) + id3.insert(n1, dx_ds3) + rd3Outer[0][n1] = dx_ds3 + + rxInner.insert(0, ix) + rd1Inner.insert(0, id1) + rd2Inner.insert(0, id2) + rd3Inner.insert(0, id3) + + zero = [ 0.0, 0.0, 0.0 ] + rv_nids = [[],[]] + for n3 in range(2): + for n2 in range(elementsCountUpRV): + nids = [] + rv_nids[n3].append(nids) + if n3 == 0: + rx = rxInner[n2] + rd1 = rd1Inner[n2] + rd2 = rd2Inner[n2] + rd3 = rd3Inner[n2] else: - if n1 == 0: - result = element.merge(elementtemplateRVSeptalWallInnerSide1) - result = element.setNodesByIdentifier(eftRVSeptalWallInnerSide1, nodeIdentifiers) - result = element.setScaleFactors(eftRVSeptalWallInnerSide1, scalefactors5) - elif n1 == (elementsCountAcrossSeptum - 1): - result = element.merge(elementtemplateRVSeptalWallInnerSide2) - result = element.setNodesByIdentifier(eftRVSeptalWallInnerSide2, nodeIdentifiers) - result = element.setScaleFactors(eftRVSeptalWallInnerSide2, scalefactors5) + rx = rxOuter[n2] + rd1 = rd1Outer[n2] + rd2 = rd2Outer[n2] + rd3 = rd3Outer[n2] + for n1 in range(elementsCountAroundRV - 1): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, rx[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, rd1[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, rd2[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, rd3[n1]) + nid = nodeIdentifier + nodeIdentifier += 1 + nids.append(nid) + + + ################# + # Create elements + ################# + + mesh = fm.findMeshByDimension(3) + tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + eft = tricubichermite.createEftNoCrossDerivatives() + + elementIdentifier = 1 + + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + norr = elementsCountAroundRV - 1 + nowr = norr*elementsCountUpRV + + # LV elements + + # calculate values used for apex scale factors + radiansPerElementAroundLVFreeWall = lvArcAroundRadians/elementsCountAroundLVFreeWall + if elementsCountAroundSeptum > 2: + radiansPerElementAroundSeptum = (rvArcAroundRadians - radiansPerElementAroundLVFreeWall)/(elementsCountAroundSeptum - 1) + radiansPerElementAroundSeptumTransition = 0.5*(radiansPerElementAroundLVFreeWall + radiansPerElementAroundSeptum) + else: + radiansPerElementAroundSeptum = radiansPerElementAroundSeptumTransition = rvArcAroundRadians/elementsCountAroundSeptum + radiansAround = -0.5*rvArcAroundRadians + + for e2 in range(elementsCountUp): + + for e1 in range(elementsCountAroundLV): + + eft1 = eft + scalefactors = None + if e2 == 0: + # create bottom apex elements, varying eft scale factor identifiers around apex + # scale factor identifiers follow convention of offsetting by 100 for each 'version' + va = e1 + vb = (e1 + 1)%elementsCountAroundLV + nids = [ 1 , 2 + va , 2 + vb, + 1 + nowl, 2 + va + nowl, 2 + vb + nowl ] + eft1 = tricubichermite.createEftShellApexBottom(va*100, vb*100) + # calculate general linear map coefficients + if e1 < elementsCountAroundSeptum: + deltaRadians = dRadians1 = dRadians2 = radiansPerElementAroundSeptum + if e1 == 0: + deltaRadians = radiansPerElementAroundSeptumTransition + dRadians1 = radiansPerElementAroundLVFreeWall + elif e1 == (elementsCountAroundSeptum - 1): + deltaRadians = radiansPerElementAroundSeptumTransition + dRadians2 = radiansPerElementAroundLVFreeWall else: - pass - #print('RV septal wall element merge', existingElementIdentifier, result) - - # RV free wall - - # RV free wall inner bottom elements - eftRVFreeWallInnerBottom = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVFreeWallInnerBottom, [1], [3, 4, 3, 4]) - for s in range(2): - n = s - ln = n + 1 - # d/dxi2 = d/ds3, d/dxi3 = -sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVFreeWallInnerBottom, n*8 + 3, ln, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftRVFreeWallInnerBottom, n*8 + 5, ln, Node.VALUE_LABEL_D_DS2, 1, [1, s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVFreeWallInnerBottom = mesh.createElementtemplate() - elementtemplateRVFreeWallInnerBottom.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVFreeWallInnerBottom.defineField(coordinates, -1, eftRVFreeWallInnerBottom) - #print('eftRVFreeWallInnerBottom', result) - - # RV free wall inner side 1 elements - eftRVFreeWallInnerSide1 = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVFreeWallInnerSide1, [1], [3, 4, 3, 4]) - for s in range(2): - n = s*2 - ln = n + 1 - # d/dxi1 = d/ds3, d/dxi3 = -sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVFreeWallInnerSide1, n*8 + 2, ln, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftRVFreeWallInnerSide1, n*8 + 5, ln, Node.VALUE_LABEL_D_DS1, 1, [1, s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVFreeWallInnerSide1 = mesh.createElementtemplate() - elementtemplateRVFreeWallInnerSide1.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVFreeWallInnerSide1.defineField(coordinates, -1, eftRVFreeWallInnerSide1) - #print('eftRVFreeWallInnerSide1', result) - - # RV free wall inner side 2 elements - eftRVFreeWallInnerSide2 = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftRVFreeWallInnerSide2, [1], [3, 4, 3, 4]) - for s in range(2): - n = s*2 + 1 - ln = n + 1 - # d/dxi1 = -d/ds3, d/dxi3 = sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftRVFreeWallInnerSide2, n*8 + 2, ln, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftRVFreeWallInnerSide2, n*8 + 5, ln, Node.VALUE_LABEL_D_DS1, 1, [s*2 + 2], Node.VALUE_LABEL_D_DS3, 1, [s*2 + 3]) - elementtemplateRVFreeWallInnerSide2 = mesh.createElementtemplate() - elementtemplateRVFreeWallInnerSide2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateRVFreeWallInnerSide2.defineField(coordinates, -1, eftRVFreeWallInnerSide2) - #print('eftRVFreeWallInnerSide2', result) + deltaRadians = dRadians1 = dradians2 = radiansPerElementAroundLVFreeWall + radiansAroundNext = radiansAround + deltaRadians + radians1 = -radiansAround + radians2 = -radiansAroundNext + scalefactors = [ -1.0, + math.cos(radians1), math.sin(radians1), dRadians1, + math.cos(radians2), math.sin(radians2), dRadians2, + math.cos(radians1), math.sin(radians1), dRadians1, + math.cos(radians2), math.sin(radians2), dRadians2 + ] + radiansAround = radiansAroundNext + else: + bnil = 2 + norl*(e2 - 1) + e1 + bnjl = 2 + norl*(e2 - 1) + (e1 + 1)%elementsCountAroundLV + nids = [ bnil , bnjl , bnil + norl , bnjl + norl, \ + bnil + nowl, bnjl + nowl, bnil + norl + nowl, bnjl + norl + nowl ] + if (e2 >= elementsCountUpApex) and (e1 < elementsCountAroundSeptum): + if (e2 == elementsCountUpApex) or (e1 == 0): + nids[4] = nidr + norr*(e2 - elementsCountUpApex) + e1 + if e1 == 0: + nids[6] = nids[4] + norr + if (e2 == elementsCountUpApex) or (e1 == (elementsCountAroundSeptum - 1)): + nids[5] = nidr + norr*(e2 - elementsCountUpApex) + e1 + 1 + if e1 == (elementsCountAroundSeptum - 1): + nids[7] = nids[5] + norr + if e2 == elementsCountUpApex: + if e1 == 0: + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 3 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 5 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 6 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + scalefactors = [ -1.0 ] + elif e1 == (elementsCountAroundSeptum - 1): + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 4 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 5 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 6 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 8 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + scalefactors = [ -1.0 ] + else: + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 1, 2 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 6 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + scalefactors = [ -1.0 ] + else: + if e1 == 0: + 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_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 7 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + scalefactors = [ -1.0 ] + elif e1 == (elementsCountAroundSeptum - 1): + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 6, 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 6, 8 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + scalefactors = [ -1.0 ] - for n2 in range(elementsCountUpRV): + result = elementtemplate1.defineField(coordinates, -1, eft1) + + element = mesh.createElement(elementIdentifier, elementtemplate1) + result2 = element.setNodesByIdentifier(eft1, nids) + if eft1.getNumberOfLocalScaleFactors() > 0: + result3 = element.setScaleFactors(eft1, scalefactors) + else: + result3 = 7 + #print('create element lv', elementIdentifier, result, result2, result3, nids) + elementIdentifier = elementIdentifier + 1 + + # RV elements + + nidl = 2 + (elementsCountUpApex - 1)*norl + + for e2 in range(elementsCountUpRV): + + for e1 in range(elementsCountAroundRV): - for n1 in range(elementsCountAcrossSeptum): - - bni = n2*(elementsCountAcrossSeptum + 1) + n1 - nodeIdentifiers = [ - rv_nidsInner[bni], rv_nidsInner[bni + 1], rv_nidsInner[bni + elementsCountAcrossSeptum + 1], rv_nidsInner[bni + elementsCountAcrossSeptum + 2], - rv_nidsOuter[bni], rv_nidsOuter[bni + 1], rv_nidsOuter[bni + elementsCountAcrossSeptum + 1], rv_nidsOuter[bni + elementsCountAcrossSeptum + 2] - ] - if n2 == 0: - if n1 == 0: - # RV free wall inner bottom side 1 elements: doubly curved but can't avoid a sharp corner at node 1 - eftTemp = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftTemp, [1], [3, 4, 3, 4]) - # node 2 d/dxi2 = d/ds3, d/dxi3 = -sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 1*8 + 3, 2, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftTemp, 1*8 + 5, 2, Node.VALUE_LABEL_D_DS2, 1, [1, 2], Node.VALUE_LABEL_D_DS3, 1, [3]) - # node 3 d/dxi1 = d/ds3, d/dxi3 = -sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 2*8 + 2, 3, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftTemp, 2*8 + 5, 3, Node.VALUE_LABEL_D_DS1, 1, [1, 4], Node.VALUE_LABEL_D_DS3, 1, [5]) - #print('inner bottom side 1 eftTemp.isValid()',eftTemp.isValid()) - elementtemplateTemp = mesh.createElementtemplate() - elementtemplateTemp.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateTemp.defineField(coordinates, -1, eftTemp) - element = mesh.createElement(elementIdentifier, elementtemplateTemp) - result = element.setNodesByIdentifier(eftTemp, nodeIdentifiers) - result = element.setScaleFactors(eftTemp, scalefactors5) - elif n1 == (elementsCountAcrossSeptum - 1): - # RV free wall inner bottom side 1 elements: doubly curved but can't avoid a sharp corner at node 1 - eftTemp = tricubichermite.createEftBasic() - setEftScaleFactorIds(eftTemp, [1], [3, 4, 3, 4]) - # node 1 d/dxi2 = d/ds3, d/dxi3 = -sa.d/ds2 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 0*8 + 3, 1, Node.VALUE_LABEL_D_DS3, 1, []) - mapEftFunction1Node2Terms(eftTemp, 0*8 + 5, 1, Node.VALUE_LABEL_D_DS2, 1, [1, 2], Node.VALUE_LABEL_D_DS3, 1, [3]) - # node 4 d/dxi1 = -d/ds3, d/dxi3 = sa.d/ds1 + sb.d/ds3 - mapEftFunction1Node1Term(eftTemp, 3*8 + 2, 4, Node.VALUE_LABEL_D_DS3, 1, [1]) - mapEftFunction1Node2Terms(eftTemp, 3*8 + 5, 4, Node.VALUE_LABEL_D_DS1, 1, [4], Node.VALUE_LABEL_D_DS3, 1, [5]) - #print('inner bottom side 2 eftTemp.isValid()',eftTemp.isValid()) - elementtemplateTemp = mesh.createElementtemplate() - elementtemplateTemp.setElementShapeType(Element.SHAPE_TYPE_CUBE) - result = elementtemplateTemp.defineField(coordinates, -1, eftTemp) - element = mesh.createElement(elementIdentifier, elementtemplateTemp) - result = element.setNodesByIdentifier(eftTemp, nodeIdentifiers) - result = element.setScaleFactors(eftTemp, scalefactors5) + eft1 = eft + scalefactors = None + + if e2 == 0: + if (e1 == 0) or (e1 == (elementsCountAroundRV - 1)): + # skip as two fewer elements on bottom row + continue + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + scalefactors = [ -1.0 ] + if e1 == 1: + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 3 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 4 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 7 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) + elif e1 == (elementsCountAroundRV - 2): + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 3 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 4 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 6, 8 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) else: - element = mesh.createElement(elementIdentifier, elementtemplateRVFreeWallInnerBottom) - result = element.setNodesByIdentifier(eftRVFreeWallInnerBottom, nodeIdentifiers) - result = element.setScaleFactors(eftRVFreeWallInnerBottom, scalefactors5) + remapEftNodeValueLabel(eft1, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 3, 4 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + bnil = nidl + norl*e2 + e1 - 1 + bnjl = bnil + 1 + bnir = nidr + norr*e2 + e1 - 1 + bnjr = bnir + 1 + nids = [ bnil , bnjl , bnir , bnjr, \ + bnil + nowl, bnjl + nowl, bnir + nowr, bnjr + nowr ] else: - if n1 == 0: - element = mesh.createElement(elementIdentifier, elementtemplateRVFreeWallInnerSide1) - result = element.setNodesByIdentifier(eftRVFreeWallInnerSide1, nodeIdentifiers) - result = element.setScaleFactors(eftRVFreeWallInnerSide1, scalefactors5) - elif n1 == (elementsCountAcrossSeptum - 1): - element = mesh.createElement(elementIdentifier, elementtemplateRVFreeWallInnerSide2) - result = element.setNodesByIdentifier(eftRVFreeWallInnerSide2, nodeIdentifiers) - result = element.setScaleFactors(eftRVFreeWallInnerSide2, scalefactors5) + if e1 == 0: + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + scalefactors = [ -1.0 ] + if e2 == 1: + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 3 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [] ) ]) + else: + remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + bnil = nidl + norl*(e2 - 1) + e1 + bnir = nidr + norr*(e2 - 1) + e1 + nids = [ bnil , bnir , bnil + norl , bnir + norr, \ + bnil + nowl, bnir + nowr, bnil + norl + nowl, bnir + norr + nowr ] + elif e1 == (elementsCountAroundRV - 1): + eft1 = tricubichermite.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft1, [1], []) + if e2 == 1: + remapEftNodeValueLabel(eft1, [ 1 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [1] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [1] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 3 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + else: + remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + bnil = nidl + norl*(e2 - 1) + elementsCountAroundSeptum + bnir = nidr + norr*(e2 - 1) + elementsCountAroundRV - 2 + nids = [ bnir , bnil , bnir + norr , bnil + norl, \ + bnir + nowr, bnil + nowl, bnir + norr + nowr, bnil + norl + nowl ] + scalefactors = [ -1.0 ] else: - element = mesh.createElement(elementIdentifier, elementtemplate) - result = element.setNodesByIdentifier(eft, nodeIdentifiers) - #print('RV free wall element create', elementIdentifier, result) - elementIdentifier += 1 + bnir = nidr + norr*(e2 - 1) + e1 - 1 + bnjr = bnir + 1 + nids = [ bnir , bnjr , bnir + norr , bnjr + norr, \ + bnir + nowr, bnjr + nowr, bnir + norr + nowr, bnjr + norr + nowr ] + + result = elementtemplate1.defineField(coordinates, -1, eft1) + + element = mesh.createElement(elementIdentifier, elementtemplate1) + result2 = element.setNodesByIdentifier(eft1, nids) + if eft1.getNumberOfLocalScaleFactors() > 0: + result3 = element.setScaleFactors(eft1, scalefactors) + else: + result3 = 7 + #print('create element rv', elementIdentifier, result, result2, result3, nids) + elementIdentifier = elementIdentifier + 1 + fm.endChange() + + + @staticmethod + def refineMesh(meshrefinement, options): + """ + Refine source mesh into separate region, with change of basis. + Stops at end of ventricles, hence can be called from ventriclesbase. + :param meshrefinement: MeshRefinement, which knows source and target region. + :param options: Dict containing options. See getDefaultOptions(). + """ + assert isinstance(meshrefinement, MeshRefinement) + elementsCountAroundLVFreeWall = options['Number of elements around LV free wall'] + elementsCountAroundSeptum = options['Number of elements around septum'] + elementsCountAroundLV = elementsCountAroundLVFreeWall + elementsCountAroundSeptum + elementsCountUpApex = options['Number of elements up apex'] + elementsCountUpSeptum = options['Number of elements up septum'] + elementsCountUp = elementsCountUpApex + elementsCountUpSeptum + refineElementsCountSurface = options['Refine number of elements surface'] + refineElementsCountThroughLVWall = options['Refine number of elements through LV wall'] + refineElementsCountThroughRVWall = options['Refine number of elements through RV wall'] + + elementsCountUpRV = elementsCountUpSeptum + 1 + elementsCountAroundRV = elementsCountAroundSeptum + 2 + startRvElementIdentifier = elementsCountAroundLV*elementsCountUp + 1 + limitRvElementIdentifier = startRvElementIdentifier + elementsCountUpRV*elementsCountAroundRV - 2 + + element = meshrefinement._sourceElementiterator.next() + while element.isValid(): + numberInXi1 = refineElementsCountSurface + numberInXi2 = refineElementsCountSurface + elementId = element.getIdentifier() + if elementId < startRvElementIdentifier: + numberInXi3 = refineElementsCountThroughLVWall + elif elementId < limitRvElementIdentifier: + numberInXi3 = refineElementsCountThroughRVWall + meshrefinement.refineElementCubeStandard3d(element, numberInXi1, numberInXi2, numberInXi3) + if elementId == (limitRvElementIdentifier - 1): + return # finish on last so can continue in ventriclesbase + element = meshrefinement._sourceElementiterator.next() + + @staticmethod + def generateMesh(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(). + """ + if not options['Refine']: + MeshType_3d_heartventricles2.generateBaseMesh(region, options) + return + baseRegion = region.createRegion() + MeshType_3d_heartventricles2.generateBaseMesh(baseRegion, options) + meshrefinement = MeshRefinement(baseRegion, region) + MeshType_3d_heartventricles2.refineMesh(meshrefinement, options) + + +def getSeptumPoints(septumArcRadians, lvRadius, radialDisplacement, elementsCountAroundLVFreeWall, elementsCountAroundSeptum, z, n3): + ''' + Symmetric 2-cubic interpolation of n septum elements around arc. + :return: x[], dx_ds1[] + ''' + radiansPerElementAroundLVFreeWall = (2.0*math.pi - septumArcRadians)/elementsCountAroundLVFreeWall + # get cubic curve with arc length scaling across to centre of septum + radiansAround = 0.5*septumArcRadians + circleArcLength = lvRadius*radiansAround + v1 = [ lvRadius - radialDisplacement, 0.0, z ] + d1 = [ 0.0, circleArcLength, 0.0 ] + v2 = [ lvRadius*math.cos(radiansAround), lvRadius*math.sin(radiansAround), z ] + d2 = [ -circleArcLength*math.sin(radiansAround), circleArcLength*math.cos(radiansAround), 0.0 ] + cubicArcLength = computeCubicHermiteArcLength(v1, d1, v2, d2, False) + scale = cubicArcLength/circleArcLength + d1 = [ d*scale for d in d1 ] + d2 = [ d*scale for d in d2 ] + elementLengthMid = (2.0*cubicArcLength - lvRadius*radiansPerElementAroundLVFreeWall)/(elementsCountAroundSeptum + n3 - 1) + #elementLengthMid = (2.0*cubicArcLength - lvRadius*radiansPerElementAroundLVFreeWall)/(elementsCountAroundSeptum + n3*0.5 - 1) + #elementLengthMid = (2.0*cubicArcLength - lvRadius*radiansPerElementAroundLVFreeWall)/(elementsCountAroundSeptum - 1) + #elementLengthEnd = 0.5*(lvRadius*radiansPerElementAroundLVFreeWall + elementLengthMid) + length = (elementsCountAroundSeptum % 2)*elementLengthMid*0.5 + x = [] + dx_ds1 = [] + for n1 in range(elementsCountAroundSeptum//2): + xi = length/cubicArcLength + pos1 = interpolateCubicHermite(v1, d1, v2, d2, xi) + x.append([ v for v in pos1 ]) + deriv1 = interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi) + scale = elementLengthMid/vector.magnitude(deriv1) + dx_ds1.append([ d*scale for d in deriv1 ]) + length += elementLengthMid + return x, dx_ds1 + + +def getRvOuterPoints(rvArcAroundRadians, lvRadius, rvAddWidthRadius, rvAddCrossRadius, elementsCountAroundRV, dEndMag, z, xi, mirror=False): + ''' + Get array of points and derivatives around half of RV from +x axis anticlockwise. + LV is assumed to be circular, centred at x, y = (0,0) + :param rvArcAroundRadians: Angle in radians around outside of RV to tangent nodes where it meets LV. + :param lvRadius: Radius of the LV. + :param rvAddWidthRadius: Radius to add to LV to get RV width at centre. + :param rvAddCrossRadius: Additional radius to add only laterally around septum. + :param elementsCountAroundRV: Number of elements around RV. + :param dEndMag: Magnitude of the external derivative at LV junction. + :param z: Z coordinate to give to all values of x[]. dx_ds1[2] is all zero. + :param xi: xi value ranging from 0 at bottom of RV to 1 at base. + :param mirror: Set to True to mirror points, otherwise only second half of RV returned. + :return: Arrays x[], dx_ds1[]. + ''' + elementsCountEnd = 0 if (xi <= 0.0) else 1 + #print('\ngetRvOuterPoints', rvArcAroundRadians, lvRadius, rvAddWidthRadius, rvAddCrossRadius, elementsCountAroundRV, dEndMag, z) + startRadians = 0.5*rvArcAroundRadians + ellipseEndRadians = startRadians - xi*0.5*rvArcAroundRadians/elementsCountAroundRV + b = lvRadius + rvAddCrossRadius # cross radius + if rvAddCrossRadius >= 0.0: + theta = math.pi - math.asin(lvRadius*math.sin(ellipseEndRadians)/b) + else: + theta = math.pi/2.0 + #print('phi',ellipseEndRadians) + #print('theta',theta) + a = (lvRadius*(math.cos(ellipseEndRadians) - 1.0) - rvAddWidthRadius)/(math.cos(theta) - 1.0) # width radius + cx = lvRadius + rvAddWidthRadius - a + #print(ellipseEndRadians,'a',a,'b',b,'cx',cx) + # get cubic curve joining half ellipse with to end point, first with unit derivatives then compute arc length + x1 = ( cx, b ) + d1 = ( -1.0, 0.0 ) + x2 = ( lvRadius*math.cos(startRadians), lvRadius*math.sin(startRadians) ) + d2 = ( -math.sin(startRadians), math.cos(startRadians) ) + cubicArcLength = computeCubicHermiteArcLength(x1, d1, x2, d2, True) + d1 = ( d1[0]*cubicArcLength, d1[1]*cubicArcLength ) + d2 = ( d2[0]*cubicArcLength, d2[1]*cubicArcLength ) + quarterEllipsePerimeter = 0.25*getApproximateEllipsePerimeter(a, b) + halfLength = cubicArcLength + quarterEllipsePerimeter + #print('getRvOuterPoints. length cubic', cubicArcLength, ', quarterEllipsePerimeter', quarterEllipsePerimeter, ', halfLength',halfLength) + + if elementsCountEnd == 0: + elementLengthTrans = elementLengthMid = (2.0*halfLength - dEndMag)/(elementsCountAroundRV - 1) + elementLengthEnd = 0.5*(elementLengthMid + dEndMag) + else: + elementLengthEnd = dEndMag + elementLengthMid = (2.0*halfLength - elementLengthEnd*3.0)/(elementsCountAroundRV - 3) + elementLengthTrans = 0.5*(elementLengthEnd + elementLengthMid) + #print('elementLengthMid',elementLengthMid,', elementLengthEnd',elementLengthEnd) + length = 0.0 + angle = 0.0 + if (elementsCountAroundRV % 2) == 1: + length += 0.5*elementLengthMid + angle = updateEllipseAngleByArcLength(a, b, angle, length) + x = [] + dx_ds1 = [] + nLimit = elementsCountAroundRV//2 + for n in range(nLimit): + dMag = elementLengthMid + if n < (nLimit - 2): + elementLength = elementLengthMid + elif n == (nLimit - 2): + elementLength = elementLengthTrans + else: + elementLength = elementLengthEnd + if elementsCountEnd > 0: + dMag = elementLengthEnd + if length <= quarterEllipsePerimeter: + cosAngle = math.cos(angle) + sinAngle = math.sin(angle) + v = ( cx + a*cosAngle, b*sinAngle ) + t = ( -a*sinAngle, b*cosAngle ) + angle = updateEllipseAngleByArcLength(a, b, angle, elementLength) + else: + xi = (length - quarterEllipsePerimeter)/cubicArcLength + v = interpolateCubicHermite(x1, d1, x2, d2, xi) + t = interpolateCubicHermiteDerivative(x1, d1, x2, d2, xi) + x.append([ v[0], v[1], z ]) + scale = dMag/math.sqrt(t[0]*t[0] + t[1]*t[1]) + dx_ds1.append([ t[0]*scale, t[1]*scale, 0.0 ]) + length += elementLength + #print('getRvOuterPoints. length',length,', overshoot', length - halfLength) + if mirror: + xm = [] + dxm_ds1 = [] + for n in range(elementsCountAroundRV//2 - 1, -(elementsCountAroundRV%2), -1): + xm.append([ x[n][0], -x[n][1], x[n][2] ]) + dxm_ds1.append([ -dx_ds1[n][0], dx_ds1[n][1], dx_ds1[n][2] ]) + x = xm + x + dx_ds1 = dxm_ds1 + dx_ds1 + return x, dx_ds1 + + +def getRVOuterSize(xiUpWidth, xiUpCross, rvWidthRadius, rvExtraCrossRadiusBase): + if xiUpWidth < 0.0: + #print('getRVOuterSize NEGATIVE xiUpWidth', xiUpWidth) + return 0.0, 0.0 + if xiUpCross < 0.0: + xiUpCross = 0.0 + xiUpFast = 1.0 - (1.0 - xiUpWidth)*(1.0 - xiUpWidth) + xiUpSlow = xiUpCross + rvAddWidthRadius = interpolateCubicHermite([0.0], [0.0], [rvWidthRadius], [0.0], xiUpFast)[0] + rvAddCrossRadius = interpolateCubicHermite([0.0], [0.0], [rvExtraCrossRadiusBase], [0.0], xiUpSlow)[0] + #print('getRVOuterSize xiUpWidth', xiUpWidth, ', addWidth', rvAddWidthRadius, ', addCross', rvAddCrossRadius) + return rvAddWidthRadius, rvAddCrossRadius diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index c4c1e50c..2f645398 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -6,27 +6,12 @@ ''' from scaffoldmaker.utils.eft_utils import * from scaffoldmaker.utils.zinc_utils import * +import scaffoldmaker.utils.vector as vector from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate from opencmiss.zinc.node import Node from opencmiss.zinc.status import OK as ZINC_OK import math -def normalise(v): - ''' - :return: vector v normalised to unit length - ''' - mag = 0.0 - for s in v: - mag += s*s - mag = math.sqrt(mag) - return [ s/mag for s in v ] - -def crossproduct3(a, b): - ''' - :return: vector 3-D cross product of a and b - ''' - return [ a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0] ] - class eftfactory_tricubichermite: ''' Factory class for creating element field templates for a 3-D mesh using tricubic Hermite basis. @@ -245,6 +230,20 @@ def createEftSplitXi1RightOut(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftSplitXi1RightOut: Failed to validate eft' return eft + def createEftSplitXi2RightStraight(self): + ''' + Create an element field template suitable for the inner elements of the + join between left and right chambers at the bottom of the RV, with xi2 bifurcating. + Straight through version. + :return: Element field template + ''' + eft = self.createEftNoCrossDerivatives() + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, [ 5, 6 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, []), (Node.VALUE_LABEL_D_DS3, [1]) ]) + remapEftNodeValueLabel(eft, [ 7, 8 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, []), (Node.VALUE_LABEL_D_DS3, []) ]) + assert eft.validate(), 'eftfactory_tricubichermite.createEftSplitXi2RightStraight: Failed to validate eft' + return eft + def createEftTubeSeptumOuter(self): ''' Create an element field template suitable for the outer elements of @@ -539,13 +538,13 @@ def replaceElementWithInlet4(self, element, startElementId, nodetemplate, startN result, fc = coordinates.evaluateReal(cache, 3) resulta, a = coordinates.evaluateDerivative(diff1, cache, 3) resultb, b = coordinates.evaluateDerivative(diff2, cache, 3) - n = normalise(crossproduct3(a, b)) + n = vector.normalise(vector.crossproduct3(a, b)) #print(resulta, 'a =', a, ',', resultb, ' b=', b, ' fc=', fc, ' n=',n) ic = [ (fc[i] + tubeLength*n[i]) for i in range(3) ] - na = normalise(a) - nb = normalise(b) - a = normalise([ -(na[i] + nb[i]) for i in range(3) ]) - b = normalise(crossproduct3(a, n)) + na = vector.normalise(a) + nb = vector.normalise(b) + a = vector.normalise([ -(na[i] + nb[i]) for i in range(3) ]) + b = vector.normalise(vector.crossproduct3(a, n)) zero = [ 0.0, 0.0, 0.0 ] nodeIdentifier = startNodeId diff --git a/scaffoldmaker/utils/geometry.py b/scaffoldmaker/utils/geometry.py new file mode 100644 index 00000000..f60c3600 --- /dev/null +++ b/scaffoldmaker/utils/geometry.py @@ -0,0 +1,68 @@ +''' +Utility functions for geometry. +Created on Apr 11, 2018 + +@author: Richard Christie +''' + +import math + +def getApproximateEllipsePerimeter(a, b): + ''' + Get perimeter of ellipse using Ramanujan II approximation. + :param a: Major axis length. + :param b: Minor axis length. + :return: Perimeter length. + ''' + h = ((a-b)/(a+b))**2 + return math.pi*(a + b)*(1.0 + 3.0*h/(10.0 + math.sqrt(4.0 - 3.0*h))) + +def getEllipseArcLength(a, b, angle1Radians, angle2Radians): + ''' + Calculates perimeter distance between two angles by summing line segments at regular angles. + :param a: Major axis length (On x, 0 / PI). + :param b: Minor axis length.(On y, PI/2, 3PI/2). + :param angle1Radians: First angle anticlockwise from major axis. + :param angle2Radians: Second angle anticlockwise from major axis. + :return: Perimeter lenge, positive if anticlockwise, otherwise negative. + ''' + angle1 = min(angle1Radians, angle2Radians) + angle2 = max(angle1Radians, angle2Radians) + # Max 100 segments around ellipse + segmentCount = int(math.ceil(50*(angle2-angle1)/math.pi)) + length = 0.0 + for i in range(segmentCount + 1): + r = i/segmentCount + angle = r*angle1 + (1.0 - r)*angle2 + x = ( a*math.cos(angle), b*math.sin(angle) ) + if i > 0: + delta = math.sqrt((x[0] - lastX[0])*(x[0] - lastX[0]) + (x[1] - lastX[1])*(x[1] - lastX[1])) + length += delta + lastX = x + if angle1Radians < angle2Radians: + return length + else: + return -length + +def updateEllipseAngleByArcLength(a, b, inAngleRadians, arcLength): + ''' + Update angle around ellipse to subtend arcLength around the perimeter. + Iterates using Newton's method. + :param inAngleRadians: Initial angle anticlockwise from major axis. + :param arcLength: Arc length to traverse. Positive=anticlockwise, negative=clockwise. + :param a: Major axis length (On x, 0 / PI). + :param b: Minor axis length.(On y, PI/2, 3PI/2). + :return: New angle, in radians. + ''' + angle = inAngleRadians + lengthMoved = 0.0 + lengthTol = (a + b)*1.0E-4 # broader tolerance due to reliance on inexact getEllipseArcLength() + #print('inAngleRadians', inAngleRadians, ', arcLength', arcLength) + while math.fabs(arcLength - lengthMoved) > lengthTol: + t = ( -a*math.sin(angle), b*math.cos(angle) ) + dlength_dangle = math.sqrt(t[0]*t[0] + t[1]*t[1]) + angle += (arcLength - lengthMoved)/dlength_dangle + lengthMoved = getEllipseArcLength(a, b, inAngleRadians, angle) + #print('lengthMoved', lengthMoved) + #print('updateEllipseAngleByArcLength a', a, 'b', b, ', angle', inAngleRadians, ', arcLength', arcLength, ' -> ', angle) + return angle diff --git a/scaffoldmaker/utils/interpolation.py b/scaffoldmaker/utils/interpolation.py index a9095a0d..55aaeca0 100644 --- a/scaffoldmaker/utils/interpolation.py +++ b/scaffoldmaker/utils/interpolation.py @@ -5,6 +5,12 @@ @author: Richard Christie ''' +import math +import scaffoldmaker.utils.vector as vector + +gaussXi3 = ( (-math.sqrt(0.6)+1.0)/2.0, 0.5, (+math.sqrt(0.6)+1.0)/2.0 ) +gaussWt3 = ( 5.0/18.0, 4.0/9.0, 5.0/18.0 ) + def interpolateCubicHermite(v1, d1, v2, d2, xi): """ Return cubic Hermite interpolated value of tuples v1, d1 (end 1) to v2, d2 (end 2) for xi in [0,1] @@ -16,10 +22,7 @@ def interpolateCubicHermite(v1, d1, v2, d2, xi): f2 = xi - 2.0*xi2 + xi3 f3 = 3.0*xi2 - 2.0*xi3 f4 = -xi2 + xi3 - result = [] - for i in range(len(v1)): - result.append(f1*v1[i] + f2*d1[i] + f3*v2[i] + f4*d2[i]) - return tuple(result) + return tuple([ (f1*v1[i] + f2*d1[i] + f3*v2[i] + f4*d2[i]) for i in range(len(v1)) ]) def interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi): """ @@ -31,7 +34,69 @@ def interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi): f2 = 1.0 - 4.0*xi + 3.0*xi2 f3 = 6.0*xi - 6.0*xi2 f4 = -2.0*xi + 3.0*xi2 - result = [] - for i in range(len(v1)): - result.append(f1*v1[i] + f2*d1[i] + f3*v2[i] + f4*d2[i]) - return tuple(result) + return tuple([ (f1*v1[i] + f2*d1[i] + f3*v2[i] + f4*d2[i]) for i in range(len(v1)) ]) + +def interpolateCubicHermiteSecondDerivative(v1, d1, v2, d2, xi): + """ + Return cubic Hermite interpolated second derivatives of tuples v1, d1 (end 1) to v2, d2 (end 2) for xi in [0,1] + :return: tuple containing result + """ + f1 = -6.0 + 12.0*xi + f2 = -4.0 + 6.0*xi + f3 = 6.0 - 12.0*xi + f4 = -2.0 + 6.0*xi + return tuple([ (f1*v1[i] + f2*d1[i] + f3*v2[i] + f4*d2[i]) for i in range(len(v1)) ]) + +def computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives): + """ + Compute arc length between v1 and v2, scaling unit d1 and d2. + Iterative; not optimised. + :param d1: Initial derivative at v1. + :param d2: Initial derivative at v2. + :param rescaleDerivatives: If True, rescale initial d1 and d2 to |v2 - v| + :return: Arc length. + """ + if rescaleDerivatives: + lastArcLength = math.sqrt(sum((v2[i] - v1[i])*(v2[i] - v1[i]) for i in range(len(v1)))) + else: + lastArcLength = getCubicHermiteArcLength(v1, d1, v2, d2) + d1 = vector.normalise(d1) + d2 = vector.normalise(d2) + tol = 1.0E-6 + for iters in range(100): + #print('iter',iters,'=',lastArcLength) + d1s = [lastArcLength*d for d in d1] + d2s = [lastArcLength*d for d in d2] + arcLength = getCubicHermiteArcLength(v1, d1s, v2, d2s) + if iters > 9: + arcLength = 0.8*arcLength + 0.2*lastArcLength + if math.fabs(arcLength - lastArcLength) < tol*arcLength: + #print('computeCubicHermiteArcLength converged at iter',iters,'=',arcLength,', closeness', math.fabs(arcLength - lastArcLength)) + return arcLength + lastArcLength = arcLength + print('computeCubicHermiteArcLength Max iters reached:',iters,'=',arcLength,', closeness', math.fabs(arcLength - lastArcLength)) + return arcLength + +def getCubicHermiteArcLength(v1, d1, v2, d2): + ''' + :return: Arc length of cubic curve using 3 point Gaussian quadrature. + ''' + arcLength = 0.0 + for i in range(3): + dm = interpolateCubicHermiteDerivative(v1, d1, v2, d2, gaussXi3[i]) + arcLength += gaussWt3[i]*math.sqrt(sum(d*d for d in dm)) + return arcLength + +def getCubicHermiteCurvature(v1, d1, v2, d2, radialVector, xi): + """ + :param radialVector: Radial direction, assumed unit normal to curve tangent at point. + :return: Scalar curvature (1/R) of the 1-D cubic Hermite curve. + """ + tangent = interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi) + dTangent = interpolateCubicHermiteSecondDerivative(v1, d1, v2, d2, xi) + #tangentVector = vector.normalise(tangent) + #tangentCurvature = vector.dotproduct(dTangent, tangentVector) + radialCurvature = vector.dotproduct(dTangent, radialVector) + magTangent = vector.magnitude(tangent) + curvature = radialCurvature/(magTangent*magTangent) + return curvature diff --git a/scaffoldmaker/utils/vector.py b/scaffoldmaker/utils/vector.py new file mode 100644 index 00000000..7aaca964 --- /dev/null +++ b/scaffoldmaker/utils/vector.py @@ -0,0 +1,33 @@ +''' +Utility functions for vectors. +Created on Apr 13, 2018 + +@author: Richard Christie +''' + +import math + +def crossproduct3(a, b): + ''' + :return: vector 3-D cross product of a and b + ''' + return [ a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0] ] + +def dotproduct(a, b): + ''' + :return: vector dot (inner) product of a and b + ''' + return sum(a[i]*b[i] for i in range(len(a))) + +def magnitude(v): + ''' + return: scalar magnitude of vector v + ''' + return math.sqrt(sum(c*c for c in v)) + +def normalise(v): + ''' + :return: vector v normalised to unit length + ''' + mag = math.sqrt(sum(c*c for c in v)) + return [ c/mag for c in v ]