From 379651c310834a53a244e26356dcddc52d37a45a Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 19 Jun 2018 15:58:42 +1200 Subject: [PATCH 1/9] Add initial solid sphere --- .../meshtypes/meshtype_3d_solidsphere1.py | 341 ++++++++++++++++++ scaffoldmaker/scaffoldmaker.py | 5 +- 2 files changed, 345 insertions(+), 1 deletion(-) create mode 100644 scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py new file mode 100644 index 00000000..85318fa0 --- /dev/null +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -0,0 +1,341 @@ +""" +Generates a 3-D unit solid sphere mesh with variable numbers of elements +around, up and through the thickness. +""" + +from __future__ import division +import math +from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite +from scaffoldmaker.utils.zinc_utils import * +from scaffoldmaker.utils.interpolation import * +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate +from opencmiss.zinc.field import Field +from opencmiss.zinc.node import Node + +class MeshType_3d_solidsphere1: + ''' + classdocs + ''' + @staticmethod + def getName(): + return '3D Solid Sphere 1' + + @staticmethod + def getDefaultOptions(): + return { + 'Number of elements around' : 8, + 'Number of elements up' : 8, + 'Number of elements radial' : 1, + 'Diameter' : 1, + 'Use cross derivatives' : False, + 'Refine' : False, + 'Refine number of elements around' : 1, + 'Refine number of elements up' : 1, + 'Refine number of elements radial' : 1 + } + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Number of elements around', + 'Number of elements up', + 'Number of elements radial', + 'Diameter', + 'Use cross derivatives', + 'Refine', + 'Refine number of elements around', + 'Refine number of elements up', + 'Refine number of elements radial' + ] + + @staticmethod + def checkOptions(options): + for key in [ + 'Number of elements radial', + 'Refine number of elements around', + 'Refine number of elements up', + 'Refine number of elements radial']: + if options[key] < 1: + options[key] = 1 + 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['Diameter'] < 0: + options['Diameter'] = 1 + + + @staticmethod + 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 + """ + elementsCountAround = options['Number of elements around'] + elementsCountUp = options['Number of elements up'] + elementsCountRadial = options['Number of elements radial'] + useCrossDerivatives = options['Use cross derivatives'] + diameter = options['Diameter'] + radius = diameter/2 + + fm = region.getFieldmodule() + fm.beginChange() + coordinates = getOrCreateCoordinateField(fm) + + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplateApex = nodes.createNodetemplate() + nodetemplateApex.defineField(coordinates) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + if useCrossDerivatives: + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) + else: + nodetemplate = nodetemplateApex + + mesh = fm.findMeshByDimension(3) + + tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + eft = tricubichermite.createEftBasic() + + tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + + cache = fm.createFieldcache() + + # Create nodes + nodeIdentifier = 1 + radiansPerElementAround = 2.0*math.pi/elementsCountAround + radiansPerElementUp = math.pi/elementsCountUp + + x = [ 0.0, 0.0, 0.0 ] + dx_ds1 = [ 0.0, 0.0, 0.0 ] + dx_ds2 = [ 0.0, 0.0, 0.0 ] + dx_ds3 = [ 0.0, 0.0, 0.0 ] + zero = [ 0.0, 0.0, 0.0 ] + + cubicArcLengthList = [0]*(elementsCountUp+1) + + # Pre-calculate cubicArcLength along elementsCountUp + for n2 in range(1,elementsCountUp + 1): + radiansUp = n2*radiansPerElementUp + cosRadiansUp = math.cos(radiansUp) + sinRadiansUp = math.sin(radiansUp) + + # Calculate cubic hermite arclength linking point on axis to surface on sphere + v1 = [0,0,-radius+n2*2*radius/elementsCountUp] + d1 = [0,1,0] + v2 = [ + radius*math.cos(math.pi/2)*sinRadiansUp, + radius*math.sin(math.pi/2)*sinRadiansUp, + -radius*cosRadiansUp + ] + d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] + cubicArcLengthList[n2] = computeCubicHermiteArcLength(v1, d1, v2, d2, True) + + # Create node for bottom pole + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, -radius*2/elementsCountUp]) # height of element along central axis + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + + # Create nodes along axis between top and bottom poles + for n2 in range(1,elementsCountUp): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius+n2*2*radius/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -cubicArcLengthList[n2]/elementsCountRadial, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ cubicArcLengthList[n2]/elementsCountRadial, 0.0, 0.0 ]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + + # Create nodes for top pole + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, radius]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + + # Create other nodes + for n3 in range(1,elementsCountRadial+1): + xi = 1/elementsCountRadial*n3 + radiansUpArcOriginList = [0]*(elementsCountUp) + + # Pre-calculate RC for points on vertical arc running between top and bottom poles + pt = [0, radius*xi, 0] + arcOrigin = (radius*radius - pt[2]*pt[2] - pt[1]*pt[1])/(-2*pt[1]) + RC = math.sqrt(arcOrigin*arcOrigin + radius*radius) + + radiansUpArcOriginList[0] = math.acos(-radius/RC) + + # Identify nodes on the vertical arc using radiansAround = pi/2 + for n2 in range(1,elementsCountUp): + radiansUp = n2*radiansPerElementUp + cosRadiansUp = math.cos(radiansUp) + sinRadiansUp = math.sin(radiansUp) + + # Calculate node coordinates on arc using cubic hermite interpolation + cubicArcLength = cubicArcLengthList[n2] + v1 = [0,0,-radius+n2*2*radius/elementsCountUp] + d1 = [math.cos(math.pi/2),math.sin(math.pi/2),0] + d1 = vector.normalise(d1) + d1 = [d*cubicArcLength for d in d1] + v2 = [ + radius*math.cos(math.pi/2)*sinRadiansUp, + radius*math.sin(math.pi/2)*sinRadiansUp, + -radius*cosRadiansUp + ] + d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] + d2 = vector.normalise(d2) + d2 = [d*cubicArcLength for d in d2] + x = list(interpolateCubicHermite(v1, d1, v2, d2, xi)) + + # Calculate radiansUp for each point wrt arcOrigin + radiansUpArcOriginList[n2] = math.acos(x[2]/RC) + + + for n2 in range(1,elementsCountUp): + radiansUp = n2*radiansPerElementUp + cosRadiansUp = math.cos(radiansUp) + sinRadiansUp = math.sin(radiansUp) + + for n1 in range(elementsCountAround): + radiansAround = n1*radiansPerElementAround + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + cubicArcLength = cubicArcLengthList[n2] + + # Calculate node coordinates on arc using cubic hermite interpolation + v1 = [0,0,-radius+n2*2*radius/elementsCountUp] + d1 = [cosRadiansAround,sinRadiansAround,0] + d1 = vector.normalise(d1) + d1 = [d*cubicArcLength for d in d1] + v2 = [ + radius*cosRadiansAround*sinRadiansUp, + radius*sinRadiansAround*sinRadiansUp, + -radius*cosRadiansUp + ] + d2 = [cosRadiansAround*sinRadiansUp,sinRadiansAround*sinRadiansUp,-cosRadiansUp] + d2 = vector.normalise(d2) + d2 = [d*cubicArcLength for d in d2] + x = list(interpolateCubicHermite(v1, d1, v2, d2, xi)) + + # For dx_ds1 - Calculate radius wrt origin where interpolated points lie on + orthoRadius = vector.magnitude(x) + orthoRadiansUp = math.pi - math.acos(x[2]/orthoRadius) + sinOrthoRadiansUp = math.sin(orthoRadiansUp) + cosOrthoRadiansUp = math.cos(orthoRadiansUp) + + # For dx_ds2 - Assign radiansUp from radiansUpArcOriginList and calculate diff between radiansUp as we move up + radiansUpArcOrigin = radiansUpArcOriginList[n2] + sinRadiansUpArcOrigin = math.sin(radiansUpArcOrigin) + cosRadiansUpArcOrigin = math.cos(radiansUpArcOrigin) + radiansPerElementUpArcOrigin = radiansUpArcOriginList[n2]-radiansUpArcOriginList[n2-1] + + dx_ds1 = [ + orthoRadius*-sinRadiansAround*sinOrthoRadiansUp*radiansPerElementAround, + orthoRadius*cosRadiansAround*sinOrthoRadiansUp*radiansPerElementAround, + 0.0 + ] + + dx_ds2 = [ + RC*cosRadiansAround*cosRadiansUpArcOrigin*radiansPerElementUpArcOrigin, + RC*sinRadiansAround*cosRadiansUpArcOrigin*radiansPerElementUpArcOrigin, + -RC*sinRadiansUpArcOrigin*radiansPerElementUpArcOrigin + ] + + dx_ds3 = list(interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi)) + dx_ds3 = vector.normalise(dx_ds3) + dx_ds3 = [d*cubicArcLength/elementsCountRadial for d in dx_ds3] + + 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) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + nodeIdentifier = nodeIdentifier + 1 + + # create elements + elementIdentifier = 1 + no2 = elementsCountAround + no3 = elementsCountAround*(elementsCountUp - 1) + rni = (1 + elementsCountUp) - no3 - no2 + 1 # regular node identifier + for e3 in range(1, elementsCountRadial): + for e2 in range(1, elementsCountUp - 1): + bni = rni + e3*no3 + e2*no2 + for e1 in range(elementsCountAround): + element = mesh.createElement(elementIdentifier, elementtemplate) + na = e1 + nb = (e1 + 1)%elementsCountAround + nodeIdentifiers = [ + bni + na, bni + nb, bni + no2 + na, bni + no2 + nb, + bni + no3 + na, bni + no3 + nb, bni + no3 + no2 + na, bni + no3 + no2 + nb + ] + result = element.setNodesByIdentifier(eft, nodeIdentifiers) + # print('regular element', elementIdentifier, result, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + + fm.endChange() + + @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_solidsphere1.generateBaseMesh(region, options) + return + + refineElementsCountAround = options['Refine number of elements around'] + refineElementsCountUp = options['Refine number of elements up'] + refineElementsCountRadial = options['Refine number of elements radial'] + + baseRegion = region.createRegion() + MeshType_3d_solidsphere1.generateBaseMesh(baseRegion, options) + + meshrefinement = MeshRefinement(baseRegion, region) + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountUp, refineElementsCountRadial) diff --git a/scaffoldmaker/scaffoldmaker.py b/scaffoldmaker/scaffoldmaker.py index 4ff18994..e549b8cf 100644 --- a/scaffoldmaker/scaffoldmaker.py +++ b/scaffoldmaker/scaffoldmaker.py @@ -20,6 +20,8 @@ from scaffoldmaker.meshtypes.meshtype_3d_sphereshellseptum1 import MeshType_3d_sphereshellseptum1 from scaffoldmaker.meshtypes.meshtype_3d_tube1 import MeshType_3d_tube1 from scaffoldmaker.meshtypes.meshtype_3d_tubeseptum1 import MeshType_3d_tubeseptum1 +from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 + class Scaffoldmaker(object): @@ -41,7 +43,8 @@ def __init__(self): MeshType_3d_sphereshell1, MeshType_3d_sphereshellseptum1, MeshType_3d_tube1, - MeshType_3d_tubeseptum1 + MeshType_3d_tubeseptum1, + MeshType_3d_solidsphere1 ] def getMeshTypes(self): From 599bdcd41f4e2fb6fa44e52dc2bb31b5fb106e95 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 24 Jul 2018 13:21:56 +1200 Subject: [PATCH 2/9] Create solid sphere with 1 radial layer --- .../meshtypes/meshtype_3d_solidsphere1.py | 361 +++++++++++++----- scaffoldmaker/scaffoldmaker.py | 6 +- scaffoldmaker/utils/eft_utils.py | 4 +- .../utils/eftfactory_tricubichermite.py | 191 +++++++++ 4 files changed, 459 insertions(+), 103 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 85318fa0..6ba84df7 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -1,14 +1,14 @@ """ Generates a 3-D unit solid sphere mesh with variable numbers of elements -around, up and through the thickness. +around, up the central axis, and radially. """ from __future__ import division import math from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite -from scaffoldmaker.utils.zinc_utils import * from scaffoldmaker.utils.interpolation import * 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 @@ -25,14 +25,14 @@ def getName(): def getDefaultOptions(): return { 'Number of elements around' : 8, - 'Number of elements up' : 8, + 'Number of elements up' : 8, 'Number of elements radial' : 1, 'Diameter' : 1, 'Use cross derivatives' : False, 'Refine' : False, 'Refine number of elements around' : 1, 'Refine number of elements up' : 1, - 'Refine number of elements radial' : 1 + 'Refine number of elements radial' : 1 } @staticmethod @@ -40,21 +40,21 @@ def getOrderedOptionNames(): return [ 'Number of elements around', 'Number of elements up', - 'Number of elements radial', + 'Number of elements radial', 'Diameter', 'Use cross derivatives', 'Refine', 'Refine number of elements around', - 'Refine number of elements up', + 'Refine number of elements up', 'Refine number of elements radial' ] @staticmethod def checkOptions(options): - for key in [ + for key in [ 'Number of elements radial', 'Refine number of elements around', - 'Refine number of elements up', + 'Refine number of elements up', 'Refine number of elements radial']: if options[key] < 1: options[key] = 1 @@ -64,8 +64,7 @@ def checkOptions(options): options['Number of elements around'] = 4 if options['Diameter'] < 0: options['Diameter'] = 1 - - + @staticmethod def generateBaseMesh(region, options): """ @@ -77,10 +76,10 @@ def generateBaseMesh(region, options): elementsCountAround = options['Number of elements around'] elementsCountUp = options['Number of elements up'] elementsCountRadial = options['Number of elements radial'] - useCrossDerivatives = options['Use cross derivatives'] + useCrossDerivatives = options['Use cross derivatives'] diameter = options['Diameter'] radius = diameter/2 - + fm = region.getFieldmodule() fm.beginChange() coordinates = getOrCreateCoordinateField(fm) @@ -106,168 +105,160 @@ def generateBaseMesh(region, options): else: nodetemplate = nodetemplateApex - mesh = fm.findMeshByDimension(3) - - tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) - eft = tricubichermite.createEftBasic() - - tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) - - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate.defineField(coordinates, -1, eft) - cache = fm.createFieldcache() + ################# # Create nodes + ################# + nodeIdentifier = 1 radiansPerElementAround = 2.0*math.pi/elementsCountAround - radiansPerElementUp = math.pi/elementsCountUp - + radiansPerElementUp = math.pi/elementsCountUp + x = [ 0.0, 0.0, 0.0 ] dx_ds1 = [ 0.0, 0.0, 0.0 ] dx_ds2 = [ 0.0, 0.0, 0.0 ] dx_ds3 = [ 0.0, 0.0, 0.0 ] zero = [ 0.0, 0.0, 0.0 ] - + cubicArcLengthList = [0]*(elementsCountUp+1) - - # Pre-calculate cubicArcLength along elementsCountUp - for n2 in range(1,elementsCountUp + 1): + + # Pre-calculate cubicArcLength along elementsCountUp + for n2 in range(1,elementsCountUp + 1): radiansUp = n2*radiansPerElementUp cosRadiansUp = math.cos(radiansUp) sinRadiansUp = math.sin(radiansUp) # Calculate cubic hermite arclength linking point on axis to surface on sphere - v1 = [0,0,-radius+n2*2*radius/elementsCountUp] + v1 = [0,0,-radius+n2*2*radius/elementsCountUp] d1 = [0,1,0] v2 = [ - radius*math.cos(math.pi/2)*sinRadiansUp, + radius*math.cos(math.pi/2)*sinRadiansUp, radius*math.sin(math.pi/2)*sinRadiansUp, -radius*cosRadiansUp - ] + ] d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] cubicArcLengthList[n2] = computeCubicHermiteArcLength(v1, d1, v2, d2, True) - - # Create node for bottom pole + + # Create node for bottom pole node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, -radius*2/elementsCountUp]) # height of element along central axis - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) nodeIdentifier = nodeIdentifier + 1 - - # Create nodes along axis between top and bottom poles + + # Create nodes along axis between top and bottom poles for n2 in range(1,elementsCountUp): node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius+n2*2*radius/elementsCountUp]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -cubicArcLengthList[n2]/elementsCountRadial, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ cubicArcLengthList[n2]/elementsCountRadial, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ cubicArcLengthList[n2]/elementsCountRadial, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, cubicArcLengthList[n2]/elementsCountRadial, 0.0 ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) nodeIdentifier = nodeIdentifier + 1 - - # Create nodes for top pole + + # Create nodes for top pole node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, radius]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ 0.0, -radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) - nodeIdentifier = nodeIdentifier + 1 - - # Create other nodes - for n3 in range(1,elementsCountRadial+1): - xi = 1/elementsCountRadial*n3 + nodeIdentifier = nodeIdentifier + 1 + + # Create other nodes + for n3 in range(1,elementsCountRadial+1): + xi = 1/elementsCountRadial*n3 radiansUpArcOriginList = [0]*(elementsCountUp) - - # Pre-calculate RC for points on vertical arc running between top and bottom poles + + # Pre-calculate RC for points on vertical arc running between top and bottom poles pt = [0, radius*xi, 0] - arcOrigin = (radius*radius - pt[2]*pt[2] - pt[1]*pt[1])/(-2*pt[1]) - RC = math.sqrt(arcOrigin*arcOrigin + radius*radius) - + arcOrigin = (radius*radius - pt[2]*pt[2] - pt[1]*pt[1])/(-2*pt[1]) + RC = math.sqrt(arcOrigin*arcOrigin + radius*radius) + radiansUpArcOriginList[0] = math.acos(-radius/RC) - - # Identify nodes on the vertical arc using radiansAround = pi/2 + + # Identify nodes on the vertical arc using radiansAround = pi/2 for n2 in range(1,elementsCountUp): radiansUp = n2*radiansPerElementUp cosRadiansUp = math.cos(radiansUp) sinRadiansUp = math.sin(radiansUp) - - # Calculate node coordinates on arc using cubic hermite interpolation - cubicArcLength = cubicArcLengthList[n2] + + # Calculate node coordinates on arc using cubic hermite interpolation + cubicArcLength = cubicArcLengthList[n2] v1 = [0,0,-radius+n2*2*radius/elementsCountUp] d1 = [math.cos(math.pi/2),math.sin(math.pi/2),0] d1 = vector.normalise(d1) - d1 = [d*cubicArcLength for d in d1] + d1 = [d*cubicArcLength for d in d1] v2 = [ radius*math.cos(math.pi/2)*sinRadiansUp, radius*math.sin(math.pi/2)*sinRadiansUp, -radius*cosRadiansUp - ] + ] d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] d2 = vector.normalise(d2) d2 = [d*cubicArcLength for d in d2] x = list(interpolateCubicHermite(v1, d1, v2, d2, xi)) - - # Calculate radiansUp for each point wrt arcOrigin + + # Calculate radiansUp for each point wrt arcOrigin radiansUpArcOriginList[n2] = math.acos(x[2]/RC) - - + for n2 in range(1,elementsCountUp): radiansUp = n2*radiansPerElementUp cosRadiansUp = math.cos(radiansUp) - sinRadiansUp = math.sin(radiansUp) - + sinRadiansUp = math.sin(radiansUp) + for n1 in range(elementsCountAround): - radiansAround = n1*radiansPerElementAround + radiansAround = n1*radiansPerElementAround cosRadiansAround = math.cos(radiansAround) - sinRadiansAround = math.sin(radiansAround) - cubicArcLength = cubicArcLengthList[n2] - - # Calculate node coordinates on arc using cubic hermite interpolation + sinRadiansAround = math.sin(radiansAround) + cubicArcLength = cubicArcLengthList[n2] + + # Calculate node coordinates on arc using cubic hermite interpolation v1 = [0,0,-radius+n2*2*radius/elementsCountUp] d1 = [cosRadiansAround,sinRadiansAround,0] d1 = vector.normalise(d1) - d1 = [d*cubicArcLength for d in d1] + d1 = [d*cubicArcLength for d in d1] v2 = [ radius*cosRadiansAround*sinRadiansUp, radius*sinRadiansAround*sinRadiansUp, -radius*cosRadiansUp - ] + ] d2 = [cosRadiansAround*sinRadiansUp,sinRadiansAround*sinRadiansUp,-cosRadiansUp] d2 = vector.normalise(d2) d2 = [d*cubicArcLength for d in d2] x = list(interpolateCubicHermite(v1, d1, v2, d2, xi)) - - # For dx_ds1 - Calculate radius wrt origin where interpolated points lie on + + # For dx_ds1 - Calculate radius wrt origin where interpolated points lie on orthoRadius = vector.magnitude(x) - orthoRadiansUp = math.pi - math.acos(x[2]/orthoRadius) + orthoRadiansUp = math.pi - math.acos(x[2]/orthoRadius) sinOrthoRadiansUp = math.sin(orthoRadiansUp) - cosOrthoRadiansUp = math.cos(orthoRadiansUp) - + cosOrthoRadiansUp = math.cos(orthoRadiansUp) + # For dx_ds2 - Assign radiansUp from radiansUpArcOriginList and calculate diff between radiansUp as we move up - radiansUpArcOrigin = radiansUpArcOriginList[n2] + radiansUpArcOrigin = radiansUpArcOriginList[n2] sinRadiansUpArcOrigin = math.sin(radiansUpArcOrigin) cosRadiansUpArcOrigin = math.cos(radiansUpArcOrigin) radiansPerElementUpArcOrigin = radiansUpArcOriginList[n2]-radiansUpArcOriginList[n2-1] - + dx_ds1 = [ orthoRadius*-sinRadiansAround*sinOrthoRadiansUp*radiansPerElementAround, orthoRadius*cosRadiansAround*sinOrthoRadiansUp*radiansPerElementAround, @@ -279,11 +270,11 @@ def generateBaseMesh(region, options): RC*sinRadiansAround*cosRadiansUpArcOrigin*radiansPerElementUpArcOrigin, -RC*sinRadiansUpArcOrigin*radiansPerElementUpArcOrigin ] - + dx_ds3 = list(interpolateCubicHermiteDerivative(v1, d1, v2, d2, xi)) dx_ds3 = vector.normalise(dx_ds3) dx_ds3 = [d*cubicArcLength/elementsCountRadial for d in dx_ds3] - + node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) @@ -291,32 +282,202 @@ def generateBaseMesh(region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) if useCrossDerivatives: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) nodeIdentifier = nodeIdentifier + 1 - - # create elements + + ################# + # Create elements + ################# + + mesh = fm.findMeshByDimension(3) + + tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + eft = tricubichermite.createEftBasic() + + tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE) + + # Regular elements + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + + # Bottom tetrahedon elements + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + # Axial elements + elementtemplate2 = mesh.createElementtemplate() + elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + # Top tetrahedron elements + elementtemplate3 = mesh.createElementtemplate() + elementtemplate3.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + # Bottom pyramid elements + elementtemplate4 = mesh.createElementtemplate() + elementtemplate4.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + # Top pyramid elements + elementtemplate5 = mesh.createElementtemplate() + elementtemplate5.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementIdentifier = 1 + no2 = elementsCountAround no3 = elementsCountAround*(elementsCountUp - 1) - rni = (1 + elementsCountUp) - no3 - no2 + 1 # regular node identifier - for e3 in range(1, elementsCountRadial): - for e2 in range(1, elementsCountUp - 1): - bni = rni + e3*no3 + e2*no2 - for e1 in range(elementsCountAround): - element = mesh.createElement(elementIdentifier, elementtemplate) - na = e1 - nb = (e1 + 1)%elementsCountAround - nodeIdentifiers = [ - bni + na, bni + nb, bni + no2 + na, bni + no2 + nb, - bni + no3 + na, bni + no3 + nb, bni + no3 + no2 + na, bni + no3 + no2 + nb + rni = (1 + elementsCountUp) - no3 - no2 + 1 # regular node identifier + radiansPerElementAround = 2.0*math.pi/elementsCountAround + for e3 in range(elementsCountRadial): + # Create elements on bottom pole + if e3 == 0: + # Create tetrahedron elements on the bottom pole + bni1 = elementsCountUp + 2 + for e1 in range(elementsCountAround): + va = e1 + vb = (e1 + 1)%elementsCountAround + eft1 = tricubichermite.createEftTetrahedronBottom(va*100, vb*100) + elementtemplate1.defineField(coordinates, -1, eft1) + element = mesh.createElement(elementIdentifier, elementtemplate1) + nodeIdentifiers = [ 1, 2, bni1 + va, bni1 + vb ] + result1 = element.setNodesByIdentifier(eft1, nodeIdentifiers) + # print('Tetrahedron bottom element', nodeIdentifiers) + # set general linear map coefficients + radiansAround = va*radiansPerElementAround + radiansAroundNext = vb*radiansPerElementAround + scalefactors = [ + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround ] - result = element.setNodesByIdentifier(eft, nodeIdentifiers) - # print('regular element', elementIdentifier, result, nodeIdentifiers) + result2 = element.setScaleFactors(eft1, scalefactors) + print('Tetrahedron Bottom element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 - + + else: + # Create pyramid elements on the bottom pole + bni4 = elementsCountUp + 1 + (e3-1)*no3 + 1 + for e1 in range(elementsCountAround): + va = e1 + vb = (e1 + 1)%elementsCountAround + # eft4 = tricubichermite.createEftPyramidBottom(va*100, vb*100) + # elementtemplate4.defineField(coordinates, -1, eft4) + # element = mesh.createElement(elementIdentifier, elementtemplate4) + nodeIdentifiers = [ 1, bni4 + va, bni4 + vb, bni4 + no3 + va, bni4 + no3 + vb ] + # result1 = element.setNodesByIdentifier(eft4, nodeIdentifiers) + #print('Pyramid bottom element', elementIdentifier, nodeIdentifiers) + # set general linear map coefficients + # radiansAround = va*radiansPerElementAround + # radiansAroundNext = vb*radiansPerElementAround + # scalefactors = [ + # -1.0, + # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + #] + #result2 = element.setScaleFactors(eft4, scalefactors) + #print('pyramid bottom element', elementIdentifier, result1, result2, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + + # create regular radial elements + for e2 in range(1, elementsCountUp-1): + if e3 == 0: + for e1 in range(elementsCountAround): + # create central radial elements: 6 node wedges + va = e1 + vb = (e1 + 1)%elementsCountAround + eft2 = tricubichermite.createEftWedgeRadial(va*100, vb*100) + elementtemplate2.defineField(coordinates, -1, eft2) + element = mesh.createElement(elementIdentifier, elementtemplate2) + bni2 = elementsCountUp + 1 + (e2-1) * no2 + 1 + nodeIdentifiers = [ e3 + e2 + 1, e3 + e2 + 2, bni2 + va, bni2 + vb, bni2 + va + elementsCountAround, bni2 + vb + elementsCountAround ] + result1 = element.setNodesByIdentifier(eft2, nodeIdentifiers) + # set general linear map coefficients + radiansAround = va*radiansPerElementAround + radiansAroundNext = vb*radiansPerElementAround + scalefactors = [ + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + ] + result2 = element.setScaleFactors(eft2, scalefactors) + print('axis element', elementIdentifier, result1, result2, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + else: + # Regular elements + bni = rni + e3*no3 + e2*no2 + for e1 in range(elementsCountAround): + element = mesh.createElement(elementIdentifier, elementtemplate) + na = e1 + nb = (e1 + 1)%elementsCountAround + nodeIdentifiers = [ + bni + na, bni + nb, bni + no2 + na, bni + no2 + nb, + bni + no3 + na, bni + no3 + nb, bni + no3 + no2 + na, bni + no3 + no2 + nb + ] + result = element.setNodesByIdentifier(eft, nodeIdentifiers) + print('regular element', elementIdentifier, result, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + + # Create elements on top pole + if e3 == 0: + # Create tetrahedron elements on the top pole + bni3 = elementsCountUp + 1 + (elementsCountUp-2) * no2 + 1 + for e1 in range(elementsCountAround): + va = e1 + vb = (e1 + 1)%elementsCountAround + eft3 = tricubichermite.createEftTetrahedronTop(va*100, vb*100) + elementtemplate3.defineField(coordinates, -1, eft3) + element = mesh.createElement(elementIdentifier, elementtemplate3) + nodeIdentifiers = [ elementsCountUp, elementsCountUp + 1, bni3 + va, bni3 + vb ] + result1 = element.setNodesByIdentifier(eft3, nodeIdentifiers) + # print('Tetrahedron top element', nodeIdentifiers) + # set general linear map coefficients + radiansAround = va*radiansPerElementAround + radiansAroundNext = vb*radiansPerElementAround + scalefactors = [ + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + ] + result2 = element.setScaleFactors(eft3, scalefactors) + print('Tetrahedron top element', elementIdentifier, result1, result2, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + + else: + # Create pyramid elements on the top pole + bni5 = elementsCountUp + 1 + (e3-1)*no3 + (elementsCountUp-2) * no2 + 1 + for e1 in range(elementsCountAround): + # va = e1 + # vb = (e1 + 1)%elementsCountAround + # eft5 = tricubichermite.createEftPyramidTop(va*100, vb*100) + # elementtemplate1.defineField(coordinates, -1, eft5) + # element = mesh.createElement(elementIdentifier, elementtemplate5) + nodeIdentifiers = [ elementsCountUp + 1, bni5 + va, bni5 + vb, bni5 + no3 + va, bni5 + no3 + vb ] + # result1 = element.setNodesByIdentifier(eft5, nodeIdentifiers) + # # print('Pyramid top element', elementIdentifier, nodeIdentifiers) + # # set general linear map coefficients + # radiansAround = va*radiansPerElementAround + # radiansAroundNext = vb*radiansPerElementAround + # scalefactors = [ + # -1.0, + # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + # ] + # result2 = element.setScaleFactors(eft5, scalefactors) + # print('pyramid top element', elementIdentifier, result1, result2, nodeIdentifiers) + elementIdentifier = elementIdentifier + 1 + fm.endChange() @staticmethod diff --git a/scaffoldmaker/scaffoldmaker.py b/scaffoldmaker/scaffoldmaker.py index e549b8cf..0a27faad 100644 --- a/scaffoldmaker/scaffoldmaker.py +++ b/scaffoldmaker/scaffoldmaker.py @@ -21,6 +21,8 @@ from scaffoldmaker.meshtypes.meshtype_3d_tube1 import MeshType_3d_tube1 from scaffoldmaker.meshtypes.meshtype_3d_tubeseptum1 import MeshType_3d_tubeseptum1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 +from scaffoldmaker.meshtypes.meshtype_3d_truncatedsphere1 import MeshType_3d_truncatedsphere1 +from scaffoldmaker.meshtypes.meshtype_3d_lens1 import MeshType_3d_lens1 class Scaffoldmaker(object): @@ -44,7 +46,9 @@ def __init__(self): MeshType_3d_sphereshellseptum1, MeshType_3d_tube1, MeshType_3d_tubeseptum1, - MeshType_3d_solidsphere1 + MeshType_3d_solidsphere1, + MeshType_3d_truncatedsphere1, + MeshType_3d_lens1 ] def getMeshTypes(self): diff --git a/scaffoldmaker/utils/eft_utils.py b/scaffoldmaker/utils/eft_utils.py index d493419a..faa8684a 100644 --- a/scaffoldmaker/utils/eft_utils.py +++ b/scaffoldmaker/utils/eft_utils.py @@ -98,8 +98,8 @@ def remapEftNodeValueLabel(eft, localNodeIndexes, fromValueLabel, expressionTerm :param expressionTerms: List of (valueLabel, scaleFactorIndexesList ) to remap to. e.g. [ (Node.VALUE_LABEL_D_DS2, []), (Node.VALUE_LABEL_D_DS3, [5, 6]) ] ''' - functionCount = eft.getNumberOfFunctions() - for f in range(1, functionCount + 1): + functionCount = eft.getNumberOfFunctions() + for f in range(1, functionCount + 1): if eft.getFunctionNumberOfTerms(f) == 1: localNodeIndex = eft.getTermLocalNodeIndex(f, 1) if (localNodeIndex in localNodeIndexes) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and (not getEftTermScaling(eft, f, 1)): diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index ba2ef0bc..02bd3ce5 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -173,6 +173,197 @@ def createEftShellApexTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): assert eft.validate(), 'eftfactory_tricubichermite.createEftShellApexTop: Failed to validate eft' return eft + def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for the central axis of a solid + cylinder, with xi1 collapsed, xi2 up and xi3 out radially. + Each collapsed node has 3 scale factors giving the cos, sin coefficients + of the radial line from global derivatives, plus the arc subtended by + the element in radians, so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic to remap D2_DS1DS2 at apex + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + if not self._useCrossDerivatives: + for n in [ 4, 5, 6, 7 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) + eft.setFunctionNumberOfTerms(n*8 + 6, 0) + eft.setFunctionNumberOfTerms(n*8 + 7, 0) + eft.setFunctionNumberOfTerms(n*8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + # remap parameters before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) + for layer in range(2): + so = layer*6 + 1 + ln = layer*2 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer*2 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 1, 2, 2, 3, 4, 5, 6 ] + remapEftLocalNodes(eft, 6, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeRadial: Failed to validate eft' + return eft + + def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for a solid tetrahedron for the axis of a + solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and + xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients + of the radial line from global derivatives, plus the arc subtended by + the element in radians, so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic to remap D2_DS1DS2 at apex + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 4, 5, 6, 7 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 + if not self._useCrossDerivatives: + eft.setFunctionNumberOfTerms(n*8 + 6, 0) # D2_DS1DS3 + eft.setFunctionNumberOfTerms(n*8 + 7, 0) # D2_DS2DS3 + eft.setFunctionNumberOfTerms(n*8 + 8, 0) # D3_DS1DS2DS3 + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + + # remap parameters on xi3 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) + # reverse D_DS2 at bottom + remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + for layer in range(2): + so = layer*6 + 1 + ln = layer*2 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer*2 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + # remap parameters on xi3 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) + remapEftNodeValueLabel(eft, [ 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [])]) + + ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' + return eft + + def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + ''' + Create a tricubic hermite element field for a solid tetrahedron for the axis of a + solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and + xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients + of the radial line from global derivatives, plus the arc subtended by + the element in radians, so the circumferential direction is rounded. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :return: Element field template + ''' + # start with full tricubic to remap D2_DS1DS2 at apex + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 4, 5, 6, 7 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 + if not self._useCrossDerivatives: + eft.setFunctionNumberOfTerms(n*8 + 6, 0) # D2_DS1DS3 + eft.setFunctionNumberOfTerms(n*8 + 7, 0) # D2_DS2DS3 + eft.setFunctionNumberOfTerms(n*8 + 8, 0) # D3_DS1DS2DS3 + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + + # remap parameters on xi3 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) + # reverse D_DS2 at bottom + # remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + for layer in range(2): + so = layer*6 + 1 + ln = layer*2 + 1 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer*2 + 2 + # 2 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + # remap parameters on xi3 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) + remapEftNodeValueLabel(eft, [ 7, 8 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1])]) + + ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] + remapEftLocalNodes(eft, 4, ln_map) + + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' + return eft + def createEftSplitXi1LeftStraight(self): ''' Create an element field template suitable for the inner elements of the From e7af20b918a7c9e8b7dfb1dbaa89a5b9c32192a2 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 7 Aug 2018 14:33:01 +1200 Subject: [PATCH 3/9] Complete sphere with pyramid elements for multiple radial layers --- .../meshtypes/meshtype_3d_solidsphere1.py | 113 +++++---- .../utils/eftfactory_tricubichermite.py | 237 ++++++++++++++---- 2 files changed, 256 insertions(+), 94 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 6ba84df7..ce8c68af 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -175,7 +175,7 @@ def generateBaseMesh(region, options): cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, radius]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) @@ -330,58 +330,62 @@ def generateBaseMesh(region, options): no3 = elementsCountAround*(elementsCountUp - 1) rni = (1 + elementsCountUp) - no3 - no2 + 1 # regular node identifier radiansPerElementAround = 2.0*math.pi/elementsCountAround + for e3 in range(elementsCountRadial): - # Create elements on bottom pole + # Create elements on bottom pole if e3 == 0: # Create tetrahedron elements on the bottom pole - bni1 = elementsCountUp + 2 + bni1 = elementsCountUp + 2 + radiansInclineNext = math.pi*0.5/elementsCountRadial + for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1)%elementsCountAround - eft1 = tricubichermite.createEftTetrahedronBottom(va*100, vb*100) + eft1 = tricubichermite.createEftTetrahedronBottom(va*100, vb*100, 10000 + 2) elementtemplate1.defineField(coordinates, -1, eft1) element = mesh.createElement(elementIdentifier, elementtemplate1) nodeIdentifiers = [ 1, 2, bni1 + va, bni1 + vb ] result1 = element.setNodesByIdentifier(eft1, nodeIdentifiers) - # print('Tetrahedron bottom element', nodeIdentifiers) # set general linear map coefficients radiansAround = va*radiansPerElementAround radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround ] result2 = element.setScaleFactors(eft1, scalefactors) - print('Tetrahedron Bottom element', elementIdentifier, result1, result2, nodeIdentifiers) + # print('Tetrahedron Bottom element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 else: # Create pyramid elements on the bottom pole bni4 = elementsCountUp + 1 + (e3-1)*no3 + 1 + radiansIncline = math.pi*0.5*e3/elementsCountRadial + radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial + for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1)%elementsCountAround - # eft4 = tricubichermite.createEftPyramidBottom(va*100, vb*100) - # elementtemplate4.defineField(coordinates, -1, eft4) - # element = mesh.createElement(elementIdentifier, elementtemplate4) + eft4 = tricubichermite.createEftPyramidBottom(va*100, vb*100, 100000 + e3*2) + elementtemplate4.defineField(coordinates, -1, eft4) + element = mesh.createElement(elementIdentifier, elementtemplate4) nodeIdentifiers = [ 1, bni4 + va, bni4 + vb, bni4 + no3 + va, bni4 + no3 + vb ] - # result1 = element.setNodesByIdentifier(eft4, nodeIdentifiers) - #print('Pyramid bottom element', elementIdentifier, nodeIdentifiers) + result1 = element.setNodesByIdentifier(eft4, nodeIdentifiers) # set general linear map coefficients - # radiansAround = va*radiansPerElementAround - # radiansAroundNext = vb*radiansPerElementAround - # scalefactors = [ - # -1.0, - # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, - # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround - #] - #result2 = element.setScaleFactors(eft4, scalefactors) - #print('pyramid bottom element', elementIdentifier, result1, result2, nodeIdentifiers) + radiansAround = va*radiansPerElementAround + radiansAroundNext = vb*radiansPerElementAround + scalefactors = [ + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) + ] + result2 = element.setScaleFactors(eft4, scalefactors) + # print('pyramid bottom element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 # create regular radial elements @@ -408,11 +412,12 @@ def generateBaseMesh(region, options): math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround ] result2 = element.setScaleFactors(eft2, scalefactors) - print('axis element', elementIdentifier, result1, result2, nodeIdentifiers) + # print('axis element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 else: # Regular elements bni = rni + e3*no3 + e2*no2 + for e1 in range(elementsCountAround): element = mesh.createElement(elementIdentifier, elementtemplate) na = e1 @@ -422,22 +427,23 @@ def generateBaseMesh(region, options): bni + no3 + na, bni + no3 + nb, bni + no3 + no2 + na, bni + no3 + no2 + nb ] result = element.setNodesByIdentifier(eft, nodeIdentifiers) - print('regular element', elementIdentifier, result, nodeIdentifiers) + # print('regular element', elementIdentifier, result, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 # Create elements on top pole if e3 == 0: - # Create tetrahedron elements on the top pole + # # Create tetrahedron elements on the top pole bni3 = elementsCountUp + 1 + (elementsCountUp-2) * no2 + 1 + radiansInclineNext = math.pi*0.5/elementsCountRadial + for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1)%elementsCountAround - eft3 = tricubichermite.createEftTetrahedronTop(va*100, vb*100) + eft3 = tricubichermite.createEftTetrahedronTop(va*100, vb*100, 100000 + 2) elementtemplate3.defineField(coordinates, -1, eft3) element = mesh.createElement(elementIdentifier, elementtemplate3) nodeIdentifiers = [ elementsCountUp, elementsCountUp + 1, bni3 + va, bni3 + vb ] result1 = element.setNodesByIdentifier(eft3, nodeIdentifiers) - # print('Tetrahedron top element', nodeIdentifiers) # set general linear map coefficients radiansAround = va*radiansPerElementAround radiansAroundNext = vb*radiansPerElementAround @@ -445,36 +451,39 @@ def generateBaseMesh(region, options): -1.0, math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) ] result2 = element.setScaleFactors(eft3, scalefactors) - print('Tetrahedron top element', elementIdentifier, result1, result2, nodeIdentifiers) + # print('Tetrahedron top element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 else: # Create pyramid elements on the top pole bni5 = elementsCountUp + 1 + (e3-1)*no3 + (elementsCountUp-2) * no2 + 1 + radiansIncline = math.pi*0.5*e3/elementsCountRadial + radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial + for e1 in range(elementsCountAround): - # va = e1 - # vb = (e1 + 1)%elementsCountAround - # eft5 = tricubichermite.createEftPyramidTop(va*100, vb*100) - # elementtemplate1.defineField(coordinates, -1, eft5) - # element = mesh.createElement(elementIdentifier, elementtemplate5) - nodeIdentifiers = [ elementsCountUp + 1, bni5 + va, bni5 + vb, bni5 + no3 + va, bni5 + no3 + vb ] - # result1 = element.setNodesByIdentifier(eft5, nodeIdentifiers) - # # print('Pyramid top element', elementIdentifier, nodeIdentifiers) - # # set general linear map coefficients - # radiansAround = va*radiansPerElementAround - # radiansAroundNext = vb*radiansPerElementAround - # scalefactors = [ - # -1.0, - # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, - # math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - # math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround - # ] - # result2 = element.setScaleFactors(eft5, scalefactors) + va = e1 + vb = (e1 + 1)%elementsCountAround + eft5 = tricubichermite.createEftPyramidTop(va*100, vb*100, 100000 + e3*2) + + elementtemplate5.defineField(coordinates, -1, eft5) + element = mesh.createElement(elementIdentifier, elementtemplate5) + nodeIdentifiers = [ bni5 + va, bni5 + vb, elementsCountUp + 1, bni5 + no3 + va, bni5 + no3 + vb ] + result1 = element.setNodesByIdentifier(eft5, nodeIdentifiers) + # set general linear map coefficients + radiansAround = va*radiansPerElementAround + radiansAroundNext = vb*radiansPerElementAround + scalefactors = [ + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) + ] + result2 = element.setScaleFactors(eft5, scalefactors) # print('pyramid top element', elementIdentifier, result1, result2, nodeIdentifiers) elementIdentifier = elementIdentifier + 1 @@ -499,4 +508,4 @@ def generateMesh(region, options): MeshType_3d_solidsphere1.generateBaseMesh(baseRegion, options) meshrefinement = MeshRefinement(baseRegion, region) - meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountUp, refineElementsCountRadial) + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountUp, refineElementsCountRadial) \ No newline at end of file diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index 02bd3ce5..2cacadb8 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -187,7 +187,7 @@ def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 :return: Element field template ''' - # start with full tricubic to remap D2_DS1DS2 at apex + # start with full tricubic eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) if not self._useCrossDerivatives: for n in [ 4, 5, 6, 7 ]: @@ -202,12 +202,13 @@ def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + # remap parameters before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) for layer in range(2): so = layer*6 + 1 ln = layer*2 + 1 - # 2 terms for d/dxi2 via general linear map: + # 2 terms for d/dxi3 via general linear map: remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) @@ -217,7 +218,7 @@ def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln = layer*2 + 2 - # 2 terms for d/dxi2 via general linear map: + # 2 terms for d/dxi3 via general linear map: remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) @@ -232,21 +233,23 @@ def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeRadial: Failed to validate eft' return eft - def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients of the radial line from global derivatives, plus the arc subtended by - the element in radians, so the circumferential direction is rounded. - Need to create a new template for each sector around axis giving common - nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and - add 100 for each radial line around axis. + the element in radians, so the circumferential direction is rounded. Collapsed + node on xi2 = 0 has additional 2 scale factors giving the cos, sin coefficients of the + angle of inclination from global z-axis. Need to create a new template for each sector + around axis giving common nodeScaleFactorOffset values on common faces. + Suggestion is to start at 0 and add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases :return: Element field template ''' - # start with full tricubic to remap D2_DS1DS2 at apex + # start with full tricubic eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) for n in [ 4, 5, 6, 7 ]: eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 @@ -257,32 +260,47 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) - + # remap parameters on xi3 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) # reverse D_DS2 at bottom remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + for layer in range(2): - so = layer*6 + 1 + so = layer*10 + 1 ln = layer*2 + 1 - # 2 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + if layer == 0: + # 3 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) + else: + # 2 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln = layer*2 + 2 - # 2 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + if layer == 0: + # 3 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) + else: + # 2 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -290,7 +308,7 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [ 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [ 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, []) ]) ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] remapEftLocalNodes(eft, 4, ln_map) @@ -298,21 +316,23 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' return eft - def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): + def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a - solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and + solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients - of the radial line from global derivatives, plus the arc subtended by - the element in radians, so the circumferential direction is rounded. - Need to create a new template for each sector around axis giving common + of the radial line from global derivatives, plus the arc subtended by the element in + radians, so the circumferential direction is rounded. Collapsed node on xi2 = 1 has + additional 2 scale factors giving the cos, sin coefficients of the angle of inclination + from global z-axis. Need to create a new template for each sector around axis giving common nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases :return: Element field template ''' - # start with full tricubic to remap D2_DS1DS2 at apex + # start with full tricubic eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) for n in [ 4, 5, 6, 7 ]: eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 @@ -325,30 +345,41 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 setEftScaleFactorIds(eft, [1], [ nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) - + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2 ]) + # remap parameters on xi3 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) - # reverse D_DS2 at bottom - # remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + for layer in range(2): so = layer*6 + 1 ln = layer*2 + 1 - # 2 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + if layer == 0: + # 2 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + else: + # 3 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [ so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [ 1, so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln = layer*2 + 2 - # 2 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + if layer == 0: + # 2 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + else: + # 3 terms for d/dxi3 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 6, so + 10]), (Node.VALUE_LABEL_D_DS3, [ so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [ 1, so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -356,7 +387,7 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [ 7, 8 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, [ 7, 8 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] remapEftLocalNodes(eft, 4, ln_map) @@ -364,6 +395,128 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' return eft + def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): + ''' + Create a tricubic hermite element field for a solid pyramid for the axis of a + solid sphere pole, with xi1 and xi3 collapsed on xi2 = 0. Each collapsed node + has 5 scale factors giving the cos, sin coefficients of the radial line from + global derivatives, plus the arc subtended by the element in radians, so the + circumferential direction is rounded, cos and sin coefficients of the inclination + angle. Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :return: Element field template + ''' + # start with full tricubic + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) + eft.setFunctionNumberOfTerms(n*8 + 6, 0) + eft.setFunctionNumberOfTerms(n*8 + 7, 0) + eft.setFunctionNumberOfTerms(n*8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4 ]) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) + + for layer in range(2): + so = layer*10 + 1 + ln = layer*4 + 1 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer*4 + 2 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS3, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ] + remapEftLocalNodes(eft, 5, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' + return eft + + def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): + ''' + Create a tricubic hermite element field for a solid pyramid for the axis of a + solid sphere pole, with xi1 and xi3 collapsed on xi2 = 1. Each collapsed node + has 5 scale factors giving the cos, sin coefficients of the radial line from + global derivatives, plus the arc subtended by the element in radians, so the + circumferential direction is rounded, cos and sin coefficients of the inclination + angle. Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. + :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 + :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 + :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :return: Element field template + ''' + # start with full tricubic + eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) + for n in [ 0, 1, 4, 5 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) + eft.setFunctionNumberOfTerms(n*8 + 6, 0) + eft.setFunctionNumberOfTerms(n*8 + 7, 0) + eft.setFunctionNumberOfTerms(n*8 + 8, 0) + + # GRC: allow scale factor identifier for global -1.0 to be prescribed + setEftScaleFactorIds(eft, [1], [ + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4 ]) + + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS3, []) + + for layer in range(2): + so = layer*10 + 1 + ln = layer*4 + 3 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 1, so + 3, so + 5]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln = layer*4 + 4 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS3, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 6, so + 8, so + 10]) ]) + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [ 1, 2, 3, 3, 4, 5, 3, 3 ] + remapEftLocalNodes(eft, 5, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' + return eft + def createEftSplitXi1LeftStraight(self): ''' Create an element field template suitable for the inner elements of the From e25653ffc9fdcde1482737e3a9cc31215f323e2d Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 7 Aug 2018 15:22:58 +1200 Subject: [PATCH 4/9] Switch D_DS2 and D_DS3 at poles --- .../meshtypes/meshtype_3d_solidsphere1.py | 9 ++--- .../utils/eftfactory_tricubichermite.py | 40 ++++++++++--------- 2 files changed, 26 insertions(+), 23 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index ce8c68af..20663160 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -145,8 +145,8 @@ def generateBaseMesh(region, options): cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, -radius*2/elementsCountUp]) # height of element along central axis - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, -radius*2/elementsCountUp]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) @@ -160,7 +160,6 @@ def generateBaseMesh(region, options): cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius+n2*2*radius/elementsCountUp]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ cubicArcLengthList[n2]/elementsCountRadial, 0.0, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, cubicArcLengthList[n2]/elementsCountRadial, 0.0 ]) if useCrossDerivatives: @@ -175,8 +174,8 @@ def generateBaseMesh(region, options): cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, radius]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, radius*2/elementsCountUp ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index 2cacadb8..d2c7fd05 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -267,17 +267,18 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs # remap parameters on xi3 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) - # reverse D_DS2 at bottom - remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + # switch D_DS2 and D_DS and negate at bottom + remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, [1]) ]) for layer in range(2): so = layer*10 + 1 ln = layer*2 + 1 if layer == 0: # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) + else: # 2 terms for d/dxi3 via general linear map: remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) @@ -292,9 +293,10 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs ln = layer*2 + 2 if layer == 0: # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) + else: # 2 terms for d/dxi3 via general linear map: remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) @@ -361,9 +363,10 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) else: # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [ so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [ 1, so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [ so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [ 1, so + 4]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, []) ]) # switch D_DS2 and D_DS3 at top # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -377,9 +380,10 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) else: # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 6, so + 10]), (Node.VALUE_LABEL_D_DS3, [ so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [ 1, so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 6, so + 10]), (Node.VALUE_LABEL_D_DS2, [ so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [ 1, so + 9]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, []) ]) # switch D_DS2 and D_DS3 at top # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -433,9 +437,9 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, so = layer*10 + 1 ln = layer*4 + 1 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -443,9 +447,9 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, ln = layer*4 + 2 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS3, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -494,9 +498,9 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no so = layer*10 + 1 ln = layer*4 + 3 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 4]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 1, so + 3, so + 5]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -504,9 +508,9 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no ln = layer*4 + 4 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS3, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 9]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 6, so + 8, so + 10]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) From a511d1c5b977f272144386c195be5b8e48ae309c Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 7 Aug 2018 16:47:29 +1200 Subject: [PATCH 5/9] Remove links to truncated sphere, lens --- scaffoldmaker/scaffoldmaker.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/scaffoldmaker/scaffoldmaker.py b/scaffoldmaker/scaffoldmaker.py index 0a27faad..ff7234c4 100644 --- a/scaffoldmaker/scaffoldmaker.py +++ b/scaffoldmaker/scaffoldmaker.py @@ -21,8 +21,6 @@ from scaffoldmaker.meshtypes.meshtype_3d_tube1 import MeshType_3d_tube1 from scaffoldmaker.meshtypes.meshtype_3d_tubeseptum1 import MeshType_3d_tubeseptum1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 -from scaffoldmaker.meshtypes.meshtype_3d_truncatedsphere1 import MeshType_3d_truncatedsphere1 -from scaffoldmaker.meshtypes.meshtype_3d_lens1 import MeshType_3d_lens1 class Scaffoldmaker(object): @@ -46,9 +44,7 @@ def __init__(self): MeshType_3d_sphereshellseptum1, MeshType_3d_tube1, MeshType_3d_tubeseptum1, - MeshType_3d_solidsphere1, - MeshType_3d_truncatedsphere1, - MeshType_3d_lens1 + MeshType_3d_solidsphere1 ] def getMeshTypes(self): From 69b827af91be68cc438bb5204184209b7c9d4f58 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Fri, 10 Aug 2018 16:03:11 +1200 Subject: [PATCH 6/9] Make changes from review --- .../meshtypes/meshtype_3d_solidsphere1.py | 62 +++++++++---------- scaffoldmaker/scaffoldmaker.py | 6 +- scaffoldmaker/utils/eft_utils.py | 4 +- .../utils/eftfactory_tricubichermite.py | 12 ++-- 4 files changed, 44 insertions(+), 40 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 20663160..71b3037d 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -27,7 +27,7 @@ def getDefaultOptions(): 'Number of elements around' : 8, 'Number of elements up' : 8, 'Number of elements radial' : 1, - 'Diameter' : 1, + 'Diameter' : 1.0, 'Use cross derivatives' : False, 'Refine' : False, 'Refine number of elements around' : 1, @@ -41,7 +41,7 @@ def getOrderedOptionNames(): 'Number of elements around', 'Number of elements up', 'Number of elements radial', - 'Diameter', + 'Diameter', 'Use cross derivatives', 'Refine', 'Refine number of elements around', @@ -62,11 +62,11 @@ def checkOptions(options): options['Number of elements up'] = 4 if options['Number of elements around'] < 4: options['Number of elements around'] = 4 - if options['Diameter'] < 0: - options['Diameter'] = 1 + if options['Diameter'] < 0.0: + options['Diameter'] = 1.0 - @staticmethod - def generateBaseMesh(region, options): + @classmethod + def generateBaseMesh(cls, region, options): """ Generate the base tricubic Hermite mesh. See also generateMesh(). :param region: Zinc region to define model in. Must be empty. @@ -78,7 +78,7 @@ def generateBaseMesh(region, options): elementsCountRadial = options['Number of elements radial'] useCrossDerivatives = options['Use cross derivatives'] diameter = options['Diameter'] - radius = diameter/2 + radius = diameter/2.0 fm = region.getFieldmodule() fm.beginChange() @@ -121,7 +121,7 @@ def generateBaseMesh(region, options): dx_ds3 = [ 0.0, 0.0, 0.0 ] zero = [ 0.0, 0.0, 0.0 ] - cubicArcLengthList = [0]*(elementsCountUp+1) + cubicArcLengthList = [0.0]*(elementsCountUp+1) # Pre-calculate cubicArcLength along elementsCountUp for n2 in range(1,elementsCountUp + 1): @@ -130,14 +130,14 @@ def generateBaseMesh(region, options): sinRadiansUp = math.sin(radiansUp) # Calculate cubic hermite arclength linking point on axis to surface on sphere - v1 = [0,0,-radius+n2*2*radius/elementsCountUp] - d1 = [0,1,0] + v1 = [0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp] + d1 = [0.0, 1.0, 0.0] v2 = [ - radius*math.cos(math.pi/2)*sinRadiansUp, - radius*math.sin(math.pi/2)*sinRadiansUp, + radius*math.cos(math.pi/2.0)*sinRadiansUp, + radius*math.sin(math.pi/2.0)*sinRadiansUp, -radius*cosRadiansUp ] - d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] + d2 = [math.cos(math.pi/2.0)*sinRadiansUp,math.sin(math.pi/2.0)*sinRadiansUp,-cosRadiansUp] cubicArcLengthList[n2] = computeCubicHermiteArcLength(v1, d1, v2, d2, True) # Create node for bottom pole @@ -146,7 +146,7 @@ def generateBaseMesh(region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, -radius*2/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, -radius*2.0/elementsCountUp]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) @@ -158,9 +158,9 @@ def generateBaseMesh(region, options): for n2 in range(1,elementsCountUp): node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius+n2*2*radius/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ cubicArcLengthList[n2]/elementsCountRadial, 0.0, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2/elementsCountUp]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, 0.0, radius*2.0/elementsCountUp]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, cubicArcLengthList[n2]/elementsCountRadial, 0.0 ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) @@ -175,7 +175,7 @@ def generateBaseMesh(region, options): coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [ 0.0, 0.0, radius]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, [ radius*radiansPerElementUp, 0.0, 0.0 ]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, [ 0.0, radius*radiansPerElementUp, 0.0 ]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, radius*2/elementsCountUp ]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, [ 0.0, 0.0, radius*2.0/elementsCountUp ]) if useCrossDerivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) @@ -186,11 +186,11 @@ def generateBaseMesh(region, options): # Create other nodes for n3 in range(1,elementsCountRadial+1): xi = 1/elementsCountRadial*n3 - radiansUpArcOriginList = [0]*(elementsCountUp) + radiansUpArcOriginList = [0.0]*(elementsCountUp) # Pre-calculate RC for points on vertical arc running between top and bottom poles - pt = [0, radius*xi, 0] - arcOrigin = (radius*radius - pt[2]*pt[2] - pt[1]*pt[1])/(-2*pt[1]) + pt = [0.0, radius*xi, 0.0] + arcOrigin = (radius*radius - pt[2]*pt[2] - pt[1]*pt[1])/(-2.0*pt[1]) RC = math.sqrt(arcOrigin*arcOrigin + radius*radius) radiansUpArcOriginList[0] = math.acos(-radius/RC) @@ -203,16 +203,16 @@ def generateBaseMesh(region, options): # Calculate node coordinates on arc using cubic hermite interpolation cubicArcLength = cubicArcLengthList[n2] - v1 = [0,0,-radius+n2*2*radius/elementsCountUp] - d1 = [math.cos(math.pi/2),math.sin(math.pi/2),0] + v1 = [0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp] + d1 = [math.cos(math.pi/2.0), math.sin(math.pi/2.0), 0.0] d1 = vector.normalise(d1) d1 = [d*cubicArcLength for d in d1] v2 = [ - radius*math.cos(math.pi/2)*sinRadiansUp, - radius*math.sin(math.pi/2)*sinRadiansUp, + radius*math.cos(math.pi/2.0)*sinRadiansUp, + radius*math.sin(math.pi/2.0)*sinRadiansUp, -radius*cosRadiansUp ] - d2 = [math.cos(math.pi/2)*sinRadiansUp,math.sin(math.pi/2)*sinRadiansUp,-cosRadiansUp] + d2 = [math.cos(math.pi/2.0)*sinRadiansUp,math.sin(math.pi/2.0)*sinRadiansUp,-cosRadiansUp] d2 = vector.normalise(d2) d2 = [d*cubicArcLength for d in d2] x = list(interpolateCubicHermite(v1, d1, v2, d2, xi)) @@ -232,8 +232,8 @@ def generateBaseMesh(region, options): cubicArcLength = cubicArcLengthList[n2] # Calculate node coordinates on arc using cubic hermite interpolation - v1 = [0,0,-radius+n2*2*radius/elementsCountUp] - d1 = [cosRadiansAround,sinRadiansAround,0] + v1 = [0.0, 0.0, -radius+n2*2.0*radius/elementsCountUp] + d1 = [cosRadiansAround, sinRadiansAround, 0.0] d1 = vector.normalise(d1) d1 = [d*cubicArcLength for d in d1] v2 = [ @@ -488,15 +488,15 @@ def generateBaseMesh(region, options): fm.endChange() - @staticmethod - def generateMesh(region, options): + @classmethod + def generateMesh(cls, region, options): """ Generate base or refined mesh. :param region: Zinc region to create mesh in. Must be empty. :param options: Dict containing options. See getDefaultOptions(). """ if not options['Refine']: - MeshType_3d_solidsphere1.generateBaseMesh(region, options) + cls.generateBaseMesh(region, options) return refineElementsCountAround = options['Refine number of elements around'] @@ -504,7 +504,7 @@ def generateMesh(region, options): refineElementsCountRadial = options['Refine number of elements radial'] baseRegion = region.createRegion() - MeshType_3d_solidsphere1.generateBaseMesh(baseRegion, options) + cls.generateBaseMesh(baseRegion, options) meshrefinement = MeshRefinement(baseRegion, region) meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountUp, refineElementsCountRadial) \ No newline at end of file diff --git a/scaffoldmaker/scaffoldmaker.py b/scaffoldmaker/scaffoldmaker.py index ff7234c4..6c6a6ff3 100644 --- a/scaffoldmaker/scaffoldmaker.py +++ b/scaffoldmaker/scaffoldmaker.py @@ -16,11 +16,11 @@ from scaffoldmaker.meshtypes.meshtype_3d_heartventricles2 import MeshType_3d_heartventricles2 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase1 import MeshType_3d_heartventriclesbase1 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase2 import MeshType_3d_heartventriclesbase2 +from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 from scaffoldmaker.meshtypes.meshtype_3d_sphereshell1 import MeshType_3d_sphereshell1 from scaffoldmaker.meshtypes.meshtype_3d_sphereshellseptum1 import MeshType_3d_sphereshellseptum1 from scaffoldmaker.meshtypes.meshtype_3d_tube1 import MeshType_3d_tube1 from scaffoldmaker.meshtypes.meshtype_3d_tubeseptum1 import MeshType_3d_tubeseptum1 -from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 class Scaffoldmaker(object): @@ -40,11 +40,11 @@ def __init__(self): MeshType_3d_heartventricles2, MeshType_3d_heartventriclesbase1, MeshType_3d_heartventriclesbase2, + MeshType_3d_solidsphere1, MeshType_3d_sphereshell1, MeshType_3d_sphereshellseptum1, MeshType_3d_tube1, - MeshType_3d_tubeseptum1, - MeshType_3d_solidsphere1 + MeshType_3d_tubeseptum1 ] def getMeshTypes(self): diff --git a/scaffoldmaker/utils/eft_utils.py b/scaffoldmaker/utils/eft_utils.py index faa8684a..d493419a 100644 --- a/scaffoldmaker/utils/eft_utils.py +++ b/scaffoldmaker/utils/eft_utils.py @@ -98,8 +98,8 @@ def remapEftNodeValueLabel(eft, localNodeIndexes, fromValueLabel, expressionTerm :param expressionTerms: List of (valueLabel, scaleFactorIndexesList ) to remap to. e.g. [ (Node.VALUE_LABEL_D_DS2, []), (Node.VALUE_LABEL_D_DS3, [5, 6]) ] ''' - functionCount = eft.getNumberOfFunctions() - for f in range(1, functionCount + 1): + functionCount = eft.getNumberOfFunctions() + for f in range(1, functionCount + 1): if eft.getFunctionNumberOfTerms(f) == 1: localNodeIndex = eft.getTermLocalNodeIndex(f, 1) if (localNodeIndex in localNodeIndexes) and (eft.getTermNodeValueLabel(f, 1) == fromValueLabel) and (not getEftTermScaling(eft, f, 1)): diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index d2c7fd05..1babfe05 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -246,7 +246,8 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs Suggestion is to start at 0 and add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 - :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :param nodeScaleFactorOffsetUp: offset of first scale factor for inclination at pole, + increase by 2 in each layer away from axis. Suggest starting at 100000 on axis. :return: Element field template ''' # start with full tricubic @@ -331,7 +332,8 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 - :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :param nodeScaleFactorOffsetUp: offset of first scale factor for inclination at pole, + increase by 2 in each layer away from axis. Suggest starting at 100000 on axis. :return: Element field template ''' # start with full tricubic @@ -411,7 +413,8 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 - :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :param nodeScaleFactorOffsetUp: offset of first scale factor for inclination at pole, + increase by 2 in each layer away from axis. Suggest starting at 100000 on axis. :return: Element field template ''' # start with full tricubic @@ -472,7 +475,8 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 - :param nodeScaleFactorOffsetUp: offset of node scale factors as inclination increases + :param nodeScaleFactorOffsetUp: offset of first scale factor for inclination at pole, + increase by 2 in each layer away from axis. Suggest starting at 100000 on axis. :return: Element field template ''' # start with full tricubic From a629ba782689ae2119e2d39a9973532077f840b6 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Aug 2018 14:15:17 +1200 Subject: [PATCH 7/9] Modify tetrahedron elements to conform with pyramid elements --- .../meshtypes/meshtype_3d_solidsphere1.py | 31 +-- .../utils/eftfactory_tricubichermite.py | 191 +++++++++--------- 2 files changed, 108 insertions(+), 114 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 71b3037d..110c1e57 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -331,16 +331,19 @@ def generateBaseMesh(cls, region, options): radiansPerElementAround = 2.0*math.pi/elementsCountAround for e3 in range(elementsCountRadial): + # Create elements on bottom pole + radiansIncline = math.pi*0.5*e3/elementsCountRadial + radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial + if e3 == 0: # Create tetrahedron elements on the bottom pole bni1 = elementsCountUp + 2 - radiansInclineNext = math.pi*0.5/elementsCountRadial for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1)%elementsCountAround - eft1 = tricubichermite.createEftTetrahedronBottom(va*100, vb*100, 10000 + 2) + eft1 = tricubichermite.createEftTetrahedronBottom(va*100, vb*100, 10000) elementtemplate1.defineField(coordinates, -1, eft1) element = mesh.createElement(elementIdentifier, elementtemplate1) nodeIdentifiers = [ 1, 2, bni1 + va, bni1 + vb ] @@ -350,10 +353,12 @@ def generateBaseMesh(cls, region, options): radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround),radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround ] result2 = element.setScaleFactors(eft1, scalefactors) # print('Tetrahedron Bottom element', elementIdentifier, result1, result2, nodeIdentifiers) @@ -362,8 +367,6 @@ def generateBaseMesh(cls, region, options): else: # Create pyramid elements on the bottom pole bni4 = elementsCountUp + 1 + (e3-1)*no3 + 1 - radiansIncline = math.pi*0.5*e3/elementsCountRadial - radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial for e1 in range(elementsCountAround): va = e1 @@ -430,15 +433,17 @@ def generateBaseMesh(cls, region, options): elementIdentifier = elementIdentifier + 1 # Create elements on top pole + radiansIncline = math.pi*0.5*e3/elementsCountRadial + radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial + if e3 == 0: # # Create tetrahedron elements on the top pole bni3 = elementsCountUp + 1 + (elementsCountUp-2) * no2 + 1 - radiansInclineNext = math.pi*0.5/elementsCountRadial for e1 in range(elementsCountAround): va = e1 vb = (e1 + 1)%elementsCountAround - eft3 = tricubichermite.createEftTetrahedronTop(va*100, vb*100, 100000 + 2) + eft3 = tricubichermite.createEftTetrahedronTop(va*100, vb*100, 100000) elementtemplate3.defineField(coordinates, -1, eft3) element = mesh.createElement(elementIdentifier, elementtemplate3) nodeIdentifiers = [ elementsCountUp, elementsCountUp + 1, bni3 + va, bni3 + vb ] @@ -448,10 +453,12 @@ def generateBaseMesh(cls, region, options): radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), + math.cos(radiansAround), math.sin(radiansAround),radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround ] result2 = element.setScaleFactors(eft3, scalefactors) # print('Tetrahedron top element', elementIdentifier, result1, result2, nodeIdentifiers) @@ -460,8 +467,6 @@ def generateBaseMesh(cls, region, options): else: # Create pyramid elements on the top pole bni5 = elementsCountUp + 1 + (e3-1)*no3 + (elementsCountUp-2) * no2 + 1 - radiansIncline = math.pi*0.5*e3/elementsCountRadial - radiansInclineNext = math.pi*0.5*(e3 + 1)/elementsCountRadial for e1 in range(elementsCountAround): va = e1 diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index 1babfe05..abc68c11 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -236,14 +236,14 @@ def createEftWedgeRadial(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1): def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a - solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and - xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients - of the radial line from global derivatives, plus the arc subtended by - the element in radians, so the circumferential direction is rounded. Collapsed - node on xi2 = 0 has additional 2 scale factors giving the cos, sin coefficients of the - angle of inclination from global z-axis. Need to create a new template for each sector - around axis giving common nodeScaleFactorOffset values on common faces. - Suggestion is to start at 0 and add 100 for each radial line around axis. + solid sphere pole, with xi1 and xi3 collapsed on xi2 = 0, and xi1 collapsed on xi3 = 0. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Collapsed nodes on xi2 = 0 have 2 additional + scale factors cos and sin coefficients of the inclination angle. + Need to create a new template for each sector around axis giving common + nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and + add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 :param nodeScaleFactorOffset1: offset of node scale factors at axis on xi1=1 :param nodeScaleFactorOffsetUp: offset of first scale factor for inclination at pole, @@ -252,82 +252,74 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs ''' # start with full tricubic eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) - for n in [ 4, 5, 6, 7 ]: - eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 - if not self._useCrossDerivatives: - eft.setFunctionNumberOfTerms(n*8 + 6, 0) # D2_DS1DS3 - eft.setFunctionNumberOfTerms(n*8 + 7, 0) # D2_DS2DS3 - eft.setFunctionNumberOfTerms(n*8 + 8, 0) # D3_DS1DS2DS3 + for n in [ 2, 3, 6, 7 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) + if n > 3: + eft.setFunctionNumberOfTerms(n*8 + 6, 0) + eft.setFunctionNumberOfTerms(n*8 + 7, 0) + eft.setFunctionNumberOfTerms(n*8 + 8, 0) # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) - # remap parameters on xi3 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) - # switch D_DS2 and D_DS and negate at bottom - remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, [1]) ]) - + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) + for layer in range(2): so = layer*10 + 1 - ln = layer*2 + 1 - if layer == 0: - # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) - - else: - # 2 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) - + ln = layer*4 + 1 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) # zero other cross derivative parameters - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - ln = layer*2 + 2 - if layer == 0: - # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) - - else: - # 2 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) - + ln = layer*4 + 2 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) # zero other cross derivative parameters - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [ 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, []) ]) + remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + + # remap parameters on xi3 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 3, 4 ], Node.VALUE_LABEL_D_DS1, []) + + so = 20 + 1 + remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 4]), (Node.VALUE_LABEL_D_DS3, [so + 5]) ]) + remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) - ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] + ln_map = [ 1, 1, 2, 2, 1, 1, 3, 4 ] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronTop: Failed to validate eft' return eft def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a - solid sphere pole, with xi1 collapsed on xi3 = 0, xi2 collapsed on xi3 = 1 and - xi3 radial. Each collapsed node has 3 scale factors giving the cos, sin coefficients - of the radial line from global derivatives, plus the arc subtended by the element in - radians, so the circumferential direction is rounded. Collapsed node on xi2 = 1 has - additional 2 scale factors giving the cos, sin coefficients of the angle of inclination - from global z-axis. Need to create a new template for each sector around axis giving common + solid sphere pole, with xi1 and xi3 collapsed on xi2 = 1, and xi1 collapsed on xi3 = 0. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Collapsed nodes on xi2 = 1 have 2 additional + scale factors cos and sin coefficients of the inclination angle. + Need to create a new template for each sector around axis giving common nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and add 100 for each radial line around axis. :param nodeScaleFactorOffset0: offset of node scale factors at axis on xi1=0 @@ -338,69 +330,66 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 ''' # start with full tricubic eft = self._mesh.createElementfieldtemplate(self._tricubicHermiteBasis) - for n in [ 4, 5, 6, 7 ]: - eft.setFunctionNumberOfTerms(n*8 + 4, 0) # D2_DS1DS2 - if not self._useCrossDerivatives: - eft.setFunctionNumberOfTerms(n*8 + 6, 0) # D2_DS1DS3 - eft.setFunctionNumberOfTerms(n*8 + 7, 0) # D2_DS2DS3 - eft.setFunctionNumberOfTerms(n*8 + 8, 0) # D3_DS1DS2DS3 + for n in [ 0, 1, 4, 5 ]: + eft.setFunctionNumberOfTerms(n*8 + 4, 0) + if n > 1: + eft.setFunctionNumberOfTerms(n*8 + 6, 0) + eft.setFunctionNumberOfTerms(n*8 + 7, 0) + eft.setFunctionNumberOfTerms(n*8 + 8, 0) # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2 ]) + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) - # remap parameters on xi3 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, []) + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS3, []) for layer in range(2): - so = layer*6 + 1 - ln = layer*2 + 1 - if layer == 0: - # 2 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) - else: - # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [ so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [ 1, so + 4]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 2, so + 3, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, []) ]) # switch D_DS2 and D_DS3 at top + so = layer*10 + 1 + ln = layer*4 + 3 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 4]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 1, so + 3, so + 5]) ]) # zero other cross derivative parameters - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - ln = layer*2 + 2 - if layer == 0: - # 2 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, so + 4), (Node.VALUE_LABEL_D_DS3, so + 5) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) - else: - # 3 terms for d/dxi3 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [ so + 6, so + 10]), (Node.VALUE_LABEL_D_DS2, [ so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [ 1, so + 9]) ]) - # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [ 1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS3, []) ]) # switch D_DS2 and D_DS3 at top + ln = layer*4 + 4 + # 3 terms for d/dxi2 via general linear map: + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 9]) ]) + # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 6, so + 8, so + 10]) ]) # zero other cross derivative parameters - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [ 5, 6, 7, 8 ], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [ 7, 8 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) + remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS2, [1]) ]) - ln_map = [ 1, 1, 2, 2, 3, 4, 3, 4 ] + # remap parameters on xi3 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS1, []) + + so = 20 + 1 + remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) + remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) + remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 4]), (Node.VALUE_LABEL_D_DS3, [so + 5]) ]) + remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + + ln_map = [ 1, 1, 2, 2, 3, 4, 2, 2 ] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronBottom: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronTop: Failed to validate eft' return eft + def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, nodeScaleFactorOffsetUp): ''' Create a tricubic hermite element field for a solid pyramid for the axis of a From 22709683e1e8403333e67bbd8df2810baad40822 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Aug 2018 14:31:35 +1200 Subject: [PATCH 8/9] Modify radius calculation and set diameter to 0 for negative input --- scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index 110c1e57..ffe33eaf 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -63,7 +63,7 @@ def checkOptions(options): if options['Number of elements around'] < 4: options['Number of elements around'] = 4 if options['Diameter'] < 0.0: - options['Diameter'] = 1.0 + options['Diameter'] = 0.0 @classmethod def generateBaseMesh(cls, region, options): @@ -77,8 +77,7 @@ def generateBaseMesh(cls, region, options): elementsCountUp = options['Number of elements up'] elementsCountRadial = options['Number of elements radial'] useCrossDerivatives = options['Use cross derivatives'] - diameter = options['Diameter'] - radius = diameter/2.0 + radius = 0.5*options['Diameter'] fm = region.getFieldmodule() fm.beginChange() From 5f52b9e0f103989b8452eb96dc9ba9677294bc7e Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Mon, 13 Aug 2018 17:14:13 +1200 Subject: [PATCH 9/9] Amend to avoid scale factor repetition for tetrahedrons and pyramids --- .../meshtypes/meshtype_3d_solidsphere1.py | 42 +++--- .../utils/eftfactory_tricubichermite.py | 120 +++++++++--------- 2 files changed, 84 insertions(+), 78 deletions(-) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py index ffe33eaf..7b4c6215 100644 --- a/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py +++ b/scaffoldmaker/meshtypes/meshtype_3d_solidsphere1.py @@ -352,12 +352,12 @@ def generateBaseMesh(cls, region, options): radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAround), math.sin(radiansAround),radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansInclineNext), math.sin(radiansInclineNext) ] result2 = element.setScaleFactors(eft1, scalefactors) # print('Tetrahedron Bottom element', elementIdentifier, result1, result2, nodeIdentifiers) @@ -379,11 +379,11 @@ def generateBaseMesh(cls, region, options): radiansAround = va*radiansPerElementAround radiansAroundNext = vb*radiansPerElementAround scalefactors = [ - -1.0, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) + -1.0, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansInclineNext), math.sin(radiansInclineNext) ] result2 = element.setScaleFactors(eft4, scalefactors) # print('pyramid bottom element', elementIdentifier, result1, result2, nodeIdentifiers) @@ -452,12 +452,12 @@ def generateBaseMesh(cls, region, options): radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAround), math.sin(radiansAround),radiansPerElementAround, - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansInclineNext), math.sin(radiansInclineNext) ] result2 = element.setScaleFactors(eft3, scalefactors) # print('Tetrahedron top element', elementIdentifier, result1, result2, nodeIdentifiers) @@ -481,10 +481,10 @@ def generateBaseMesh(cls, region, options): radiansAroundNext = vb*radiansPerElementAround scalefactors = [ -1.0, - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansIncline), math.sin(radiansIncline), - math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext), - math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, math.cos(radiansInclineNext), math.sin(radiansInclineNext) + math.cos(radiansAround), math.sin(radiansAround), radiansPerElementAround, + math.cos(radiansAroundNext), math.sin(radiansAroundNext), radiansPerElementAround, + math.cos(radiansIncline), math.sin(radiansIncline), + math.cos(radiansInclineNext), math.sin(radiansInclineNext) ] result2 = element.setScaleFactors(eft5, scalefactors) # print('pyramid top element', elementIdentifier, result1, result2, nodeIdentifiers) diff --git a/scaffoldmaker/utils/eftfactory_tricubichermite.py b/scaffoldmaker/utils/eftfactory_tricubichermite.py index abc68c11..752b1058 100644 --- a/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -237,10 +237,10 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a solid sphere pole, with xi1 and xi3 collapsed on xi2 = 0, and xi1 collapsed on xi3 = 0. - Each collapsed node has 3 scale factors giving the cos, sin coefficients of - the radial line from global derivatives, plus the arc subtended by the element in radians, - so the circumferential direction is rounded. Collapsed nodes on xi2 = 0 have 2 additional - scale factors cos and sin coefficients of the inclination angle. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Collapsed nodes on xi2 = 0 have 2 additional + scale factors cos and sin coefficients of the inclination angle. Need to create a new template for each sector around axis giving common nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and add 100 for each radial line around axis. @@ -261,24 +261,25 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4]) # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) for layer in range(2): - so = layer*10 + 1 + soAround = 1 + soUp = layer*2 + 1 + 12 ln = layer*4 + 1 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 1, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 2, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [1, soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 2, soAround + 3 , soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 1, soAround + 3, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -286,9 +287,9 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs ln = layer*4 + 2 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 4, soUp + 2] ), (Node.VALUE_LABEL_D_DS2, [soAround + 5, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [1, soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 5, soAround + 6, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 4, soAround + 6, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -299,11 +300,11 @@ def createEftTetrahedronBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffs # remap parameters on xi3 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 3, 4 ], Node.VALUE_LABEL_D_DS1, []) - so = 20 + 1 - remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) - remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 4]), (Node.VALUE_LABEL_D_DS3, [so + 5]) ]) - remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + soAround = 1 + 6 + remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [soAround + 1]), (Node.VALUE_LABEL_D_DS3, [soAround + 2]) ]) + remapEftNodeValueLabel(eft, [ 3 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 2, soAround + 3]), (Node.VALUE_LABEL_D_DS3, [soAround + 1, soAround + 3]) ]) + remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [soAround + 4]), (Node.VALUE_LABEL_D_DS3, [soAround + 5]) ]) + remapEftNodeValueLabel(eft, [ 4 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 5, soAround + 6]), (Node.VALUE_LABEL_D_DS3, [soAround + 4, soAround + 6]) ]) ln_map = [ 1, 1, 2, 2, 1, 1, 3, 4 ] remapEftLocalNodes(eft, 4, ln_map) @@ -315,10 +316,10 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 ''' Create a tricubic hermite element field for a solid tetrahedron for the axis of a solid sphere pole, with xi1 and xi3 collapsed on xi2 = 1, and xi1 collapsed on xi3 = 0. - Each collapsed node has 3 scale factors giving the cos, sin coefficients of - the radial line from global derivatives, plus the arc subtended by the element in radians, - so the circumferential direction is rounded. Collapsed nodes on xi2 = 1 have 2 additional - scale factors cos and sin coefficients of the inclination angle. + Each collapsed node has 3 scale factors giving the cos, sin coefficients of + the radial line from global derivatives, plus the arc subtended by the element in radians, + so the circumferential direction is rounded. Collapsed nodes on xi2 = 1 have 2 additional + scale factors cos and sin coefficients of the inclination angle. Need to create a new template for each sector around axis giving common nodeScaleFactorOffset values on common faces. Suggestion is to start at 0 and add 100 for each radial line around axis. @@ -339,24 +340,25 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3 ]) + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4]) # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS3, []) for layer in range(2): - so = layer*10 + 1 + soAround = 1 + 6 + soUp = layer*2 + 1 + 12 ln = layer*4 + 3 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 1, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 2, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3 , soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -364,9 +366,9 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 ln = layer*4 + 4 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 4, soUp + 2] ), (Node.VALUE_LABEL_D_DS2, [1, soAround + 5, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -376,12 +378,12 @@ def createEftTetrahedronTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1 # remap parameters on xi3 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2 ], Node.VALUE_LABEL_D_DS1, []) - - so = 20 + 1 - remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 1]), (Node.VALUE_LABEL_D_DS3, [so + 2]) ]) - remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3]), (Node.VALUE_LABEL_D_DS3, [so + 1, so + 3]) ]) - remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [so + 4]), (Node.VALUE_LABEL_D_DS3, [so + 5]) ]) - remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, so + 5, so + 6]), (Node.VALUE_LABEL_D_DS3, [so + 4, so + 6]) ]) + + soAround = 1 + remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [soAround + 1]), (Node.VALUE_LABEL_D_DS3, [soAround + 2]) ]) + remapEftNodeValueLabel(eft, [ 1 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 2, soAround + 3]), (Node.VALUE_LABEL_D_DS3, [soAround + 1, soAround + 3]) ]) + remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D_DS3, [ (Node.VALUE_LABEL_D_DS1, [soAround + 4]), (Node.VALUE_LABEL_D_DS3, [soAround + 5]) ]) + remapEftNodeValueLabel(eft, [ 2 ], Node.VALUE_LABEL_D2_DS1DS3, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 5, soAround + 6]), (Node.VALUE_LABEL_D_DS3, [soAround + 4, soAround + 6]) ]) ln_map = [ 1, 1, 2, 2, 3, 4, 2, 2 ] remapEftLocalNodes(eft, 4, ln_map) @@ -416,22 +418,23 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4 ]) + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4]) # remap parameters on xi2 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS3, []) for layer in range(2): - so = layer*10 + 1 + soAround = 1 + soUp = layer*2 + 6 + 1 ln = layer*4 + 1 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [1, so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 1, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 2, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [1, soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 2, soAround + 3 , soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 1, soAround + 3, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -439,9 +442,9 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, ln = layer*4 + 2 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [1, so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 4, soUp + 2] ), (Node.VALUE_LABEL_D_DS2, [soAround + 5, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [1, soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 5, soAround + 6, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [soAround + 4, soAround + 6, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -449,6 +452,7 @@ def createEftPyramidBottom(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, ln_map = [ 1, 1, 2, 3, 1, 1, 4, 5 ] remapEftLocalNodes(eft, 5, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' return eft @@ -478,22 +482,23 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no # GRC: allow scale factor identifier for global -1.0 to be prescribed setEftScaleFactorIds(eft, [1], [ - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4 ]) + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 2, nodeScaleFactorOffset0 + 3, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 2, nodeScaleFactorOffset1 + 3, + nodeScaleFactorOffsetUp + 1, nodeScaleFactorOffsetUp + 2, + nodeScaleFactorOffsetUp + 3, nodeScaleFactorOffsetUp + 4]) # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [ 3, 4, 7, 8 ], Node.VALUE_LABEL_D_DS3, []) for layer in range(2): - so = layer*10 + 1 + soAround = 1 + soUp = layer*2 + 1 + 6 ln = layer*4 + 3 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 1, so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 2, so + 5]), (Node.VALUE_LABEL_D_DS3, [so + 4]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 1, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 2, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 2, so + 3 , so + 5]), (Node.VALUE_LABEL_D_DS2, [1, so + 1, so + 3, so + 5]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 2, soAround + 3 , soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 1, soAround + 3, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -501,9 +506,9 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no ln = layer*4 + 4 # 3 terms for d/dxi2 via general linear map: - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, so + 6, so + 10] ), (Node.VALUE_LABEL_D_DS2, [1, so + 7, so + 10]), (Node.VALUE_LABEL_D_DS3, [so + 9]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D_DS2, [ (Node.VALUE_LABEL_D_DS1, [1, soAround + 4, soUp + 2] ), (Node.VALUE_LABEL_D_DS2, [1, soAround + 5, soUp + 2]), (Node.VALUE_LABEL_D_DS3, [soUp + 1]) ]) # 2 terms for cross derivative 1 2 to correct circular apex: -sin(theta).phi, cos(theta).phi - remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [so + 7, so + 8, so + 10]), (Node.VALUE_LABEL_D_DS2, [1, so + 6, so + 8, so + 10]) ]) + remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS2, [ (Node.VALUE_LABEL_D_DS1, [soAround + 5, soAround + 6, soUp + 2]), (Node.VALUE_LABEL_D_DS2, [1, soAround + 4, soAround + 6, soUp + 2]) ]) # zero other cross derivative parameters remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [ ln ], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -511,6 +516,7 @@ def createEftPyramidTop(self, nodeScaleFactorOffset0, nodeScaleFactorOffset1, no ln_map = [ 1, 2, 3, 3, 4, 5, 3, 3 ] remapEftLocalNodes(eft, 5, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftPyramidTop: Failed to validate eft' return eft