diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py index eb16cc04..318113ca 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_heartventricles3.py @@ -23,7 +23,7 @@ from scaffoldmaker.utils.interpolation import computeCubicHermiteDerivativeScaling, getCubicHermiteArcLength, interpolateSampleCubicHermite, \ sampleCubicHermiteCurves, smoothCubicHermiteDerivativesLine from scaffoldmaker.utils.meshrefinement import MeshRefinement -from scaffoldmaker.utils.shieldmesh import ShieldMesh +from scaffoldmaker.utils.shieldmesh import ShieldMesh2D from scaffoldmaker.utils.tracksurface import TrackSurface, calculate_surface_axes @@ -492,7 +492,7 @@ def generateBaseMesh(cls, region, options): lad1 = [ [-s for s in d ] for d in lad1 ] lad3 = [ vector.setMagnitude(d, lvFreeWallThickness) for d in lad3 ] - rvShield = ShieldMesh(elementsCountAroundRVFreeWall, elementsCountUpRVFreeWall, 0) + rvShield = ShieldMesh2D(elementsCountAroundRVFreeWall, elementsCountUpRVFreeWall, 0) rvx = rvShield.px rvd1 = rvShield.pd1 rvd2 = rvShield.pd2 @@ -603,7 +603,7 @@ def generateBaseMesh(cls, region, options): # LV free wall elementsCountUpLV = elementsCountUpLVFreeWall + elementsCountUpLVApex - lvShield = ShieldMesh(elementsCountAroundLVFreeWall, elementsCountUpLV, elementsCountUpLVApex, lvTrackSurface) + lvShield = ShieldMesh2D(elementsCountAroundLVFreeWall, elementsCountUpLV, elementsCountUpLVApex, lvTrackSurface) lvx = lvShield.px lvd1 = lvShield.pd1 lvd2 = lvShield.pd2 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py new file mode 100644 index 00000000..5d59cf1c --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py @@ -0,0 +1,203 @@ +""" +Generates a solid sphere (spheroid/ellipsoid in general) using a ShieldMesh of all cube elements, + with variable numbers of elements across axes and shell directions. +""" + +from __future__ import division + +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from scaffoldmaker.utils.spheremesh import SphereMesh, SphereShape +from scaffoldmaker.utils import vector + + +class MeshType_3d_solidsphere2(Scaffold_base): + """ +Generates a solid sphere using a ShieldMesh of all cube elements, +with variable numbers of elements across axes and shell directions. + """ + + @staticmethod + def getName(): + return '3D Solid Sphere 2' + + @classmethod + def getDefaultOptions(cls, parameterSetName='Default'): + options = { + 'Number of elements across axis 1': 4, + 'Number of elements across axis 2': 4, + 'Number of elements across axis 3': 4, + 'Number of elements across shell': 0, + 'Number of elements across transition': 1, + 'Radius1': 1.0, + 'Radius2': 1.0, + 'Radius3': 1.0, + 'Shell element thickness proportion': 1.0, + 'Range of elements required in direction 1': [0, 4], + 'Range of elements required in direction 2': [0, 4], + 'Range of elements required in direction 3': [0, 4], + 'Box derivatives': [1, 2, 3], + 'Use cross derivatives': False, + 'Refine': False, + 'Refine number of elements': 1, + } + return options + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Number of elements across axis 1', + 'Number of elements across axis 2', + 'Number of elements across axis 3', + # 'Number of elements across shell', + # 'Number of elements across transition', + 'Radius1', + 'Radius2', + 'Radius3', + # 'Shell element thickness proportion', + 'Range of elements required in direction 1', + 'Range of elements required in direction 2', + 'Range of elements required in direction 3', + 'Box derivatives', + 'Refine', + 'Refine number of elements' + ] + + @classmethod + def checkOptions(cls, options): + dependentChanges = False + + min1, min2, min3 = 4, 4, 4 + co1, co2, co3 = 1, 1, 1 + + if options['Number of elements across axis 1'] < min1: + options['Number of elements across axis 1'] = min1 + if options['Number of elements across axis 2'] < min2: + options['Number of elements across axis 2'] = min2 + if options['Number of elements across axis 3'] < min3: + options['Number of elements across axis 3'] = min3 + + if options['Number of elements across axis 1'] % 2: + options['Number of elements across axis 1'] += co1 + if options['Number of elements across axis 2'] % 2: + options['Number of elements across axis 2'] += co2 + if options['Number of elements across axis 3'] % 2: + options['Number of elements across axis 3'] += co3 + + for radius in ['Radius1', 'Radius2', 'Radius3', ]: + if options[radius] <= 0: + options[radius] = 1.0 + + if not all(abs(d) in [1, 2, 3] for d in options['Box derivatives']): + options['Box derivatives'] = [1, 2, 3] + if len(options['Box derivatives']) > len(set(options['Box derivatives'])): + options['Box derivatives'] = [1, 2, 3] + + for i in range(1, 4): + if options['Range of elements required in direction {}'.format(i)][1] == \ + options['Number of elements across axis {}'.format(i)] - 1: + options['Range of elements required in direction {}'.format(i)][1] = \ + options['Number of elements across axis {}'.format(i)] + + nm = 4 + for i in range(1, nm): + if options['Range of elements required in direction {}'.format(i)][0] == 1: + options['Range of elements required in direction {}'.format(i)][0] = 0 + + maxelems = [options['Number of elements across axis 1'], + options['Number of elements across axis 2'], + options['Number of elements across axis 3']] + ranges = [options['Range of elements required in direction 1'], + options['Range of elements required in direction 2'], + options['Range of elements required in direction 3']] + for i in range(3): + if ranges[i][1] > maxelems[i] or ranges[i][1] <= max(1, ranges[i][0]): + ranges[i][1] = maxelems[i] + for i in range(3): + if ranges[i][0] >= min(maxelems[i] - 1, ranges[i][1]) or ranges[i][0] < 1: + ranges[i][0] = 0 + + options['Range of elements required in direction 1'] = ranges[0] + options['Range of elements required in direction 2'] = ranges[1] + options['Range of elements required in direction 3'] = ranges[2] + + # if options['Number of elements across transition'] < 1: + # options['Number of elements across transition'] = 1 + # Rcrit = min(options['Number of elements across major']-4, options['Number of elements across minor']-4)//2 + # if options['Number of elements across shell'] + options['Number of elements across transition'] - 1 > Rcrit: + # dependentChanges = True + # options['Number of elements across shell'] = Rcrit + # options['Number of elements across transition'] = 1 + # + # if options['Shell element thickness proportion'] < 0.15: + # options['Shell element thickness proportion'] = 1.0 + + return dependentChanges + + @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 + """ + + elementsCountAcrossAxis1 = options['Number of elements across axis 1'] + elementsCountAcrossAxis2 = options['Number of elements across axis 2'] + elementsCountAcrossAxis3 = options['Number of elements across axis 3'] + + elementsCountAcrossShell = options['Number of elements across shell'] + elementsCountAcrossTransition = options['Number of elements across transition'] + shellProportion = options['Shell element thickness proportion'] + radius = [options['Radius1'], options['Radius2'], options['Radius3']] + useCrossDerivatives = options['Use cross derivatives'] + rangeOfRequiredElements = [options['Range of elements required in direction 1'], + options['Range of elements required in direction 2'], + options['Range of elements required in direction 3']] + sphereBoxDerivatives = [-options['Box derivatives'][0], options['Box derivatives'][1], + options['Box derivatives'][2]] # To make the values more intuitive for the user but + # consistent with [back, right, up] + # sphereBoxDerivatives = [1, 3, 2] # consistent with default derivatives of cylinder mesh. + # This is the default value that is used for base sphere. + sphere_shape = SphereShape.SPHERE_SHAPE_FULL + + fm = region.getFieldmodule() + coordinates = findOrCreateFieldCoordinates(fm) + + mesh = fm.findMeshByDimension(3) + boxGroup = AnnotationGroup(region, ("box group", "")) + boxMeshGroup = boxGroup.getMeshGroup(mesh) + transitionGroup = AnnotationGroup(region, ("transition group", "")) + transitionMeshGroup = transitionGroup.getMeshGroup(mesh) + meshGroups = [boxMeshGroup, transitionMeshGroup] + annotationGroups = [boxGroup, transitionGroup] + + centre = [0.0, 0.0, 0.0] + axis1 = [1.0, 0.0, 0.0] + axis2 = [0.0, 1.0, 0.0] + axis3 = [0.0, 0.0, 1.0] + axes = [vector.scaleVector(axis1, radius[0]), + vector.scaleVector(axis2, radius[1]), + vector.scaleVector(axis3, radius[2])] + elementsCountAcross = [elementsCountAcrossAxis1, elementsCountAcrossAxis2, elementsCountAcrossAxis3] + + sphere1 = SphereMesh(fm, coordinates, centre, axes, elementsCountAcross, + elementsCountAcrossShell, elementsCountAcrossTransition, shellProportion, + sphereShape=sphere_shape, rangeOfRequiredElements=rangeOfRequiredElements, + boxDerivatives=sphereBoxDerivatives, useCrossDerivatives=False, meshGroups=meshGroups) + + return annotationGroups + + @classmethod + def refineMesh(cls, meshRefinement, options): + """ + Refine source mesh into separate region, with change of basis. + :param meshRefinement: MeshRefinement, which knows source and target region. + :param options: Dict containing options. See getDefaultOptions(). + """ + assert isinstance(meshRefinement, MeshRefinement) + refineElementsCount = options['Refine number of elements'] + meshRefinement.refineAllElementsCubeStandard3d(refineElementsCount, refineElementsCount, refineElementsCount) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index b53ef1c3..350c312d 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -36,6 +36,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_smallintestine1 import MeshType_3d_smallintestine1 from scaffoldmaker.meshtypes.meshtype_3d_solidcylinder1 import MeshType_3d_solidcylinder1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 +from scaffoldmaker.meshtypes.meshtype_3d_solidsphere2 import MeshType_3d_solidsphere2 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_stellate1 import MeshType_3d_stellate1 @@ -83,6 +84,7 @@ def __init__(self): MeshType_3d_ostium1, MeshType_3d_smallintestine1, MeshType_3d_solidsphere1, + MeshType_3d_solidsphere2, MeshType_3d_solidcylinder1, MeshType_3d_sphereshell1, MeshType_3d_sphereshellseptum1, diff --git a/src/scaffoldmaker/utils/cylindermesh.py b/src/scaffoldmaker/utils/cylindermesh.py index c93d2147..14e0cf1e 100644 --- a/src/scaffoldmaker/utils/cylindermesh.py +++ b/src/scaffoldmaker/utils/cylindermesh.py @@ -15,7 +15,7 @@ from scaffoldmaker.utils.interpolation import sampleCubicHermiteCurves, interpolateSampleCubicHermite, \ smoothCubicHermiteDerivativesLine from scaffoldmaker.utils.mirror import Mirror -from scaffoldmaker.utils.shieldmesh import ShieldMesh, ShieldShape, ShieldRimDerivativeMode +from scaffoldmaker.utils.shieldmesh import ShieldMesh2D, ShieldShape2D, ShieldRimDerivativeMode class CylinderShape(Enum): @@ -214,11 +214,11 @@ def createCylinderMesh3d(self, fieldModule, coordinates): elementsCountRim = self._elementsCountAcrossRim - shieldMode = ShieldShape.SHIELD_SHAPE_FULL if self._cylinderShape is self._cylinderShape.CYLINDER_SHAPE_FULL \ - else ShieldShape.SHIELD_SHAPE_LOWER_HALF + shieldMode = ShieldShape2D.SHIELD_SHAPE_FULL if self._cylinderShape is self._cylinderShape.CYLINDER_SHAPE_FULL \ + else ShieldShape2D.SHIELD_SHAPE_LOWER_HALF ellipseShape = EllipseShape.Ellipse_SHAPE_FULL \ if self._cylinderShape is self._cylinderShape.CYLINDER_SHAPE_FULL else EllipseShape.Ellipse_SHAPE_LOWER_HALF - self._shield = ShieldMesh(self._elementsCountAcrossMinor, self._elementsCountAcrossMajor, elementsCountRim, + self._shield = ShieldMesh2D(self._elementsCountAcrossMinor, self._elementsCountAcrossMajor, elementsCountRim, None, self._elementsCountAlong, shieldMode, shieldType=ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND) @@ -425,9 +425,9 @@ def __init__(self, centre, majorAxis, minorAxis, self.coreMajorRadius = coreMajorRadius self.coreMinorRadius = coreMinorRadius elementsCountRim = elementsCountAcrossShell + elementsCountAcrossTransition - 1 - shieldMode = ShieldShape.SHIELD_SHAPE_FULL if ellipseShape is EllipseShape.Ellipse_SHAPE_FULL\ - else ShieldShape.SHIELD_SHAPE_LOWER_HALF - shield = ShieldMesh(elementsCountAcrossMinor, elementsCountAcrossMajor, elementsCountRim, + shieldMode = ShieldShape2D.SHIELD_SHAPE_FULL if ellipseShape is EllipseShape.Ellipse_SHAPE_FULL\ + else ShieldShape2D.SHIELD_SHAPE_LOWER_HALF + shield = ShieldMesh2D(elementsCountAcrossMinor, elementsCountAcrossMajor, elementsCountRim, None, 1, shieldMode, shieldType=ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND) self.elementsCountAcrossMajor = elementsCountAcrossMajor @@ -464,6 +464,7 @@ def generate2DEllipseMesh(self): self.smoothTransitionRims() if self.ellipseShape == EllipseShape.Ellipse_SHAPE_FULL: self.generateNodesForUpperHalf() + self.calculateD2() def generateBase1DMesh(self, rx): """ @@ -750,6 +751,22 @@ def generateNodesForUpperHalf(self): self.pd3[2 * self.elementsCountUp - n2][n1] = mirror.mirrorVector( self.pd3[n2][n1]) + def calculateD2(self): + """ + :return: + """ + btx = self.px + btd2 = self.pd2 + + nte = normalToEllipse(self.majorAxis, self.minorAxis) + for n2 in range(self.elementsCountAcrossMajor + 1): + for n1 in range(self.elementsCountAcrossMinor + 1): + if btx[n2][n1]: + btd2[n2][n1] = nte + + def getShield(self): + return self.__shield + def createEllipsePerimeter(centre, majorAxis, minorAxis, elementsCountAround, height): """ diff --git a/src/scaffoldmaker/utils/shieldmesh.py b/src/scaffoldmaker/utils/shieldmesh.py index 058c0400..b202f833 100644 --- a/src/scaffoldmaker/utils/shieldmesh.py +++ b/src/scaffoldmaker/utils/shieldmesh.py @@ -1,10 +1,11 @@ -''' +""" Class for generating a shield-shaped mesh, with a regular flat top but a rounded bottom formed by having two points where 3 square elements merge to form a triangle. -''' -from __future__ import division +""" +from __future__ import division +import copy from enum import Enum from opencmiss.zinc.element import Element @@ -19,22 +20,38 @@ from scaffoldmaker.utils.tracksurface import TrackSurface, calculate_surface_axes -class ShieldShape(Enum): +class ShieldShape2D(Enum): SHIELD_SHAPE_FULL = 1 SHIELD_SHAPE_LOWER_HALF = 2 + +class ShieldShape3D(Enum): + SHIELD_SHAPE_FULL = 1 + SHIELD_SHAPE_HALF_AAP = 2 # AAP is a hemisphere where x_a3>=0 + SHIELD_SHAPE_OCTANT_PPP = 3 # PPP means positive axis1, positive axis2, positive axis3. + SHIELD_SHAPE_OCTANT_PNP = 4 + SHIELD_SHAPE_OCTANT_NNP = 5 + SHIELD_SHAPE_OCTANT_NPP = 6 + SHIELD_SHAPE_OCTANT_PPN = 7 + SHIELD_SHAPE_OCTANT_PNN = 8 + SHIELD_SHAPE_OCTANT_NNN = 9 + SHIELD_SHAPE_OCTANT_NPN = 10 + + class ShieldRimDerivativeMode(Enum): SHIELD_RIM_DERIVATIVE_MODE_AROUND = 1 # rim derivatives are d1 anticlockwise around shield, d3 outward SHIELD_RIM_DERIVATIVE_MODE_REGULAR = 2 # rim derivatives d1, d2 match interior nodes for regular elements -class ShieldMesh: - ''' + +class ShieldMesh2D: + """ Shield mesh generator. - ''' + """ - def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, trackSurface : TrackSurface=None, - elementsCountAlong=1, shieldMode=ShieldShape.SHIELD_SHAPE_LOWER_HALF, shieldType=ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR): - ''' + def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, trackSurface: TrackSurface=None, + elementsCountAlong=1, shieldMode=ShieldShape2D.SHIELD_SHAPE_LOWER_HALF, + shieldType=ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR): + """ Data structure for defining a shield-shaped mesh which is flat on the top and rounded around the bottom and/or the same mirror mirrored on top. It is represented as a regular box of elementsCountAcross x elementsCountUp @@ -56,20 +73,24 @@ def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, t Extra elements up add regular rows of nodes/elements on top, and extra non-rim elements across add regular columns of nodes/elements up the centre. :param elementsCountAcross: Number of elements across top of shield. Must be at least 4 + 2*elementsCountRim. - :param elementsCountUpFull: Number of elements up central axis of shield. Must be at least 2 + elementsCountRim if half and 4 + 2*elementsCountRim if full. - :param elementsCountAlong: Number of elements through wall for ventricle case (only 1 element is supported) and along cylinder axis in cylinder case. + :param elementsCountUpFull: Number of elements up central axis of shield. Must be at least 2 + elementsCountRim + if half and 4 + 2*elementsCountRim if full. + :param elementsCountAlong: Number of elements through wall for ventricle case (only 1 element is supported) and + along cylinder axis in cylinder case. :param elementsCountRim: Number of elements around bottom rim (not top) outside of 'triple points'. :param trackSurface: Optional trackSurface to store or restrict points to. :param shieldMode: It determines if the shield is full or just part of it. - :param shieldType: To distinguish between cylinder and ventricle type. Derivatives and directions are chosen differently for two cases. - ''' + :param shieldType: To distinguish between cylinder and ventricle type. Derivatives and directions are chosen + differently for two cases. + """ assert elementsCountRim >= 0 assert elementsCountAlong >= 1 assert elementsCountAcross >= (elementsCountRim + 4) assert elementsCountUpFull >= (elementsCountRim + 2) self.elementsCountAcross = elementsCountAcross self.elementsCountUpFull = elementsCountUpFull - elementsCountUp = elementsCountUpFull//2 if shieldMode == ShieldShape.SHIELD_SHAPE_FULL else elementsCountUpFull + elementsCountUp = elementsCountUpFull//2 if shieldMode == ShieldShape2D.SHIELD_SHAPE_FULL\ + else elementsCountUpFull self.elementsCountUp = elementsCountUp self.elementsCountRim = elementsCountRim self.elementsCountAlong = elementsCountAlong @@ -79,29 +100,30 @@ def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, t self.trackSurface = trackSurface self._mode = shieldMode self._type = shieldType - self.px = [ [] for _ in range(elementsCountAlong+1) ] - self.pd1 = [ [] for _ in range(elementsCountAlong+1) ] - self.pd2 = [ [] for _ in range(elementsCountAlong+1) ] - self.pd3 = [ [] for _ in range(elementsCountAlong+1) ] - self.nodeId = [ [] for _ in range(elementsCountAlong+1) ] + self.px = [[] for _ in range(elementsCountAlong+1)] + self.pd1 = [[] for _ in range(elementsCountAlong+1)] + self.pd2 = [[] for _ in range(elementsCountAlong+1)] + self.pd3 = [[] for _ in range(elementsCountAlong+1)] + self.nodeId = [[] for _ in range(elementsCountAlong+1)] for n3 in range(elementsCountAlong+1): for n2 in range(elementsCountUpFull + 1): - for p in [ self.px[n3], self.pd1[n3], self.pd2[n3], self.pd3[n3], self.nodeId[n3] ]: - p.append([ None ]*(elementsCountAcross + 1)) + for p in [self.px[n3], self.pd1[n3], self.pd2[n3], self.pd3[n3], self.nodeId[n3]]: + p.append([None]*(elementsCountAcross + 1)) if trackSurface: - self.pProportions = [ [ None ]*(elementsCountAcross + 1) for n2 in range(elementsCountUp + 1) ] + self.pProportions = [[None]*(elementsCountAcross + 1) for n2 in range(elementsCountUp + 1)] if shieldType == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.elementId = [ [[ None ]*elementsCountAcross for n2 in range(elementsCountUpFull)] for e3 in range(elementsCountAlong) ] + self.elementId = [[[None]*elementsCountAcross for n2 in range(elementsCountUpFull)] + for e3 in range(elementsCountAlong)] else: self.elementId = [[None] * elementsCountAcross for n2 in range(elementsCountUpFull)] def convertRimIndex(self, ix, rx=0): - ''' + """ Convert point index around the lower rim to n1, n2 across and up box. :param ix: index around from 0 to self.elementsCountAroundFull :param rx: rim index from 0 (around outside) to self.elementsCountRim :return: n1, n2 - ''' + """ assert 0 <= ix <= self.elementsCountAroundFull assert 0 <= rx <= self.elementsCountRim if ix <= self.elementsCountUpRegular: @@ -111,12 +133,11 @@ def convertRimIndex(self, ix, rx=0): return self.elementsCountAcross - rx, self.elementsCountUp - mx return self.elementsCountRim + ix - self.elementsCountUpRegular, rx - def getTriplePoints(self, n3): - ''' + """ Compute coordinates and derivatives of points where 3 square elements merge. :param n3: Index of through-wall coordinates to use. - ''' + """ n1a = self.elementsCountRim n1b = n1a + 1 n1c = n1a + 2 @@ -131,93 +152,122 @@ def getTriplePoints(self, n3): if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2a][n1c], self.px[n3][n2c][n1b]],[[(-self.pd1[n3][n2a][n1c][c] - self.pd3[n3][n2a][n1c][c]) for c in range(3)], self.pd1[n3][n2c][n1b]],2, arcLengthDerivatives=True)[0:2] + [self.px[n3][n2a][n1c], self.px[n3][n2c][n1b]], + [[(-self.pd1[n3][n2a][n1c][c] - self.pd3[n3][n2a][n1c][c]) for c in range(3)], self.pd1[n3][n2c][n1b]], + 2, arcLengthDerivatives=True)[0:2] ltx.append(tx[1]) # tx, td1 = sampleCubicHermiteCurves( - # [ self.px[n3][n2a][n1b], self.px[n3][n2c][n1c] ], [ [-self.pd3[n3][n2a][n1b][c] for c in range(3)], [ (self.pd1[n3][n2c][n1c][c] + self.pd3[n3][n2c][n1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] + # [ self.px[n3][n2a][n1b], self.px[n3][n2c][n1c] ], [ [-self.pd3[n3][n2a][n1b][c] for c in range(3)], + # [ (self.pd1[n3][n2c][n1c][c] + self.pd3[n3][n2c][n1c][c]) for c in range(3) ] ], + # 2, arcLengthDerivatives = True)[0:2] ltx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [ self.px[n3][n2c][n1a], self.px[n3][n2b][n1c] ], [ [ (self.pd1[n3][n2c][n1a][c] - self.pd3[n3][n2c][n1a][c]) for c in range(3) ], self.pd3[n3][n2b][n1c] ], 2, arcLengthDerivatives = True)[0:2] + [self.px[n3][n2c][n1a], self.px[n3][n2b][n1c]], + [[(self.pd1[n3][n2c][n1a][c] - self.pd3[n3][n2c][n1a][c]) for c in range(3)], + self.pd3[n3][n2b][n1c]], 2, arcLengthDerivatives=True)[0:2] ltx.append(tx[1]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2a][n1c], self.px[n3][n2c][n1b]], [[(-self.pd1[n3][n2a][n1c][c] + self.pd2[n3][n2a][n1c][c]) for c in range(3)],self.pd2[n3][n2c][n1b]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2a][n1c], self.px[n3][n2c][n1b]], + [[(-self.pd1[n3][n2a][n1c][c] + self.pd2[n3][n2a][n1c][c]) for c in range(3)], + self.pd2[n3][n2c][n1b]], 2, arcLengthDerivatives=True)[0: 2] ltx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2a][n1b], self.px[n3][n2c][n1c]], [self.pd2[n3][n2a][n1b], [(self.pd1[n3][n2c][n1c][c] + self.pd2[n3][n2c][n1c][c]) for c in range(3)]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2a][n1b], self.px[n3][n2c][n1c]], + [self.pd2[n3][n2a][n1b], [(self.pd1[n3][n2c][n1c][c] + self.pd2[n3][n2c][n1c][c]) for c in range(3)]], + 2, arcLengthDerivatives=True)[0: 2] ltx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2c][n1a], self.px[n3][n2b][n1c]], [[(self.pd1[n3][n2c][n1a][c] - self.pd2[n3][n2c][n1a][c]) for c in range(3)], self.pd1[n3][n2b][n1c]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2c][n1a], self.px[n3][n2b][n1c]], + [[(self.pd1[n3][n2c][n1a][c] - self.pd2[n3][n2c][n1a][c]) for c in range(3)], self.pd1[n3][n2b][n1c]], + 2, arcLengthDerivatives=True)[0: 2] ltx.append(tx[1]) - #x = [ (ltx[0][c] + ltx[1][c] + ltx[2][c])/3.0 for c in range(3) ] - x = [ (ltx[0][c] + ltx[2][c])/2.0 for c in range(3) ] + # x = [ (ltx[0][c] + ltx[1][c] + ltx[2][c])/3.0 for c in range(3) ] + x = [(ltx[0][c] + ltx[2][c])/2.0 for c in range(3)] if self.trackSurface: - p = self.trackSurface.findNearestPosition(x, startPosition=self.trackSurface.createPositionProportion(*(self.pProportions[n2b][n1c]))) + p = self.trackSurface.findNearestPosition(x, startPosition=self.trackSurface.createPositionProportion(*( + self.pProportions[n2b][n1c]))) self.pProportions[n2b][n1b] = self.trackSurface.getProportion(p) x, sd1, sd2 = self.trackSurface.evaluateCoordinates(p, derivatives=True) d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) self.pd3[n3][n2b][n1b] = d3 - self.px [n3][n2b][n1b] = x + self.px[n3][n2b][n1b] = x if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd3[n3][n2b][n1b] = [ (self.px[n3][n2b][n1c][c] - self.px[n3][n2b][n1b][c]) for c in range(3) ] - self.pd1[n3][n2b][n1b] = [ (self.px[n3][n2c][n1b][c] - self.px[n3][n2b][n1b][c]) for c in range(3) ] + self.pd3[n3][n2b][n1b] = [(self.px[n3][n2b][n1c][c] - self.px[n3][n2b][n1b][c]) for c in range(3)] + self.pd1[n3][n2b][n1b] = [(self.px[n3][n2c][n1b][c] - self.px[n3][n2b][n1b][c]) for c in range(3)] elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: self.pd1[n3][n2b][n1b] = [(self.px[n3][n2b][n1c][c] - self.px[n3][n2b][n1b][c]) for c in range(3)] self.pd2[n3][n2b][n1b] = [(self.px[n3][n2c][n1b][c] - self.px[n3][n2b][n1b][c]) for c in range(3)] if not self.trackSurface: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd2[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][n1b], self.pd1[n3][n2b][n1b])) + self.pd2[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][n1b], + self.pd1[n3][n2b][n1b])) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: - self.pd3[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][n1b], self.pd2[n3][n2b][n1b])) + self.pd3[n3][n2b][n1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][n1b], + self.pd2[n3][n2b][n1b])) # right rtx = [] if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: tx, td1 = sampleCubicHermiteCurves( - [ self.px[n3][n2a][m1c], self.px[n3][n2c][m1b] ], [ [ (self.pd1[n3][n2a][m1c][c] - self.pd3[n3][n2a][m1c][c]) for c in range(3) ], self.pd1[n3][n2c][m1b] ], 2, arcLengthDerivatives = True)[0:2] + [self.px[n3][n2a][m1c], self.px[n3][n2c][m1b]], + [[(self.pd1[n3][n2a][m1c][c] - self.pd3[n3][n2a][m1c][c]) for c in range(3)], self.pd1[n3][n2c][m1b]], + 2, arcLengthDerivatives=True)[0:2] rtx.append(tx[1]) # tx, td1 = sampleCubicHermiteCurves( - # [ self.px[n3][n2a][m1b], self.px[n3][n2c][m1c] ], [ [-self.pd3[n3][n2a][m1b][c] for c in range(3)], [ (-self.pd3[n3][n2c][m1c][c] + self.pd1[n3][n2c][m1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] + # [ self.px[n3][n2a][m1b], self.px[n3][n2c][m1c] ], [ [-self.pd3[n3][n2a][m1b][c] for c in range(3)], + # [ (-self.pd3[n3][n2c][m1c][c] + self.pd1[n3][n2c][m1c][c]) for c in range(3) ] ], + # 2, arcLengthDerivatives = True)[0:2] rtx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [ self.px[n3][n2c][m1a], self.px[n3][n2b][m1c] ], [ [ (-self.pd1[n3][n2c][m1a][c] - self.pd3[n3][n2c][m1a][c]) for c in range(3) ], [ -d for d in self.pd3[n3][n2b][m1c] ] ], 2, arcLengthDerivatives = True)[0:2] + [self.px[n3][n2c][m1a], self.px[n3][n2b][m1c]], + [[(-self.pd1[n3][n2c][m1a][c] - self.pd3[n3][n2c][m1a][c]) for c in range(3)], + [-d for d in self.pd3[n3][n2b][m1c]]], 2, arcLengthDerivatives=True)[0:2] rtx.append(tx[1]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2a][m1c], self.px[n3][n2c][m1b]], [[(self.pd1[n3][n2a][m1c][c] + self.pd2[n3][n2a][m1c][c]) for c in range(3)],self.pd2[n3][n2c][m1b]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2a][m1c], self.px[n3][n2c][m1b]], + [[(self.pd1[n3][n2a][m1c][c] + self.pd2[n3][n2a][m1c][c]) for c in range(3)], self.pd2[n3][n2c][m1b]], + 2, arcLengthDerivatives=True)[0: 2] rtx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2a][m1b], self.px[n3][n2c][m1c]], [self.pd2[n3][n2a][m1b], [(-self.pd1[n3][n2c][m1c][c] + self.pd2[n3][n2c][m1c][c]) for c in range(3)]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2a][m1b], self.px[n3][n2c][m1c]], + [self.pd2[n3][n2a][m1b], [(-self.pd1[n3][n2c][m1c][c] + self.pd2[n3][n2c][m1c][c]) for c in range(3)]], + 2, arcLengthDerivatives=True)[0: 2] rtx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( - [self.px[n3][n2c][m1a], self.px[n3][n2b][m1c]], [[(-self.pd1[n3][n2c][m1a][c] - self.pd2[n3][n2c][m1a][c]) for c in range(3)],[-d for d in self.pd1[n3][n2b][m1c]]], 2, arcLengthDerivatives = True)[0: 2] + [self.px[n3][n2c][m1a], self.px[n3][n2b][m1c]], + [[(-self.pd1[n3][n2c][m1a][c] - self.pd2[n3][n2c][m1a][c]) for c in range(3)], + [-d for d in self.pd1[n3][n2b][m1c]]], 2, arcLengthDerivatives=True)[0: 2] rtx.append(tx[1]) - #x = [ (rtx[0][c] + rtx[1][c] + rtx[2][c])/3.0 for c in range(3) ] - x = [ (rtx[0][c] + rtx[2][c])/2.0 for c in range(3) ] + # x = [ (rtx[0][c] + rtx[1][c] + rtx[2][c])/3.0 for c in range(3) ] + x = [(rtx[0][c] + rtx[2][c])/2.0 for c in range(3)] if self.trackSurface: - p = self.trackSurface.findNearestPosition(x, startPosition=self.trackSurface.createPositionProportion(*(self.pProportions[n2b][m1c]))) + p = self.trackSurface.findNearestPosition(x, startPosition=self.trackSurface.createPositionProportion(*( + self.pProportions[n2b][m1c]))) self.pProportions[n2b][m1b] = self.trackSurface.getProportion(p) x, sd1, sd2 = self.trackSurface.evaluateCoordinates(p, derivatives=True) d1, d2, d3 = calculate_surface_axes(sd1, sd2, vector.normalise(sd1)) self.pd3[n3][n2b][m1b] = d3 - self.px [n3][n2b][m1b] = x + self.px[n3][n2b][m1b] = x if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd3[n3][n2b][m1b] = [ (self.px[n3][n2b][m1b][c] - self.px[n3][n2b][m1c][c]) for c in range(3) ] - self.pd1[n3][n2b][m1b] = [ (self.px[n3][n2c][m1b][c] - self.px[n3][n2b][m1b][c]) for c in range(3) ] + self.pd3[n3][n2b][m1b] = [(self.px[n3][n2b][m1b][c] - self.px[n3][n2b][m1c][c]) for c in range(3)] + self.pd1[n3][n2b][m1b] = [(self.px[n3][n2c][m1b][c] - self.px[n3][n2b][m1b][c]) for c in range(3)] elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: self.pd1[n3][n2b][m1b] = [(self.px[n3][n2b][m1b][c] - self.px[n3][n2b][m1c][c]) for c in range(3)] self.pd2[n3][n2b][m1b] = [(self.px[n3][n2c][m1b][c] - self.px[n3][n2b][m1b][c]) for c in range(3)] if not self.trackSurface: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - self.pd2[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][m1b], self.pd1[n3][n2b][m1b])) + self.pd2[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd3[n3][n2b][m1b], + self.pd1[n3][n2b][m1b])) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: - self.pd3[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][m1b], self.pd2[n3][n2b][m1b])) - + self.pd3[n3][n2b][m1b] = vector.normalise(vector.crossproduct3(self.pd1[n3][n2b][m1b], + self.pd2[n3][n2b][m1b])) def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): - ''' + """ Smooth derivatives leading to triple points where 3 square elements merge. :param n3: Index of through-wall coordinates to use. - ''' + """ n1a = self.elementsCountRim n1b = n1a + 1 n1z = self.elementsCountAcross - self.elementsCountRim @@ -236,7 +286,7 @@ def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): if n2 < n2b: td3.append([-self.pd3[n3][n2][n1b][c] for c in range(3)]) else: - td3.append([(self.pd1[n3][n2b][n1b][c] + self.pd3[n3][n2b][n1b][c]) for c in range(3)] ) + td3.append([(self.pd1[n3][n2b][n1b][c] + self.pd3[n3][n2b][n1b][c]) for c in range(3)]) td3 = smoothCubicHermiteDerivativesLine(tx, td3, fixStartDirection=True, fixEndDirection=True) @@ -263,28 +313,32 @@ def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): tx = [] td2 = [] for n2 in range(0, n2c): - tx .append(self.px [n3][n2][n1b]) - td2.append(self.pd2[n3][n2][n1b] if (n2 < n2b) else [(self.pd1[n3][n2][n1b][c] + self.pd2[n3][n2][n1b][c]) for c in range(3)]) - td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixAllDirections=fixAllDirections, fixEndDerivative=True, magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + tx .append(self.px[n3][n2][n1b]) + td2.append(self.pd2[n3][n2][n1b] if (n2 < n2b) else + [(self.pd1[n3][n2][n1b][c] + self.pd2[n3][n2][n1b][c]) for c in range(3)]) + td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixAllDirections=fixAllDirections, fixEndDerivative=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) for n2 in range(0, n2b): self.pd2[n3][n2][n1b] = td2[n2] # right tx = [] td2 = [] for n2 in range(0, n2c): - tx .append(self.px [n3][n2][m1b]) - td2.append(self.pd2[n3][n2][m1b] if (n2 < n2b) else [(-self.pd1[n3][n2][m1b][c] + self.pd2[n3][n2][m1b][c]) for c in range(3)]) - td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixAllDirections=fixAllDirections,fixEndDerivative=True,magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) + tx .append(self.px[n3][n2][m1b]) + td2.append(self.pd2[n3][n2][m1b] if (n2 < n2b) else + [(-self.pd1[n3][n2][m1b][c] + self.pd2[n3][n2][m1b][c]) for c in range(3)]) + td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixAllDirections=fixAllDirections, fixEndDerivative=True, + magnitudeScalingMode=DerivativeScalingMode.HARMONIC_MEAN) for n2 in range(0, n2b): self.pd2[n3][n2][m1b] = td2[n2] def smoothDerivativesAroundRim(self, n3, n3d=None, rx=0): - ''' + """ Smooth derivatives around rim. :param n3: Index of through-wall coordinates to use. :param n3d: Which n3 index to copy initial derivatives from. If None, use n3. :param rx: rim index from 0 (around outside) to self.elementsCountRim - ''' + """ assert 0 <= rx <= self.elementsCountRim if not n3d: n3d = n3 @@ -298,7 +352,7 @@ def smoothDerivativesAroundRim(self, n3, n3d=None, rx=0): elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: if n2 > self.elementsCountRim: # regular rows if n1 <= self.elementsCountRim: - td1.append([ -d for d in self.pd2[n3d][n2][n1] ]) + td1.append([-d for d in self.pd2[n3d][n2][n1]]) else: td1.append(self.pd2[n3d][n2][n1]) else: @@ -316,19 +370,61 @@ def smoothDerivativesAroundRim(self, n3, n3d=None, rx=0): elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: if n2 > self.elementsCountRim: # regular rows if n1 <= self.elementsCountRim: - self.pd2[n3][n2][n1] = [ -d for d in td1[ix] ] + self.pd2[n3][n2][n1] = [-d for d in td1[ix]] else: self.pd2[n3][n2][n1] = td1[ix] else: self.pd1[n3][n2][n1] = td1[ix] + def remap_derivatives(self, squareMapping, circleMapping=None): + """ + It remaps the derivatives as indicated by squareMapping and circleMapping. + Limited to SHIELD_RIM_DERIVATIVE_MODE_AROUND. + :param squareMapping: List[up, right]. Determines what derivatives should be in the up and right directions in + the square part. Their values are in [1,2,3] range which 1, 2, 3 means d1, d2 and d3 respectively. + The negative sign reverses the direction. e.g. [-3,2] means d3 is down and d2 is right. The third derivative + is determined by the other two. RH rule applied. Assumes [1,3] initially. + :param circleMapping: List[circumferential, radial]. Determines what derivatives should be used for + circumferential and radial directions around the circle. + [-1, 3] means d1 -> clockwise and around. d3 -> outward and radial. Assumes [1,3] initially. + :return: + """ + dct = {1: self.pd1, 2: self.pd2, 3: self.pd3} + perm = {(1, 2): -3, (2, 3): -1, (3, 1): -2, (3, 2): 1, (1, 3): 2, (2, 1): 3} + + square = True + tripleRow = [self.elementsCountRim + 1, self.elementsCountUpFull - (self.elementsCountRim + 1)] + ellipseMapping = [squareMapping, circleMapping] + for mapping in ellipseMapping: + if mapping: + signs = [] + for c in mapping: + sign = 1 if c > 0 else -1 + signs.append(sign) + derv = (abs(mapping[0]), abs(mapping[1])) + sign = 1 if perm[derv] > 0 else -1 + signs.append(signs[0]*signs[1]*sign) + dervMapping = (derv[0], derv[1], abs(perm[derv])) + temp1 = copy.deepcopy(self.pd3) + temp2 = copy.deepcopy(self.pd2) + for n2 in range(self.elementsCountUpFull + 1): + for n3 in range(self.elementsCountAlong+1): + for n1 in range(self.elementsCountAcross + 1): + if self.px[n3][n2][n1]: + is_on_square = ((self.px[n3][n2][0] and self.px[n3][0][n1]) or n2 in tripleRow) + if (is_on_square and square) or (not is_on_square and not square): + dct[dervMapping[0]][n3][n2][n1] = [signs[0]*c for c in self.pd1[n3][n2][n1]] + dct[dervMapping[1]][n3][n2][n1] = [signs[1]*c for c in temp1[n3][n2][n1]] + dct[dervMapping[2]][n3][n2][n1] = [signs[2]*c for c in temp2[n3][n2][n1]] + square = False + def generateNodesForOtherHalf(self, mirrorPlane): """ Generates coordinates and derivatives for the other half by mirroring them. It keeps the d1 direction. :param mirrorPlane: plane ax+by+cz=d in form of [a,b,c,d] :return: """ - mirror=Mirror(mirrorPlane) + mirror = Mirror(mirrorPlane) for n2 in range(self.elementsCountUp): for n3 in range(self.elementsCountAlong + 1): for n1 in range(self.elementsCountAcross + 1): @@ -338,7 +434,7 @@ def generateNodesForOtherHalf(self, mirrorPlane): self.pd2[n3][2*self.elementsCountUp-n2][n1] = mirror.mirrorVector(self.pd2[n3][n2][n1]) self.pd3[n3][2*self.elementsCountUp-n2][n1] = mirror.mirrorVector(self.pd3[n3][n2][n1]) - def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier,mirrorPlane=None): + def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier, mirrorPlane=None): """ Create shield nodes from coordinates. :param fieldmodule: Zinc fieldmodule to create nodes in. Uses DOMAIN_TYPE_NODES. @@ -357,13 +453,13 @@ def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier,mirrorPlan nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) cache = fieldmodule.createFieldcache() - #for n2 in range(self.elementsCountUp, -1, -1): + # for n2 in range(self.elementsCountUp, -1, -1): # s = "" # for n1 in range(self.elementsCountAcross + 1): # s += str(n1) if self.px[1][n2][n1] else " " # print(n2, s, n2 - self.elementsCountUp - 1) - if self._mode == ShieldShape.SHIELD_SHAPE_FULL and mirrorPlane: + if self._mode == ShieldShape2D.SHIELD_SHAPE_FULL and mirrorPlane: self.generateNodesForOtherHalf(mirrorPlane) for n2 in range(self.elementsCountUpFull + 1): @@ -373,7 +469,7 @@ def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier,mirrorPlan node = nodes.createNode(nodeIdentifier, nodetemplate) self.nodeId[n3][n2][n1] = nodeIdentifier cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, self.px [n3][n2][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, self.px[n3][n2][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, self.pd1[n3][n2][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, self.pd2[n3][n2][n1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, self.pd3[n3][n2][n1]) @@ -381,7 +477,6 @@ def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier,mirrorPlan return nodeIdentifier - def generateElements(self, fieldmodule, coordinates, startElementIdentifier, meshGroups=[]): """ Create shield elements from nodes. @@ -421,22 +516,30 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes eft1 = eft scalefactors = None if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - nids = [ self.nodeId[e3][e2][e1], self.nodeId[e3][e2 + 1][e1], self.nodeId[e3+1][e2][e1], self.nodeId[e3+1][e2 + 1][e1], - self.nodeId[e3][e2][e1 + 1], self.nodeId[e3][e2 + 1][e1 + 1], self.nodeId[e3+1][e2][e1 + 1], self.nodeId[e3+1][e2 + 1][e1 + 1] ] + nids = [self.nodeId[e3][e2][e1], self.nodeId[e3][e2 + 1][e1], + self.nodeId[e3+1][e2][e1], self.nodeId[e3+1][e2 + 1][e1], + self.nodeId[e3][e2][e1 + 1], self.nodeId[e3][e2 + 1][e1 + 1], + self.nodeId[e3+1][e2][e1 + 1], self.nodeId[e3+1][e2 + 1][e1 + 1]] elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: - nids = [self.nodeId[0][e2][e1], self.nodeId[0][e2][e1 + 1], self.nodeId[0][e2 + 1][e1],self.nodeId[0][e2 + 1][e1 + 1], - self.nodeId[1][e2][e1], self.nodeId[1][e2][e1 + 1], self.nodeId[1][e2 + 1][e1], self.nodeId[1][e2 + 1][e1 + 1]] + nids = [self.nodeId[0][e2][e1], self.nodeId[0][e2][e1 + 1], + self.nodeId[0][e2 + 1][e1], self.nodeId[0][e2 + 1][e1 + 1], + self.nodeId[1][e2][e1], self.nodeId[1][e2][e1 + 1], + self.nodeId[1][e2 + 1][e1], self.nodeId[1][e2 + 1][e1 + 1]] if (e2 < e2b) or (e2 > e2y): if (e1 < e1b) or (e1 > e1y): continue # no element due to triple point closure if (e2 < e2a) or (e2 > e2z): if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e2 < e2a: - nids = [self.nodeId[e3][e2+1][e1], self.nodeId[e3][e2+1][e1+1], self.nodeId[e3+1][e2+1][e1], self.nodeId[e3+1][e2+1][e1+1], - self.nodeId[e3][e2][e1], self.nodeId[e3][e2][e1+1], self.nodeId[e3+1][e2][e1], self.nodeId[e3+1][e2][e1+1]] + nids = [self.nodeId[e3][e2+1][e1], self.nodeId[e3][e2+1][e1+1], + self.nodeId[e3+1][e2+1][e1], self.nodeId[e3+1][e2+1][e1+1], + self.nodeId[e3][e2][e1], self.nodeId[e3][e2][e1+1], + self.nodeId[e3+1][e2][e1], self.nodeId[e3+1][e2][e1+1]] elif e2 > e2z: - nids = [self.nodeId[e3][e2][e1+1], self.nodeId[e3][e2][e1], self.nodeId[e3+1][e2][e1+1], self.nodeId[e3+1][e2][e1], - self.nodeId[e3][e2+1][e1+1], self.nodeId[e3][e2+1][e1], self.nodeId[e3+1][e2+1][e1+1], self.nodeId[e3+1][e2+1][e1]] + nids = [self.nodeId[e3][e2][e1+1], self.nodeId[e3][e2][e1], + self.nodeId[e3+1][e2][e1+1], self.nodeId[e3+1][e2][e1], + self.nodeId[e3][e2+1][e1+1], self.nodeId[e3][e2+1][e1], + self.nodeId[e3+1][e2+1][e1+1], self.nodeId[e3+1][e2+1][e1]] elif (e2 == e2a) or (e2 == e2z): # bottom and top row elements if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: @@ -444,36 +547,52 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] - remapEftNodeValueLabel(eft1, [1, 3, 5, 7], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS3, [1])]) - remapEftNodeValueLabel(eft1, [1, 3, 5, 7], Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft1, [1, 3, 5, 7], + Node.VALUE_LABEL_D_DS1, [(Node.VALUE_LABEL_D_DS3, [1])]) + remapEftNodeValueLabel(eft1, [1, 3, 5, 7], + Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) if (e1 == e1b) or (e1 == e1y): # map bottom triple point element if e1 == e1b: - remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [2, 4], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [])]) else: - remapEftNodeValueLabel(eft1, [ 6, 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, [6, 8], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [1])]) elif e2 == e2z: eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] - remapEftNodeValueLabel(eft1, [2, 4, 6, 8], Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft1, [2, 4, 6, 8], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft1, [2, 4, 6, 8], Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft1, [2, 4, 6, 8], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS3, [])]) if (e1 == e1b) or (e1 == e1y): # map top triple point element if e1 == e1b: - remapEftNodeValueLabel(eft1, [1, 3], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS3, [1])]) + remapEftNodeValueLabel(eft1, [1, 3], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [1])]) else: - remapEftNodeValueLabel(eft1, [5, 7], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS1, []),(Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft1, [5, 7], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [])]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: if (e1 == e1b) or (e1 == e1y): # map bottom triple point element eft1 = tricubichermite.createEftNoCrossDerivatives() if e1 == e1b: - remapEftNodeValueLabel(eft1, [3, 7], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, []),(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft1, [3, 7], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS2, [])]) else: setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] - remapEftNodeValueLabel(eft1, [4, 8], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, [1]),(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft1, [4, 8], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1]), + (Node.VALUE_LABEL_D_DS2, [])]) elif (e2 == e2b) or (e2 == e2y): if (e1 <= e1a) or (e1 >= e1z): @@ -481,51 +600,64 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes e2r = e1 if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e2 == e2b: - nids = [self.nodeId[e3][e2c][e1+1], self.nodeId[e3][e2r+1][e1b], self.nodeId[e3+1][e2c][e1+1], self.nodeId[e3+1][e2r+1][e1b], - self.nodeId[e3][e2c][e1], self.nodeId[e3][e2r][e1b], self.nodeId[e3+1][e2c][e1], self.nodeId[e3+1][e2r][e1b]] + nids = [self.nodeId[e3][e2c][e1+1], self.nodeId[e3][e2r+1][e1b], + self.nodeId[e3+1][e2c][e1+1], self.nodeId[e3+1][e2r+1][e1b], + self.nodeId[e3][e2c][e1], self.nodeId[e3][e2r][e1b], + self.nodeId[e3+1][e2c][e1], self.nodeId[e3+1][e2r][e1b]] if e2 == e2y: e2r = 2*self.elementsCountUp - e1-1 - nids = [self.nodeId[e3][e2r][e1b], self.nodeId[e3][e2y][e1+1], self.nodeId[e3+1][e2r][e1b], self.nodeId[e3+1][e2y][e1+1], - self.nodeId[e3][e2r+1][e1b], self.nodeId[e3][e2y][e1], self.nodeId[e3+1][e2r+1][e1b], self.nodeId[e3+1][e2y][e1]] + nids = [self.nodeId[e3][e2r][e1b], self.nodeId[e3][e2y][e1+1], + self.nodeId[e3+1][e2r][e1b], self.nodeId[e3+1][e2y][e1+1], + self.nodeId[e3][e2r+1][e1b], self.nodeId[e3][e2y][e1], + self.nodeId[e3+1][e2r+1][e1b], self.nodeId[e3+1][e2y][e1]] elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] - nids[0] = self.nodeId[0][e2r ][e1b] + nids[0] = self.nodeId[0][e2r][e1b] nids[1] = self.nodeId[0][e2r + 1][e1b] - nids[4] = self.nodeId[1][e2r ][e1b] + nids[4] = self.nodeId[1][e2r][e1b] nids[5] = self.nodeId[1][e2r + 1][e1b] - setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) + remapEftNodeValueLabel(eft1, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft1, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS2, [])]) elif e1 == e1a: eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] + scalefactors = [-1.0] if e2 == e2b: nids[0] = self.nodeId[e3][e2a][e1b] nids[2] = self.nodeId[e3+1][e2a][e1b] tripleN = [5, 7] - remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [])]) elif e2 == e2y: nids[1] = self.nodeId[e3][e2z+1][e1b] nids[3] = self.nodeId[e3+1][e2z+1][e1b] tripleN = [6, 8] - remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, [1]), (Node.VALUE_LABEL_D_DS3, [])]) - remapEftNodeValueLabel(eft1, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) - remapEftNodeValueLabel(eft1, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS3, [1] ) ]) + remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS1, [1]), + (Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS3, [1])]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] + scalefactors = [-1.0] nids[0] = self.nodeId[0][e2a][e1b] nids[4] = self.nodeId[1][e2a][e1b] - remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS2, [])]) - remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS1, []),(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS2, [])]) elif e1 == e1z: eft1 = tricubichermite.createEftNoCrossDerivatives() if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: @@ -534,27 +666,35 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes scalefactors = [-1.0] nids[4] = self.nodeId[e3][e2a][e1z] nids[6] = self.nodeId[e3+1][e2a][e1z] - setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] - remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) + remapEftNodeValueLabel(eft1, [1, 3], Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS1, [1]), + (Node.VALUE_LABEL_D_DS3, [])]) elif e2 == e2y: nids[5] = self.nodeId[e3][e2z+1][e1z] nids[7] = self.nodeId[e3+1][e2z+1][e1z] - remapEftNodeValueLabel(eft1, [2, 4], Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft1, [2, 4], Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS3, [])]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] nids[1] = self.nodeId[0][e2a][e1z] nids[5] = self.nodeId[1][e2a][e1z] - remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS1, []),(Node.VALUE_LABEL_D_DS2, [1])]) - remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS1,[(Node.VALUE_LABEL_D_DS2, [1])]) - remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft1, [1, 5], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, []), + (Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [])]) elif e1 > e1z: e2r = self.elementsCountAcross - e1 if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e2 == e2b: - nids = [self.nodeId[e3][e2r][e1z], self.nodeId[e3][e2c][e1], self.nodeId[e3+1][e2r][e1z], self.nodeId[e3+1][e2c][e1], - self.nodeId[e3][e2r-1][e1z], self.nodeId[e3][e2c][e1+1], self.nodeId[e3+1][e2r-1][e1z], self.nodeId[e3+1][e2c][e1+1]] + nids = [self.nodeId[e3][e2r][e1z], self.nodeId[e3][e2c][e1], + self.nodeId[e3+1][e2r][e1z], self.nodeId[e3+1][e2c][e1], + self.nodeId[e3][e2r-1][e1z], self.nodeId[e3][e2c][e1+1], + self.nodeId[e3+1][e2r-1][e1z], self.nodeId[e3+1][e2c][e1+1]] elif e2 == e2y: e2r = e2z+e1-e1z nids[1] = self.nodeId[e3][e2r][e1z] @@ -565,26 +705,30 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] - nids[0] = self.nodeId[0][e2r ][e1z] + nids[0] = self.nodeId[0][e2r][e1z] nids[1] = self.nodeId[0][e2r - 1][e1z] - nids[4] = self.nodeId[1][e2r ][e1z] + nids[4] = self.nodeId[1][e2r][e1z] nids[5] = self.nodeId[1][e2r - 1][e1z] - setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ) ]) + remapEftNodeValueLabel(eft1, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft1, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [])]) else: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e1 < e1a: - nids = [ self.nodeId[e3][e2 + 1][e1 + 1], self.nodeId[e3][e2][e1 + 1], self.nodeId[e3+1][e2 + 1][e1 + 1], self.nodeId[e3+1][e2][e1 + 1], - self.nodeId[e3][e2 + 1][e1], self.nodeId[e3][e2][e1], self.nodeId[e3+1][e2 + 1][e1], self.nodeId[e3+1][e2][e1]] + nids = [self.nodeId[e3][e2 + 1][e1 + 1], self.nodeId[e3][e2][e1 + 1], + self.nodeId[e3+1][e2 + 1][e1 + 1], self.nodeId[e3+1][e2][e1 + 1], + self.nodeId[e3][e2 + 1][e1], self.nodeId[e3][e2][e1], + self.nodeId[e3+1][e2 + 1][e1], self.nodeId[e3+1][e2][e1]] elif e1 == e1a: # map left column elements eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) - scalefactors = [ -1.0 ] - remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS1, [(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS3, [1])]) + scalefactors = [-1.0] + remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS1, + [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft1, [1, 2, 3, 4], Node.VALUE_LABEL_D_DS3, + [(Node.VALUE_LABEL_D_DS3, [1])]) if eft1 is not eft: elementtemplate1.defineField(coordinates, -1, eft1) @@ -596,7 +740,7 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes result3 = element.setScaleFactors(eft1, scalefactors) else: result3 = 7 - #print('create element shield', elementIdentifier, result2, result3, nids) + # print('create element shield', elementIdentifier, result2, result3, nids) if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: self.elementId[e3][e2][e1] = elementIdentifier else: @@ -607,3 +751,1439 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes meshGroup.addElement(element) return elementIdentifier + + +class ShieldMesh3D: + """ + Generates a 3D shield mesh. + """ + + def __init__(self, elementsCountAcross, elementsCountRim, box_derivatives=None): + """ + 3D shield structure can be used for a sphere mesh. 3 hex mesh merges to one box in the corner located in the + centre. + The structure is a 3D version of the 2D shield structure. It has a 'quadruple point', a unique node on the + surface where the elements merge. 'Triple curves' meet at the quadruple point and a fourth curve that connects + to another quadruple point which is inside. + The structure is like a elementsCountAcross[0] X elementsCountAcross[1] X elementsCountAcross[2] box where + some nodes does not exist and stored as None. The quadruple point is stored at [n3z][0][n1z] (top, front and + right most node) where n3z is top most and n1z is the right most indexes. + Triple curves and surfaces connecting to the quadruple point divide the exterior surface and interior region + into 3 regions (top, left and right). Triple curve 1 connects the quadruple point to the plane 2-3. + Similarly, triple curves 2 and 3 connect the quadruple point to the planes 1-3 and 1-2, respectively. + It is the same for the inside quadruple. Any point on region left has n2 = 0. Similarly on region + right -> n1 = n1z and on top -> n3 = n3z. + There is a gap between node indexes of a node on a triple curve and the next nodes on the surface. + + :param elementsCountAcross: number of elements as + a list [elementsCountAcrossAxis1, elementsCountAcrossAxis2, elementsCountAcrossAxis3] + :param elementsCountRim: + :param box_derivatives: It is a list of [deriv1,deriv2,deriv3]. It is used to set the derivative directions. + default is [1, 3, 2] which means it makes -d1, d3 and d2 in direction of axis1, axis2 and axis3. To make + d1, d2 and d3 directions in direction of axis1, axis2 and axis3 use [-1, 2, 3]. + """ + assert elementsCountRim >= 0 + # assert elementsCountAcross >= (elementsCountRim + 4) + # assert elementsCountUpFull >= (elementsCountRim + 2) + self.elementsCountAcross = elementsCountAcross + self.elementsCountRim = elementsCountRim + # self.elementsCountUpRegular = elementsCountUp - 2 - elementsCountRim + # elementsCountAcrossNonRim = self.elementsCountAcross - 2*elementsCountRim + # self.elementsCountAroundFull = 2*self.elementsCountUpRegular + elementsCountAcrossNonRim + self._boxDerivatives = box_derivatives + self._boxMapping = None + self._element_needs_scale_factor = False + self._xi_mapping = None + self._xi_signs = None + self._box_deriv_mapping = None + self._box_deriv_signs = None + + self.px = [[] for _ in range(elementsCountAcross[2] + 1)] + self.pd1 = [[] for _ in range(elementsCountAcross[2] + 1)] + self.pd2 = [[] for _ in range(elementsCountAcross[2] + 1)] + self.pd3 = [[] for _ in range(elementsCountAcross[2] + 1)] + self.nodeId = [[] for _ in range(elementsCountAcross[2] + 1)] + for n3 in range(elementsCountAcross[2] + 1): + for n2 in range(elementsCountAcross[0] + 1): + for p in [self.px[n3], self.pd1[n3], self.pd2[n3], self.pd3[n3], self.nodeId[n3]]: + p.append([None]*(elementsCountAcross[1] + 1)) + + self.elementId = [[[None]*elementsCountAcross[1] for n2 in range(elementsCountAcross[0])] + for e3 in range(elementsCountAcross[2])] + + # node types + CORNER_1 = 1 + CORNER_2 = 2 + CORNER_3 = 3 + + QUADRUPLE_DOWN_LEFT = 4 + QUADRUPLE_RIGHT = 5 + QUADRUPLE_UP = 6 + QUADRUPLE0_DOWN_LEFT = 7 + QUADRUPLE0_RIGHT = 8 + QUADRUPLE0_UP = 9 + + TRIPLE_12_LEFT = 10 + TRIPLE_12_RIGHT = 11 + TRIPLE_13_DOWN = 12 + TRIPLE_13_UP = 13 + TRIPLE_23_UP = 14 + TRIPLE_23_DOWN = 15 + TRIPLE0_12_LEFT = 16 + TRIPLE0_12_RIGHT = 17 + TRIPLE0_13_DOWN = 18 + TRIPLE0_13_Up = 19 + TRIPLE0_23_DOWN = 20 + TRIPLE0_23_UP = 21 + + BOUNDARY_12_LEFT = 22 + BOUNDARY_12_RIGHT = 23 + BOUNDARY_13_DOWN = 24 + BOUNDARY_13_UP = 25 + BOUNDARY_23_UP = 26 + BOUNDARY_23_DOWN = 27 + + TRIPLE_CURVE_1_DOWN = 28 + TRIPLE_CURVE_1_UP = 29 + TRIPLE_CURVE_2_DOWN = 30 + TRIPLE_CURVE_2_UP = 31 + TRIPLE_CURVE_3_LEFT = 32 + TRIPLE_CURVE_3_RIGHT = 33 + TRIPLE_CURVE0_1_UP = 34 + TRIPLE_CURVE0_1_DOWN = 35 + TRIPLE_CURVE0_2_DOWN = 36 + TRIPLE_CURVE0_2_UP = 37 + TRIPLE_CURVE0_3_LEFT = 38 + TRIPLE_CURVE0_3_RIGHT = 39 + + SURFACE_REGULAR_DOWN_LEFT = 40 + SURFACE_REGULAR_DOWN_RIGHT = 41 + SURFACE_REGULAR_UP = 42 + REGULAR = 43 + + # element types + ELEMENT_REGULAR = 1 + ELEMENT_QUADRUPLE_DOWN_LEFT = 2 + ELEMENT_QUADRUPLE_DOWN_RIGHT = 3 + ELEMENT_QUADRUPLE_UP_LEFT = 4 + ELEMENT_QUADRUPLE_DOWN = 5 + ELEMENT_QUADRUPLE_UP = 6 + ELEMENT_QUADRUPLE_LEFT = 7 + ELEMENT_QUADRUPLE_RIGHT = 8 + ELEMENT_DOWN_RIGHT = 9 + ELEMENT_DOWN_LEFT = 10 + + def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier, rangeOfRequiredElements): + """ + Create shield nodes from coordinates. + :param fieldmodule: Zinc fieldmodule to create nodes in. Uses DOMAIN_TYPE_NODES. + :param coordinates: Coordinate field to define. + :param startNodeIdentifier: First node identifier to use. + :param rangeOfRequiredElements: Only the elements and nodes for the given ragne is generated. + :return: next nodeIdentifier. + """ + nodeIdentifier = startNodeIdentifier + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + # In case we swap derivatives in the when we use remap function. + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + cache = fieldmodule.createFieldcache() + + for n2 in range(self.elementsCountAcross[0] + 1): + for n3 in range(self.elementsCountAcross[2] + 1): + for n1 in range(self.elementsCountAcross[1] + 1): + if n3 > rangeOfRequiredElements[2][1] or n3 < rangeOfRequiredElements[2][0]\ + or n2 > rangeOfRequiredElements[0][1] or n2 < rangeOfRequiredElements[0][0]\ + or n1 > rangeOfRequiredElements[1][1] or n1 < rangeOfRequiredElements[1][0]: + continue + if self.px[n3][n2][n1]: + node = nodes.createNode(nodeIdentifier, nodetemplate) + self.nodeId[n3][n2][n1] = nodeIdentifier + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, self.px[n3][n2][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, self.pd1[n3][n2][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, self.pd2[n3][n2][n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, self.pd3[n3][n2][n1]) + nodeIdentifier += 1 + + return nodeIdentifier + + def set_derivatives(self, box_derivatives): + """ + Set the derivatives as specified by box_derivatives. + :param box_derivatives: List[d_axis1_negative, d_axis2, d_axis3]. Determines what derivatives should be + for back, right and up directions in the box region. Their values are in [1,2,3] range which 1, 2, 3 means + d1, d2 and d3 respectively. + The negative sign reverses the direction. e.g. [1, -3,2] means for the right direction use -d3. + Default is [1,3, 2]. + :return: + """ + self._boxMapping = box_derivatives + + dct = {1: self.pd1, 2: self.pd2, 3: self.pd3} + perm = {(1, 2): -3, (2, 3): -1, (3, 1): -2, (3, 2): 1, (1, 3): 2, (2, 1): 3} + + if box_derivatives: + signs = [1 if c > 0 else -1 for c in box_derivatives] + dervMapping = (abs(box_derivatives[0]), abs(box_derivatives[1]), abs(box_derivatives[2])) + + temp1 = copy.deepcopy(self.pd3) + temp2 = copy.deepcopy(self.pd2) + for n3 in range(self.elementsCountAcross[2] + 1): + for n2 in range(self.elementsCountAcross[0] + 1): + for n1 in range(self.elementsCountAcross[1] + 1): + if self.px[n3][n2][n1]: + if self.is_interior_regular_nodes(n3, n2, n1): + dct[dervMapping[0]][n3][n2][n1] = [signs[0] * c for c in self.pd1[n3][n2][n1]] + dct[dervMapping[1]][n3][n2][n1] = [signs[1] * c for c in temp1[n3][n2][n1]] + dct[dervMapping[2]][n3][n2][n1] = [signs[2] * c for c in temp2[n3][n2][n1]] + + def is_interior_regular_nodes(self, n3, n2, n1): + """ + Determine if a node indicated by [n3,n2,n1], is inside the sphere. + :return: True, if the node is inside. + """ + n3z = self.elementsCountAcross[2] + n1z = self.elementsCountAcross[1] + n2z = self.elementsCountAcross[0] + n3y = n3z - 1 + n1y = n1z - 1 + + return n3 <= n3y and n2 >= 1 and n1 <= n1y + + def set_derivatives_for_irregular_nodes(self, octant_number): + """ + For irregular nodes such as corners that are common between octants, change local derivatives to be the same + directions of the octant_PPP. This is because the remapping is done for octant_PPP only and for others we need + to make them similar, so similar remapping can be performed. Therefore, always same vectors is used for + remapping. + e.g. [1, 0, 0] generates curve 1 in octant1. For one octant d1 is [1,0,0] and for another octant -d2. + :param octant_number: see get_element_type + :return: Derivatives for irregular nodes. + """ + default = [1, 2, 3] + + corner1derivs = default + corner2derivs = default + corner3derivs = default + boundary12leftderivs = default + boundary12rightderivs = default + triple12leftderivs = default + triple12rightderivs = default + + if octant_number == 1: + corner3derivs = [-2, 1, 3] + elif octant_number == 2: + corner3derivs = [-1, -2, 3] + elif octant_number == 3: + corner3derivs = [2, -1, 3] + elif octant_number == 4: + corner3derivs = [1, 2, 3] + elif octant_number == 5: + corner3derivs = [-1, -2, 3] + elif octant_number == 6: + corner3derivs = [-2, 1, 3] + elif octant_number == 7: + corner3derivs = [1, 2, 3] + elif octant_number == 8: + corner3derivs = [2, -1, 3] + else: + corner3derivs = [1, 2, 3] + + if octant_number in [5, 6, 7, 8]: + corner1derivs = [-1, -2, 3] + corner2derivs = [-1, -2, 3] + boundary12leftderivs = [-1, -2, 3] + boundary12rightderivs = [-1, -2, 3] + triple12leftderivs = [-1, -2, 3] + triple12rightderivs = [-1, -2, 3] + + return corner1derivs, corner2derivs, corner3derivs, boundary12leftderivs, boundary12rightderivs, \ + triple12leftderivs, triple12rightderivs + + def get_box_mapping_for_other_octants(self, octant_number, octant1_box_derivatives=None): + """ + Find mapping for box nodes derivatives to make each octant similar to octant1. + :param octant1_box_derivatives: default is [1, 3, 2] for octant 1 box derivatives. + :return: boxMapping + """ + if not octant1_box_derivatives: + octant1_box_derivatives = [1, 3, 2] + + swap12 = [octant1_box_derivatives[1], octant1_box_derivatives[0], octant1_box_derivatives[2]] + + signs = self.get_octant_signs(octant_number) + + if octant_number in [2, 4, 5, 7]: + boxMapping = [swap12[c] * signs[c] for c in range(3)] + else: + boxMapping = [octant1_box_derivatives[c] * signs[c] for c in range(3)] + + return boxMapping + + def get_octant_signs(self, octant_number): + """ + Box mapping for each octant is different from octant1. for some of them 1 and 2 swaps. Also, the sign of them + might change. + :return: + """ + if octant_number == 1: + signs = [1, 1, 1] + elif octant_number == 2: + # y<0 + signs = [1, -1, 1] + elif octant_number == 3: + # x<0, y<0 + signs = [-1, -1, 1] + elif octant_number == 4: + # x<0 + signs = [-1, 1, 1] + elif octant_number == 5: + # z<0 + signs = [-1, -1, -1] + elif octant_number == 6: + # y<0, z<0 + signs = [1, -1, -1] + elif octant_number == 7: + # x<0, y<0, z<0 + signs = [1, 1, -1] + elif octant_number == 8: + # x<0, z<0 + signs = [-1, 1, -1] + else: + signs = [1, 1, 1] + + return signs + + def get_element_node_identifiers(self, octant_number, e3, e2, e1, e3zo, e1zo): + """ + Find node identifiers for given element represented by its indexes (e3,e2,e1). It uses default order of nodes + for [1,3,2] default element axes and finds the mapping between octant default local node numbers and + the sphere. + :param octant_number: + :return: nids, node s for given element represented by its indexes (e3,e2,e1). + """ + nids_default = [self.getNodeId(octant_number, e3, e2, e1, e3zo, e1zo), + self.getNodeId(octant_number, e3, e2+1, e1, e3zo, e1zo), + self.getNodeId(octant_number, e3+1, e2, e1, e3zo, e1zo), + self.getNodeId(octant_number, e3+1, e2+1, e1, e3zo, e1zo), + self.getNodeId(octant_number, e3, e2, e1+1, e3zo, e1zo), + self.getNodeId(octant_number, e3, e2+1, e1+1, e3zo, e1zo), + self.getNodeId(octant_number, e3+1, e2, e1+1, e3zo, e1zo), + self.getNodeId(octant_number, e3+1, e2+1, e1+1, e3zo, e1zo)] + + # change the order of nodes according to the box derivatives. + boxMappingOctant1 = self.get_box_mapping_for_other_octants(1, octant1_box_derivatives=self._boxDerivatives) + lnm_octant1 = self.local_node_mapping(boxMappingOctant1) + nids = [1] * 8 + for ln in range(1, 9): + nids[lnm_octant1[ln] - 1] = nids_default[ln - 1] + + return nids + + def generateElements(self, fieldmodule, coordinates, startElementIdentifier, rangeOfRequiredElements, + meshGroups=[]): + """ + Create shield elements from nodes. + :param fieldmodule: Zinc fieldmodule to create elements in. + :param coordinates: Coordinate field to define. + :param startElementIdentifier: First element identifier to use. + :param rangeOfRequiredElements: Only the elements and nodes for the given range is generated. + :param meshGroups: Zinc mesh groups to add elements to. + :return: next elementIdentifier. + """ + elementIdentifier = startElementIdentifier + useCrossDerivatives = False + mesh = fieldmodule.findMeshByDimension(3) + + tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives) + eft = tricubichermite.createEftNoCrossDerivatives() + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + + elementtemplate1 = mesh.createElementtemplate() + elementtemplate1.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + for e3 in range(self.elementsCountAcross[2]): + for e2 in range(self.elementsCountAcross[0]): + for e1 in range(self.elementsCountAcross[1]): + if e3 >= rangeOfRequiredElements[2][1] or e3 < rangeOfRequiredElements[2][0] or\ + e2 >= rangeOfRequiredElements[0][1] or e2 < rangeOfRequiredElements[0][0]\ + or e1 >= rangeOfRequiredElements[1][1] or e1 < rangeOfRequiredElements[1][0]: + continue + + scalefactors = None + self._element_needs_scale_factor = False + + octant_number, element_type, e3zo, e2zo, e1zo = self.get_element_type(e3, e2, e1) + e3o, e2o, e1o = self.get_local_element_index(octant_number, e3, e2, e1) + e3yo, e2bo, e1yo = e3zo - 1, 1, e1zo - 1 + + eft1 = eft if element_type == self.ELEMENT_REGULAR else\ + tricubichermite.createEftNoCrossDerivatives() + + nids = self.get_element_node_identifiers(octant_number, e3, e2, e1, e3zo, e1zo) + + boxMapping = self.get_box_mapping_for_other_octants(octant_number, + octant1_box_derivatives=self._boxDerivatives) + lnm = self.local_node_mapping(boxMapping) + corner1derivs, corner2derivs, corner3derivs, boundary12leftderivs, boundary12rightderivs,\ + triple12leftderivs, triple12rightderivs =\ + self.set_derivatives_for_irregular_nodes(octant_number) + + if element_type == self.ELEMENT_REGULAR: + pass + elif element_type == self.ELEMENT_QUADRUPLE_DOWN_LEFT: + self.remap_eft_node_value_label(eft1, [lnm[7]], self.QUADRUPLE_DOWN_LEFT) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.QUADRUPLE0_DOWN_LEFT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_12_LEFT, + derivatives=triple12leftderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE0_12_LEFT) + else: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_3_LEFT) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE0_3_LEFT) + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE0_13_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.CORNER_1, + derivatives=corner1derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_13_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE0_2_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_DOWN) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_DOWN_LEFT) + + elif element_type == self.ELEMENT_QUADRUPLE_DOWN: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE0_13_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_2_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE0_2_DOWN) + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.CORNER_1, + derivatives=corner1derivs) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_13_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.SURFACE_REGULAR_DOWN_LEFT) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[5]], self.SURFACE_REGULAR_DOWN_LEFT) + + elif element_type == self.ELEMENT_QUADRUPLE_DOWN_RIGHT: + if e2o == e2bo: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.QUADRUPLE0_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.QUADRUPLE_RIGHT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE0_12_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_12_RIGHT, + derivatives=triple12rightderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE0_3_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_3_RIGHT) + + if e2bo == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE0_23_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_23_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, + derivatives=corner2derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_23_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE0_1_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE_1_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.SURFACE_REGULAR_DOWN_RIGHT) + + elif e2o == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE0_1_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE0_23_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_1_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_23_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, + derivatives=corner2derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_23_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.TRIPLE_CURVE0_1_DOWN) + if e3yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[7], lnm[8]], self.TRIPLE_CURVE_1_DOWN) + + elif element_type == self.ELEMENT_QUADRUPLE_UP_LEFT: + if e2o == e2bo: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE0_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.QUADRUPLE0_UP) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.QUADRUPLE_UP) + if e2bo == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE0_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_23_UP) + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.CORNER_3, + derivatives=corner3derivs) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE0_1_UP) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE_1_UP) + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_13_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.SURFACE_REGULAR_UP) + elif e2o == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE0_1_UP) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE0_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_1_UP) + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.CORNER_3, + derivatives=corner3derivs) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.BOUNDARY_13_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_UP) + else: + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.TRIPLE_CURVE0_1_UP) + self.remap_eft_node_value_label(eft1, [lnm[7], lnm[8]], self.TRIPLE_CURVE_1_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.TRIPLE_CURVE0_1_UP) + self.remap_eft_node_value_label(eft1, [lnm[7], lnm[8]], self.TRIPLE_CURVE_1_UP) + + elif element_type == self.ELEMENT_QUADRUPLE_UP: + if e2o == e2bo: + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE0_13_Up) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE0_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE0_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_2_UP) + if e2bo == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.BOUNDARY_23_UP) + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.CORNER_3, + derivatives=corner3derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_23_UP) + + else: + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_13_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.SURFACE_REGULAR_UP) + elif e2o == e2zo: + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.CORNER_3, + derivatives=corner3derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.BOUNDARY_23_UP) + else: + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[7], lnm[8]], self.SURFACE_REGULAR_UP) + else: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4], lnm[7], lnm[8]], + self.SURFACE_REGULAR_UP) + + elif element_type == self.ELEMENT_QUADRUPLE_LEFT: + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_3_LEFT) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE0_3_LEFT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_12_LEFT, + derivatives=triple12leftderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE0_12_LEFT) + else: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_3_LEFT) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE0_3_LEFT) + if e1yo == 0: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.BOUNDARY_13_DOWN) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.CORNER_1, + derivatives=corner1derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_13_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_LEFT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_DOWN_LEFT) + + elif element_type == self.ELEMENT_QUADRUPLE_RIGHT: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE0_3_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.TRIPLE_CURVE_3_RIGHT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE0_12_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_12_RIGHT, + derivatives=triple12rightderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE0_3_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_3_RIGHT) + if e2bo == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.BOUNDARY_23_DOWN) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, + derivatives=corner2derivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_23_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.SURFACE_REGULAR_DOWN_RIGHT) + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + else: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.SURFACE_REGULAR_DOWN_RIGHT) + + elif element_type == self.ELEMENT_DOWN_RIGHT: + if e3o == 0: + if e2o == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[7]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, + derivatives=corner2derivs) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.BOUNDARY_23_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[7]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[8]], self.SURFACE_REGULAR_DOWN_RIGHT) + else: + self.remap_eft_node_value_label(eft1, [lnm[7]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[5]], self.SURFACE_REGULAR_DOWN_RIGHT) + if e2o == e2zo: + self.remap_eft_node_value_label(eft1, [lnm[6], lnm[8]], self.BOUNDARY_23_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[6], lnm[8]], self.SURFACE_REGULAR_DOWN_RIGHT) + + elif element_type == self.ELEMENT_DOWN_LEFT: + if e3o == 0: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.BOUNDARY_12_LEFT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[7]], self.SURFACE_REGULAR_DOWN_LEFT) + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.CORNER_1, + derivatives=corner1derivs) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.BOUNDARY_13_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_LEFT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_LEFT) + else: + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[7]], self.SURFACE_REGULAR_DOWN_LEFT) + if e1o == 0: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[3]], self.BOUNDARY_13_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[3]], self.SURFACE_REGULAR_DOWN_LEFT) + + else: + continue + + if self._element_needs_scale_factor: + scalefactors = [-1.0] + + if eft1 is not eft: + elementtemplate1.defineField(coordinates, -1, eft1) + element = mesh.createElement(elementIdentifier, elementtemplate1) + else: + element = mesh.createElement(elementIdentifier, elementtemplate) + result2 = element.setNodesByIdentifier(eft1, nids) + if scalefactors: + result3 = element.setScaleFactors(eft1, scalefactors) + else: + result3 = 7 + + self.elementId[e3][e2][e1] = elementIdentifier + + elementIdentifier += 1 + + if element_type == self.ELEMENT_REGULAR: + for meshGroup in meshGroups[:1]: + meshGroup.addElement(element) + else: + for meshGroup in meshGroups[1:2]: + meshGroup.addElement(element) + + return elementIdentifier + + def local_node_mapping(self, boxMapping): + """ + Find local node mapping between default local node number of octant and global local node number of sphere. + Obtain xi mapping and box irregular derivatives mapping. they are used for elements in octant that locally + might have different xi directions from their global sphere xi directions. Then they are used in + remap_eft_node_value_label. + :return: lnm, local node mapping between local node number of octant and local node number of sphere. + """ + boxMapping_default = [1, 3, 2] + + xi_node_map = {(0, 0, 0): 1, (1, 0, 0): 2, (0, 1, 0): 3, (1, 1, 0): 4, + (0, 0, 1): 5, (1, 0, 1): 6, (0, 1, 1): 7, (1, 1, 1): 8} + node_xi_map = {1: (0, 0, 0), 2: (1, 0, 0), 3: (0, 1, 0), 4: (1, 1, 0), + 5: (0, 0, 1), 6: (1, 0, 1), 7: (0, 1, 1), 8: (1, 1, 1)} + + xi_dm = {1: 'xi1', 2: 'xi2', 3: 'xi3'} + xi_default_to_global_map = {xi_dm[abs(boxMapping_default[0])]: abs(boxMapping[0]), + xi_dm[abs(boxMapping_default[1])]: abs(boxMapping[1]), + xi_dm[abs(boxMapping_default[2])]: abs(boxMapping[2])} + xi_global_to_default_map = {xi_dm[abs(boxMapping[0])]: abs(boxMapping_default[0]), + xi_dm[abs(boxMapping[1])]: abs(boxMapping_default[1]), + xi_dm[abs(boxMapping[2])]: abs(boxMapping_default[2])} + signs = [1 if boxMapping[c]*boxMapping_default[c] > 0 else -1 for c in range(3)] + xi_sign_change = [signs[abs(boxMapping_default[0])-1], signs[abs(boxMapping_default[1])-1], + signs[abs(boxMapping_default[2])-1]] + + lnm = {} + signs = [xi_sign_change[xi_global_to_default_map['xi1'] - 1], + xi_sign_change[xi_global_to_default_map['xi2'] - 1], + xi_sign_change[xi_global_to_default_map['xi3'] - 1]] + for ln in range(1, 9): + xi_l = [] + xi_default = node_xi_map[ln] + xi = [xi_default[xi_global_to_default_map['xi1'] - 1], + xi_default[xi_global_to_default_map['xi2'] - 1], + xi_default[xi_global_to_default_map['xi3'] - 1]] + + for i in range(3): + if signs[i] > 0: + xi_l.append(xi[i]) + else: + xi_l.append(1 - xi[i]) + lnm[ln] = xi_node_map[tuple(xi_l)] + + deriv_default_to_global_map = {xi_dm[abs(boxMapping_default[0])]: abs(boxMapping[0]), + xi_dm[abs(boxMapping_default[1])]: abs(boxMapping[1]), + xi_dm[abs(boxMapping_default[2])]: abs(boxMapping[2])} + signs = [1 if boxMapping[c]*boxMapping_default[c] > 0 else -1 for c in range(3)] + deriv_sign_change = [signs[abs(boxMapping_default[0])-1], signs[abs(boxMapping_default[1])-1], + signs[abs(boxMapping_default[2])-1]] + + self._xi_mapping = {1: xi_default_to_global_map['xi1'], 2: xi_default_to_global_map['xi2'], + 3: xi_default_to_global_map['xi3']} + self._xi_signs = {1: xi_sign_change[0], 2: xi_sign_change[1], 3: xi_sign_change[2]} + self._box_deriv_mapping = {1: deriv_default_to_global_map['xi1'], + 2: deriv_default_to_global_map['xi2'], + 3: deriv_default_to_global_map['xi3']} + self._box_deriv_signs = {1: deriv_sign_change[0], 2: deriv_sign_change[1], 3: deriv_sign_change[2]} + + return lnm + + def get_element_type(self, e3, e2, e1): + """ + Find element type. + :return:octant number, element type and extends for the octant. + """ + if e3 >= self.elementsCountAcross[2] // 2: + if e2 < self.elementsCountAcross[0] // 2: + if e1 >= self.elementsCountAcross[1] // 2: + octant_number = 1 + else: + octant_number = 2 + else: + if e1 < self.elementsCountAcross[1] // 2: + octant_number = 3 + else: + octant_number = 4 + else: + if e2 < self.elementsCountAcross[0] // 2: + if e1 >= self.elementsCountAcross[1] // 2: + octant_number = 5 + else: + octant_number = 6 + else: + if e1 < self.elementsCountAcross[1] // 2: + octant_number = 7 + else: + octant_number = 8 + + if octant_number in [2, 4, 5, 7]: + e3zo = self.elementsCountAcross[2] // 2 - 1 + e2zo = self.elementsCountAcross[1] // 2 - 1 + e1zo = self.elementsCountAcross[0] // 2 - 1 + else: + e3zo = self.elementsCountAcross[2] // 2 - 1 + e2zo = self.elementsCountAcross[0] // 2 - 1 + e1zo = self.elementsCountAcross[1] // 2 - 1 + + e3yo, e2bo, e1yo = e3zo - 1, 1, e1zo - 1 + e3o, e2o, e1o = self.get_local_element_index(octant_number, e3, e2, e1) + + if e3o <= e3yo and e2o >= e2bo and e1o <= e1yo: + element_type = self.ELEMENT_REGULAR + elif e3o == e3yo and e2o == 0 and e1o == e1yo: + element_type = self.ELEMENT_QUADRUPLE_DOWN_LEFT + elif e3o == e3yo and e2o >= e2bo and e1o == e1zo: + element_type = self.ELEMENT_QUADRUPLE_DOWN_RIGHT + elif e3o == e3zo and e2o >= e2bo and e1o == e1yo: + element_type = self.ELEMENT_QUADRUPLE_UP_LEFT + elif e3o == e3yo and e2o == 0 and e1o < e1yo: + element_type = self.ELEMENT_QUADRUPLE_DOWN + elif e3o == e3zo and e2o >= e2bo and e1o < e1yo: + element_type = self.ELEMENT_QUADRUPLE_UP + elif e3o < e3yo and e2o == 0 and e1o == e1yo: + element_type = self.ELEMENT_QUADRUPLE_LEFT + elif e3o < e3yo and e2o == e2bo and e1o == e1zo: + element_type = self.ELEMENT_QUADRUPLE_RIGHT + elif e3o < e3yo and e2o > e2bo and e1o == e1zo: + element_type = self.ELEMENT_DOWN_RIGHT + elif e3o < e3yo and e2o == 0 and e1o < e1yo: + element_type = self.ELEMENT_DOWN_LEFT + else: + element_type = 0 + + return octant_number, element_type, e3zo, e2zo, e1zo + + def get_local_element_index(self, octant_number, e3, e2, e1): + """ + Get octant element index from sphere element index. + :return:e3o, e2o, e1o, local octant element indexes. + """ + if octant_number == 1: + e3o, e2o, e1o = e3 - self.elementsCountAcross[2] // 2, e2, e1 - self.elementsCountAcross[1] // 2 + elif octant_number == 2: + e3o, e2o, e1o = e3 - self.elementsCountAcross[2] // 2, e1, self.elementsCountAcross[0] // 2 - 1 - e2 + elif octant_number == 3: + e3o, e2o, e1o = e3 - self.elementsCountAcross[2] // 2, self.elementsCountAcross[0] - 1 - e2,\ + self.elementsCountAcross[1] // 2 - 1 - e1 + elif octant_number == 4: + e3o, e2o, e1o = e3 - self.elementsCountAcross[2] // 2, self.elementsCountAcross[1] - 1 - e1,\ + e2 - self.elementsCountAcross[0] // 2 + elif octant_number == 5: + e3o, e2o, e1o = self.elementsCountAcross[2]//2 - 1 - e3, self.elementsCountAcross[1] - 1 - e1,\ + self.elementsCountAcross[0]//2 - 1 - e2 + elif octant_number == 6: + e3o, e2o, e1o = self.elementsCountAcross[2]//2 - 1 - e3, e2, self.elementsCountAcross[1] // 2 - 1 - e1 + elif octant_number == 7: + e3o, e2o, e1o = self.elementsCountAcross[2]//2 - 1 - e3, e1, e2 - self.elementsCountAcross[0]//2 + elif octant_number == 8: + e3o, e2o, e1o = self.elementsCountAcross[2]//2 - 1 - e3, self.elementsCountAcross[0] - 1 - e2,\ + e1 - self.elementsCountAcross[1]//2 + + return e3o, e2o, e1o + + def get_global_node_index(self, octant_number, n3o, n2o, n1o): + """ + Get octant element index from sphere element index. + :return: n3, n2, n1, sphere element indexes. + """ + if octant_number == 1: + n3, n2, n1 = n3o + self.elementsCountAcross[2]//2, n2o, n1o + self.elementsCountAcross[1] // 2 + elif octant_number == 2: + n3, n2, n1 = n3o + self.elementsCountAcross[2]//2, self.elementsCountAcross[0] // 2 - n1o, n2o + elif octant_number == 3: + n3, n2, n1 = n3o + self.elementsCountAcross[2]//2, self.elementsCountAcross[0] - n2o,\ + self.elementsCountAcross[1] // 2 - n1o + elif octant_number == 4: + n3, n2, n1 = n3o + self.elementsCountAcross[2]//2, self.elementsCountAcross[0] // 2 + n1o,\ + self.elementsCountAcross[1] - n2o + elif octant_number == 5: + n3, n2, n1 = self.elementsCountAcross[2] // 2 - n3o, self.elementsCountAcross[0]//2 - n1o,\ + self.elementsCountAcross[1] - n2o + elif octant_number == 6: + n3, n2, n1 = self.elementsCountAcross[2] // 2 - n3o, n2o, self.elementsCountAcross[1]//2 - n1o + elif octant_number == 7: + n3, n2, n1 = self.elementsCountAcross[2] // 2 - n3o, n1o + self.elementsCountAcross[0]//2, n2o + elif octant_number == 8: + n3, n2, n1 = self.elementsCountAcross[2] // 2 - n3o, self.elementsCountAcross[0] - n2o,\ + self.elementsCountAcross[1]//2 + n1o + + return n3, n2, n1 + + def get_octant_node_index(self, octant_number, n3, n2, n1): + """ + Get local octant node index from sphere node index. + :return: n3o, n2o, n1o, local octant node indexes. + """ + if octant_number == 1: + n3o, n2o, n1o = n3 - self.elementsCountAcross[2] // 2, n2, n1 - self.elementsCountAcross[1] // 2 + elif octant_number == 2: + n3o, n2o, n1o = n3 - self.elementsCountAcross[2] // 2, n1, self.elementsCountAcross[0] // 2 - n2 + elif octant_number == 3: + n3o, n2o, n1o = n3 - self.elementsCountAcross[2] // 2, self.elementsCountAcross[0] - n2,\ + self.elementsCountAcross[1] // 2 - n1 + elif octant_number == 4: + n3o, n2o, n1o = n3 - self.elementsCountAcross[2] // 2, self.elementsCountAcross[1] - n1,\ + n2 - self.elementsCountAcross[0]//2 + elif octant_number == 5: + n3o, n2o, n1o = self.elementsCountAcross[2] // 2 - n3, self.elementsCountAcross[1] - n1,\ + self.elementsCountAcross[0]//2 - n2 + elif octant_number == 6: + n3o, n2o, n1o = self.elementsCountAcross[2] // 2 - n3, n2, self.elementsCountAcross[1]//2 - n1 + elif octant_number == 7: + n3o, n2o, n1o = self.elementsCountAcross[2] // 2 - n3, n1, n2 - self.elementsCountAcross[0]//2 + elif octant_number == 8: + n3o, n2o, n1o = self.elementsCountAcross[2] // 2 - n3, self.elementsCountAcross[0] - n2,\ + n1 - self.elementsCountAcross[1]//2 + + return n3o, n2o, n1o + + def getNodeId(self, octant_number, n3, n2, n1, n3yo, n1yo): + """ + Get node identifier. Use node indexes n3, n2, n1 to find it. If the indexes represent removed nodes in the + lattice, find the right indexes first. + :return: Node identifier. + """ + n1zo = n1yo + 1 + n3zo = n3yo + 1 + n3o, n2o, n1o = self.get_octant_node_index(octant_number, n3, n2, n1) + if n2o == 0 and n1o == n1yo: + if n3o == n3yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o+1, n2o, n1o+1) + elif n3o == n3yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o-1, n1o + 1) + else: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o, n1o + 1) + elif n2o == 1 and n1o == n1zo: + if n3o == n3yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o+1, n2o-1, n1o) + elif n3o == n3yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o-1, n1o+1) + else: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o-1, n1o) + elif n3o == n3yo and (n2o == 0 or n1o == n1zo): + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o+1, n2o, n1o) + elif n3o == n3zo and n2o == 1 and n1o == n1yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o-1, n1o+1) + elif n3o == n3zo and n2o == 1: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o - 1, n1o) + elif n3o == n3zo and n2o > 1 and n1o == n1yo: + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o, n2o, n1o+1) + else: + n3r, n2r, n1r = n3, n2, n1 + + return self.nodeId[n3r][n2r][n1r] + + def remap_eft_node_value_label(self, eft, localNodeIndexes, node_type, derivatives=None): + """ + Remaps derivatives for common types of nodes. It also checks if each node needs scale factors. If one of the + nodes of an element needs it, then flags it. It assumes the directions of derivatives in a certain way + (based on default [1, 3, 2] xi directions). + If directions are different, then we can use the 'derivatives' argument to fix this. This part is done for + irregular nodes, in the box and for other nodes, derivatives should be given. + :param node_type: e.g. corner 1 represents an end node on the axis1. + :param derivatives: derivatives to be changed to. + :return: + """ + label = {1: Node.VALUE_LABEL_D2_DS2DS3, 2: Node.VALUE_LABEL_D2_DS1DS3, 3: Node.VALUE_LABEL_D2_DS1DS2} + expressionLabel = {1: Node.VALUE_LABEL_D_DS1, 2: Node.VALUE_LABEL_D_DS2, 3: Node.VALUE_LABEL_D_DS3} + sf = {1: [], -1: [1]} + + xim = self._xi_mapping + xis = self._xi_signs + if derivatives: + dm = {1: abs(derivatives[0]), 2: abs(derivatives[1]), 3: abs(derivatives[2])} + signs = [1 if c > 0 else -1 for c in derivatives] + ds = {1: signs[0], 2: signs[1], 3: signs[2]} + else: + dm = self._box_deriv_mapping + ds = self._box_deriv_signs + + # use cross derivatives for swapping derivatives. + remapEftNodeValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS1, [(Node.VALUE_LABEL_D2_DS2DS3, [])]) + remapEftNodeValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D2_DS1DS3, [])]) + remapEftNodeValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D2_DS1DS2, [])]) + + if node_type == self.CORNER_1: + if any(c != 1 for c in [-1 * xis[1] * ds[3], 1 * xis[2] * ds[2], 1 * xis[3] * ds[1]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[3]], + sf[-1 * xis[1] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[1]], + sf[1 * xis[3] * ds[1]])]) + elif node_type == self.CORNER_2: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + elif node_type == self.CORNER_3: + if any(c != 1 for c in [-1 * xis[3] * ds[1], -1 * xis[1] * ds[2], 1 * xis[2] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[2]], + sf[-1 * xis[1] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[3]], + sf[1 * xis[2] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[1]], + sf[-1 * xis[3] * ds[1]])]) + + elif node_type == self.QUADRUPLE_DOWN_LEFT: + if any(c != 1 for c in [-1 * xis[1], 1 * xis[2], 1 * xis[2], -1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[2]]), (Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + elif node_type == self.QUADRUPLE_RIGHT: + if any(c != 1 for c in [1 * xis[1], 1 * xis[3], 1 * xis[2], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[2]]), (Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + elif node_type == self.QUADRUPLE_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + elif node_type == self.QUADRUPLE0_DOWN_LEFT: + if any(c != 1 for c in [1 * xis[3] * ds[3], 1 * xis[2] * ds[2], 1 * xis[1] * ds[1], -1 * xis[1] * ds[2], + -1 * xis[1] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], + [(expressionLabel[dm[1]], sf[1 * xis[1] * ds[1]]), + (expressionLabel[dm[2]], sf[-1 * xis[1] * ds[2]]), + (expressionLabel[dm[3]], sf[-1 * xis[1] * ds[3]])]) + elif node_type == self.QUADRUPLE0_RIGHT: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], -1 * xis[3] * ds[1], 1 * xis[3] * ds[2], + 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], + [(expressionLabel[dm[1]], sf[-1 * xis[3] * ds[1]]), + (expressionLabel[dm[2]], sf[1 * xis[3] * ds[2]]), + (expressionLabel[dm[3]], sf[1 * xis[3] * ds[3]])]) + elif node_type == self.QUADRUPLE0_UP: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[3] * ds[3], -1 * xis[2] * ds[1], 1 * xis[2] * ds[2], + 1 * xis[2] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(expressionLabel[dm[1]], sf[-1 * xis[2] * ds[1]]), + (expressionLabel[dm[2]], sf[1 * xis[2] * ds[2]]), + (expressionLabel[dm[3]], sf[1 * xis[2] * ds[3]])]) + + elif node_type == self.TRIPLE_12_LEFT: + if any(c != 1 for c in [1 * xis[2] * ds[2], -1 * xis[1] * ds[3], 1 * xis[3] * ds[1]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[3]], + sf[-1 * xis[1] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[1]], + sf[1 * xis[3] * ds[1]])]) + elif node_type == self.TRIPLE_12_RIGHT: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + elif node_type == self.TRIPLE_13_DOWN: + if any(c != 1 for c in [-1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + elif node_type == self.TRIPLE_13_UP: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[1]])]) + # remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D2_DS2DS3, [])]) + # temporary to enable swap + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + # remapEftNodeValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D2_DS2DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + # finish swap + elif node_type == self.TRIPLE_23_DOWN: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + elif node_type == self.TRIPLE_23_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + elif node_type == self.TRIPLE0_12_LEFT: + if any(c != 1 for c in [1 * xis[3] * ds[3], 1 * xis[2] * ds[2], 1 * xis[1] * ds[1], -1 * xis[1] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], + [(expressionLabel[dm[1]], sf[1 * xis[1] * ds[1]]), + (expressionLabel[dm[3]], sf[-1 * xis[1] * ds[3]])]) + elif node_type == self.TRIPLE0_12_RIGHT: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], -1 * xis[3] * ds[1], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], + [(expressionLabel[dm[1]], sf[-1 * xis[3] * ds[1]]), (expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + elif node_type == self.TRIPLE0_13_DOWN: + if any(c != 1 for c in [1 * xis[3] * ds[3], 1 * xis[2] * ds[2], 1 * xis[1] * ds[1], -1 * xis[1] * ds[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], + [(expressionLabel[dm[1]], sf[1 * xis[1] * ds[1]]), (expressionLabel[dm[2]], + sf[-1 * xis[1] * ds[2]])]) + elif node_type == self.TRIPLE0_13_Up: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[3] * ds[3], -1 * xis[2] * ds[1], 1 * xis[2] * ds[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(expressionLabel[dm[1]], sf[-1 * xis[2] * ds[1]]), + (expressionLabel[dm[2]], sf[1 * xis[2] * ds[2]])]) + elif node_type == self.TRIPLE0_23_DOWN: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], 1 * xis[3] * ds[2], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], + [(expressionLabel[dm[2]], sf[1 * xis[3] * ds[2]]), + (expressionLabel[dm[3]], sf[1 * xis[3] * ds[3]])]) + elif node_type == self.TRIPLE0_23_UP: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[3] * ds[3], 1 * xis[2] * ds[2], 1 * xis[2] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(expressionLabel[dm[2]], sf[1 * xis[2] * ds[2]]), + (expressionLabel[dm[3]], sf[1 * xis[2] * ds[3]])]) + + elif node_type == self.BOUNDARY_12_LEFT: + if any(c != 1 for c in [1 * xis[2] * ds[2], -1 * xis[1] * ds[3], 1 * xis[3] * ds[1]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[3]], + sf[-1 * xis[1] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[1]], + sf[1 * xis[3] * ds[1]])]) + elif node_type == self.BOUNDARY_12_RIGHT: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + elif node_type == self.BOUNDARY_13_DOWN: + if any(c != 1 for c in [-1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + elif node_type == self.BOUNDARY_13_UP: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[1]])]) + # remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D2_DS2DS3, [])]) + # temporary to enable swap + # remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, [])]) + # remapEftNodeValueLabel(eft, localNodeIndexes, Node.VALUE_LABEL_D2_DS2DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + # finish swaps + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + elif node_type == self.BOUNDARY_23_DOWN: + if any(c != 1 for c in [1 * xis[3], 1 * xis[1], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + elif node_type == self.BOUNDARY_23_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + + elif node_type == self.TRIPLE_CURVE_1_DOWN: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + elif node_type == self.TRIPLE_CURVE_1_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + elif node_type == self.TRIPLE_CURVE_2_DOWN: + if any(c != 1 for c in [-1 * xis[1], 1 * xis[2], -1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + elif node_type == self.TRIPLE_CURVE_2_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + elif node_type == self.TRIPLE_CURVE_3_LEFT: + if any(c != 1 for c in [1 * xis[2], -1 * xis[1], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + elif node_type == self.TRIPLE_CURVE_3_RIGHT: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + elif node_type == self.TRIPLE_CURVE0_1_UP: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[3] * ds[3], 1 * xis[2] * ds[2], 1 * xis[2] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(expressionLabel[dm[2]], sf[1 * xis[2] * ds[2]]), (expressionLabel[dm[3]], + sf[1 * xis[2] * ds[3]])]) + elif node_type == self.TRIPLE_CURVE0_1_DOWN: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], 1 * xis[3] * ds[2], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], + [(expressionLabel[dm[2]], sf[1 * xis[3] * ds[2]]), (expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + elif node_type == self.TRIPLE_CURVE0_2_DOWN: + if any(c != 1 for c in [1 * xis[2] * ds[2], 1 * xis[3] * ds[3], 1 * xis[1] * ds[1], -1 * xis[1] * ds[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], + [(expressionLabel[dm[1]], sf[1 * xis[1] * ds[1]]), (expressionLabel[dm[2]], + sf[-1 * xis[1] * ds[2]])]) + elif node_type == self.TRIPLE_CURVE0_2_UP: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[3] * ds[3], -1 * xis[2] * ds[1], 1 * xis[2] * ds[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], + [(expressionLabel[dm[1]], sf[-1 * xis[2] * ds[1]]), (expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + elif node_type == self.TRIPLE_CURVE0_3_LEFT: + if any(c != 1 for c in [1 * xis[2] * ds[2], 1 * xis[3] * ds[3], 1 * xis[1] * ds[1], -1 * xis[1] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], + [(expressionLabel[dm[1]], sf[1 * xis[1] * ds[1]]), (expressionLabel[dm[3]], + sf[-1 * xis[1] * ds[3]])]) + elif node_type == self.TRIPLE_CURVE0_3_RIGHT: + if any(c != 1 for c in [1 * xis[1] * ds[1], 1 * xis[2] * ds[2], -1 * xis[3] * ds[1], 1 * xis[3] * ds[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(expressionLabel[dm[1]], + sf[1 * xis[1] * ds[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(expressionLabel[dm[2]], + sf[1 * xis[2] * ds[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], + [(expressionLabel[dm[1]], sf[-1 * xis[3] * ds[1]]), (expressionLabel[dm[3]], + sf[1 * xis[3] * ds[3]])]) + + elif node_type == self.SURFACE_REGULAR_DOWN_LEFT: + if any(c != 1 for c in [1 * xis[2], -1 * xis[1], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, + sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS3, + sf[-1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[3]])]) + elif node_type == self.SURFACE_REGULAR_DOWN_RIGHT: + if any(c != 1 for c in [1 * xis[1], 1 * xis[2], 1 * xis[3]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + elif node_type == self.SURFACE_REGULAR_UP: + if any(c != 1 for c in [1 * xis[1], -1 * xis[3], 1 * xis[2]]): + if not self._element_needs_scale_factor: + self._element_needs_scale_factor = True + setEftScaleFactorIds(eft, [1], []) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS2, sf[-1 * xis[3]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[2]])]) + elif node_type == self.REGULAR: + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[1]], [(Node.VALUE_LABEL_D_DS1, sf[1 * xis[1]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[2]], [(Node.VALUE_LABEL_D_DS2, sf[1 * xis[2]])]) + remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[3]], [(Node.VALUE_LABEL_D_DS3, sf[1 * xis[3]])]) + else: + raise ValueError("Remapping for derivatives of this 'node type' is not implemented") + + # for ci in range(1, 4): + # expressionTerms = [] + # for di in range(1, 4): + # if default_sign[ci][di-1]: + # expressionTerms.append((expressionLabel[dm[di]], sf[default_sign[ci][di-1] * xis[ci] * ds[di]])) + # remapEftNodeValueLabel(eft, localNodeIndexes, label[xim[ci]], expressionTerms) diff --git a/src/scaffoldmaker/utils/spheremesh.py b/src/scaffoldmaker/utils/spheremesh.py new file mode 100644 index 00000000..bd0d4125 --- /dev/null +++ b/src/scaffoldmaker/utils/spheremesh.py @@ -0,0 +1,1221 @@ +""" +Utility functions for generating a solid sphere/spheroid (ellipsoid in general). +""" + +from enum import Enum + +import math + +from opencmiss.utils.zinc.finiteelement import getMaximumNodeIdentifier, getMaximumElementIdentifier +from opencmiss.zinc.field import Field +from scaffoldmaker.utils import vector +from scaffoldmaker.utils.interpolation import sampleCubicHermiteCurves, smoothCubicHermiteDerivativesLine +from scaffoldmaker.utils.cylindermesh import Ellipse2D, EllipseShape +from scaffoldmaker.utils.shieldmesh import ShieldMesh3D + + +class SphereShape(Enum): + SPHERE_SHAPE_FULL = 1 + SPHERE_SHAPE_HALF_AAP = 2 # AAP is a hemisphere where x_a3>=0 + SPHERESHIELD_SHAPE_OCTANT_PPP = 3 # positive axis1, positive axis2, positive axis3 + SPHERESHIELD_SHAPE_OCTANT_NPP = 4 + SPHERESHIELD_SHAPE_OCTANT_PNP = 5 + SPHERESHIELD_SHAPE_OCTANT_NNP = 6 + SPHERESHIELD_SHAPE_OCTANT_PPN = 7 + SPHERESHIELD_SHAPE_OCTANT_NPN = 8 + SPHERESHIELD_SHAPE_OCTANT_PNN = 9 + SPHERESHIELD_SHAPE_OCTANT_NNN = 10 + + +class SphereMesh: + """ + Sphere mesh generator. + """ + + def __init__(self, fieldModule, coordinates, centre, axes, elementsCountAcross, + elementsCountAcrossShell, elementsCountAcrossTransition, shellProportion, + sphereShape=SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP, rangeOfRequiredElements=None, + boxDerivatives=None, useCrossDerivatives=False, meshGroups=[]): + """ + :param fieldModule: Zinc fieldModule to create elements in. + :param coordinates: Coordinate field to define. + :param centre, axes: centre and axes of the sphere. + :param elementsCountAcross: [elementsCountAcrossAxis1, elementsCountAcrossAxis2, elementsCountAcrossAxis3] + Total number of elements across the sphere axes. + :param elementsCountAcrossShell, elementsCountAcrossTransition: Total number of elements across each axis + consists of regular elements in the middle cube, transition elements from cube to a sphere (core boundary) + and shell elements around it. Shell nodes and derivatives are similar to the core boundary and don't need + remapping. The topology of the shield structure is extended to 3D with quadruple points. + :param sphereShape: A value from enum SphereShape specifying the shape of sphere. Octant_PPP, for example is + the octant in positive axis1, positive axis2 and positive axis3. SPHERE_SHAPE_HALF_AAP is a hemisphere on + the positive side of axis3. + :param rangeOfRequiredElements: Specifies the range of elements required to be created. It can be used to + create the part of sphere required. If None or same as elementsCountAcross the whole part will be created. + :param boxDerivatives: It is a list of [deriv1,deriv2,deriv3]. It is used to set the derivative directions. + default is [1, 3, 2] which means it makes -d1, d3 and d2 in direction of axis1, axis2 and axis3. To make + d1, d2 and d3 directions in direction of axis1, axis2 and axis3 use [-1, 2, 3]. + """ + + self._axes = axes + self._radius = [vector.magnitude(axis) for axis in axes] + self._coreRadius = [] + self._shield = None + self._elementsCount = elementsCountAcross + self._elementsCountAcrossShell = elementsCountAcrossShell + self._elementsCountAcrossTransition = elementsCountAcrossTransition + self._elementsCountAcrossRim = self._elementsCountAcrossShell + self._elementsCountAcrossTransition - 1 + self._shellProportion = shellProportion + + self._startNodeIdentifier = 1 + self._startElementIdentifier = 1 + self._endNodeIdentifier = 1 + self._endElementIdentifier = 1 + self._sphereShape = sphereShape + self._rangeOfRequiredElements = rangeOfRequiredElements + + self._useCrossDerivatives = useCrossDerivatives + + self._centre = centre + + self._boxDerivatives = boxDerivatives if boxDerivatives else [1, 3, 2] + self._meshGroups = meshGroups + self._shield3D = None + + for i in range(3): + elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - shellProportion) + self._coreRadius.append( + (1 - shellProportion * elementsCountAcrossShell / elementsAxis) * self._radius[i]) + + # generate the mesh + self.createSphereMesh3d(fieldModule, coordinates) + + def createSphereMesh3d(self, fieldModule, coordinates): + """ + Create a sphere mesh based on the shield topology. + :param fieldModule: Zinc fieldModule to create elements in. + :param coordinates: Coordinate field to define. + :return: Final values of nextNodeIdentifier, nextElementIdentifier. + """ + if 'OCTANT' in str(self._sphereShape): + minimum = [2, 2, 2] + era = 0 + elif 'SPHERE_SHAPE_HALF_AAP' in str(self._sphereShape): + minimum = [4, 4, 2] + era = 2 + else: + minimum = [4, 4, 4] + era = 3 + + for i in range(3): + assert (self._elementsCount[i] >= minimum[i]), 'createSphereMesh3d: Invalid number of elements' + + for i in range(era): + assert (self._elementsCount[i] % 2 == 0), 'createSphereMesh3d: number of across elements' \ + ' is not an even number' + + assert (self._sphereShape in [SphereShape.SPHERE_SHAPE_FULL, + SphereShape.SPHERE_SHAPE_HALF_AAP, + SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP]), \ + 'createSphereMesh3d: Invalid sphere mode.' + + elementsCountRim = self._elementsCountAcrossRim + + self._shield3D = ShieldMesh3D(self._elementsCount, elementsCountRim, box_derivatives=self._boxDerivatives) + + nodes = fieldModule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + mesh = fieldModule.findMeshByDimension(3) + + elementsCountAcrossShell = self._elementsCountAcrossShell + elementsCountAcrossTransition = self._elementsCountAcrossTransition + shellProportion = self._shellProportion + + # order of octants is important because common nodes will overwrite + OctantVariationsAll = [SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPN, SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNN, + SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNN, SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPN, + SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP, SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNP, + SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNP, SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPP] + + if self._sphereShape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + OctantVariationsType = [SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP] + elif self._sphereShape == SphereShape.SPHERE_SHAPE_HALF_AAP: + OctantVariationsType = OctantVariationsAll[4:9] + else: + OctantVariationsType = OctantVariationsAll + + for octantType in OctantVariationsType: + axes, elementsCountAcross, boxDerivatives = self.get_octant_axes_and_elements_count(octantType) + octant = OctantMesh(self._centre, axes, elementsCountAcross, + elementsCountAcrossShell, elementsCountAcrossTransition, shellProportion, + sphereShape=SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP, + useCrossDerivatives=False, boxDerivatives=boxDerivatives) + self.copy_octant_nodes_to_sphere_shield(octant, octantType) + + self.modify_octant_common_nodes() + self.sphere_to_spheroid() + + self.generateNodes(nodes, fieldModule, coordinates) + self.generateElements(mesh, fieldModule, coordinates) + + def get_octant_axes_and_elements_count(self, octant_shape): + """ + Get the octants axes, elements count across and box mapping for 8 octants in 8 different regions. + :param octant_shape: A value from enum SphereShape specifying the shape of sphere + :return: axes, elementsCountAcross, boxDerivatives for the octant. + """ + boxDerivativesV1 = self._boxDerivatives + boxDerivativesV2 = [self._boxDerivatives[1], self._boxDerivatives[0], self._boxDerivatives[2]] + axesV1 = [self._axes[0], self._axes[1], self._axes[2]] + axesV2 = [self._axes[1], self._axes[0], self._axes[2]] + elementsCountAcrossV1 = [c for c in self._elementsCount] + elementsCountAcrossV2 = [self._elementsCount[1], self._elementsCount[0], self._elementsCount[2]] + + if octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + octant_number = 1 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNP: + octant_number = 2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNP: + octant_number = 3 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPP: + octant_number = 4 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPN: + octant_number = 5 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNN: + octant_number = 6 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNN: + octant_number = 7 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPN: + octant_number = 8 + else: + raise ValueError("Not implemented.") + + axesSignsV1 = [1 if str(octant_shape).split('_')[-1][c] == 'P' else -1 for c in range(3)] + axesSignsV2 = [axesSignsV1[1], axesSignsV1[0], axesSignsV1[2]] + + signs = self._shield3D.get_octant_signs(octant_number) + if octant_number in [2, 4, 5, 7]: + boxDerivatives = [boxDerivativesV2[c] * signs[c] for c in range(3)] + elementsCountAcross = elementsCountAcrossV2 + axes = [vector.scaleVector(axesV2[c], axesSignsV2[c]) for c in range(3)] + else: + boxDerivatives = [boxDerivativesV1[c] * signs[c] for c in range(3)] + elementsCountAcross = elementsCountAcrossV1 + axes = [vector.scaleVector(axesV1[c], axesSignsV1[c]) for c in range(3)] + + if self._sphereShape == SphereShape.SPHERE_SHAPE_HALF_AAP: + elementsCountAcross[0] = elementsCountAcross[0]//2 + elementsCountAcross[1] = elementsCountAcross[1]//2 + elif self._sphereShape == SphereShape.SPHERE_SHAPE_FULL: + elementsCountAcross[0] = elementsCountAcross[0] // 2 + elementsCountAcross[1] = elementsCountAcross[1] // 2 + elementsCountAcross[2] = elementsCountAcross[2] // 2 + + axes = [vector.normalise(v) for v in axes] + return axes, elementsCountAcross, boxDerivatives + + def copy_octant_nodes_to_sphere_shield(self, octant, octant_shape): + """ + Copy octant nodes to sphere shield. + :param octant: octant instance. + :param octant_shape: A value from enum SphereShape specifying the shape of sphere + :return: + """ + n3a = 0 + n3z = self._elementsCount[2] + n2a = 0 + n2z = self._elementsCount[0] + n1a = 0 + n1z = self._elementsCount[1] + if self._sphereShape == SphereShape.SPHERE_SHAPE_FULL: + if octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + octant_number = 1 + n3a = self._elementsCount[2]//2 + n2z = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNP: + octant_number = 2 + n3a = self._elementsCount[2]//2 + n2z = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNP: + octant_number = 3 + n3a = self._elementsCount[2]//2 + n2a = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPP: + octant_number = 4 + n3a = self._elementsCount[2]//2 + n2a = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPN: + octant_number = 5 + n3z = self._elementsCount[2]//2 + n2z = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNN: + octant_number = 6 + n3z = self._elementsCount[2]//2 + n2z = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNN: + octant_number = 7 + n3z = self._elementsCount[2]//2 + n2a = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPN: + octant_number = 8 + n3z = self._elementsCount[2]//2 + n2a = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + else: + raise ValueError("Not implemented.") + elif self._sphereShape == SphereShape.SPHERE_SHAPE_HALF_AAP: + if octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + octant_number = 1 + n2z = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PNP: + octant_number = 2 + n2z = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NNP: + octant_number = 3 + n2a = self._elementsCount[0]//2 + n1z = self._elementsCount[1]//2 + elif octant_shape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_NPP: + octant_number = 4 + n2a = self._elementsCount[0]//2 + n1a = self._elementsCount[1]//2 + elif self._sphereShape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + octant_number = 1 + else: + raise ValueError("Not implemented.") + + shield3D = octant.get_shield() + for n3 in range(n3a, n3z + 1): + for n2 in range(n2a, n2z + 1): + for n1 in range(n1a, n1z + 1): + n3o, n2o, n1o = self._shield3D.get_octant_node_index(octant_number, n3, n2, n1) + self._shield3D.px[n3][n2][n1] = shield3D.px[n3o][n2o][n1o] + self._shield3D.pd1[n3][n2][n1] = shield3D.pd1[n3o][n2o][n1o] + self._shield3D.pd2[n3][n2][n1] = shield3D.pd2[n3o][n2o][n1o] + self._shield3D.pd3[n3][n2][n1] = shield3D.pd3[n3o][n2o][n1o] + + def modify_octant_common_nodes(self): + """ + + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n1y = n1z - 1 + n2z = self._elementsCount[0] + n3z = self._elementsCount[2] + n3y = n3z - 1 + + # modify pole on highest z. + if self._sphereShape == SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP: + btd1[n3z][n2z][0] = [-c for c in btd1[n3z][n2z][0]] + btd2[n3z][n2z][0] = [-c for c in btd2[n3z][n2z][0]] + elif self._sphereShape == SphereShape.SPHERE_SHAPE_HALF_AAP: + temp = btd1[n3z][n2z//2][n1z//2] + btd1[n3z][n2z//2][n1z//2] = btd2[n3z][n2z//2][n1z//2] + btd2[n3z][n2z//2][n1z//2] = [-c for c in temp] + elif self._sphereShape == SphereShape.SPHERE_SHAPE_FULL: + temp = btd1[n3z][n2z//2][n1z//2] + btd1[n3z][n2z//2][n1z//2] = btd2[n3z][n2z//2][n1z//2] + btd2[n3z][n2z//2][n1z//2] = [-c for c in temp] + + def sphere_to_spheroid(self): + """ + Using the radius in each direction,transform the sphere to ellipsoid. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + for n3 in range(self._elementsCount[2] + 1): + for n2 in range(self._elementsCount[0] + 1): + for n1 in range(self._elementsCount[1] + 1): + if btx[n3][n2][n1]: + btx[n3][n2][n1] = [self._radius[c] * btx[n3][n2][n1][c] for c in range(3)] + btd1[n3][n2][n1] = [self._radius[c] * btd1[n3][n2][n1][c] for c in range(3)] + btd2[n3][n2][n1] = [self._radius[c] * btd2[n3][n2][n1][c] for c in range(3)] + btd3[n3][n2][n1] = [self._radius[c] * btd3[n3][n2][n1][c] for c in range(3)] + + def generateNodes(self, nodes, fieldModule, coordinates): + """ + Create cylinder nodes from coordinates. + :param nodes: nodes from coordinates. + :param fieldModule: Zinc fieldmodule to create nodes in. Uses DOMAIN_TYPE_NODES. + :param coordinates: Coordinate field to define. + """ + nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) + self._startNodeIdentifier = nodeIdentifier + nodeIdentifier = self._shield3D.generateNodes(fieldModule, coordinates, nodeIdentifier, + self._rangeOfRequiredElements) + self._endNodeIdentifier = nodeIdentifier + + def generateElements(self, mesh, fieldModule, coordinates): + """ + Create cylinder elements from nodes. + :param mesh: + :param fieldModule: Zinc fieldmodule to create nodes in. Uses DOMAIN_TYPE_NODES. + :param coordinates: Coordinate field to define. + """ + elementIdentifier = max(1, getMaximumElementIdentifier(mesh) + 1) + self._startElementIdentifier = elementIdentifier + elementIdentifier = self._shield3D.generateElements(fieldModule, coordinates, elementIdentifier, + self._rangeOfRequiredElements, self._meshGroups) + self._endElementIdentifier = elementIdentifier + + +class OctantMesh: + """ + Octant mesh generator. + """ + + def __init__(self, centre, axes, elementsCountAcross, + elementsCountAcrossShell, elementsCountAcrossTransition, shellProportion, + sphereShape=SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP, useCrossDerivatives=False, boxDerivatives=None): + """ + :param centre, axes: centre and axes of the octant. + :param elementsCountAcross: [elementsCountAcrossAxis1, elementsCountAcrossAxis2, elementsCountAcrossAxis3] Total + number of elements across the octant axes. + :param elementsCountAcrossShell, elementsCountAcrossTransition: Total number of elements across each axis + consists of regular elements in the middle cube, transition elements from cube to a sphere (core boundary) + and shell elements around it. Shell nodes and derivatives are similar to the core boundary and don't need + remapping. The topology of the shield structure is extended to 3D with a quadruple points. + :param sphereShape: A value from enum SphereShape specifying one of the 8 octant regions. Octant_PPP for + example, is the octant in positive axis1, positive axis2 and positive axis3. + """ + self._axes = axes + self._radius = [vector.magnitude(axis) for axis in axes] + self._coreRadius = [] + self._shield = None + self._elementsCount = elementsCountAcross + self._elementsCountAcrossShell = elementsCountAcrossShell + self._elementsCountAcrossTransition = elementsCountAcrossTransition + self._elementsCountAcrossRim = self._elementsCountAcrossShell + self._elementsCountAcrossTransition - 1 + self._shellProportion = shellProportion + + self._sphereShape = sphereShape + + self._useCrossDerivatives = useCrossDerivatives + + self._centre = centre + + self._boxDerivatives = boxDerivatives if boxDerivatives else [1, 3, 2] + self._shield3D = None + + for i in range(3): + elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - shellProportion) + self._coreRadius.append( + (1 - shellProportion * elementsCountAcrossShell / elementsAxis) * self._radius[i]) + + # generate the mesh + self.createOctantMesh3d() + + def createOctantMesh3d(self): + """ + Create an octant mesh based on the shield topology. + """ + elementsCountRim = self._elementsCountAcrossRim + + self._shield3D = ShieldMesh3D(self._elementsCount, elementsCountRim) + + self.create_boundary_ellipses_nodes() + self.create_surface_and_interior_nodes() + self._shield3D.set_derivatives(self._boxDerivatives) + + def create_boundary_ellipses_nodes(self): + """ + Use ellipse class to create 3 ellipses needed for the octant. + :return: + """ + centre = self._centre + elementsCountAcrossShell = self._elementsCountAcrossShell + elementsCountAcrossTransition = self._elementsCountAcrossTransition + shellProportion = self._shellProportion + + ellipseAxes = [[self._axes[0], self._axes[1], self._axes[2]], + [[-c for c in self._axes[2]], self._axes[1], self._axes[0]], + [self._axes[0], [-c for c in self._axes[2]], self._axes[1]]] + + coreRadius = [(self._coreRadius[0], self._coreRadius[1]), + (self._coreRadius[2], self._coreRadius[1]), + (self._coreRadius[0], self._coreRadius[2])] + + elementsCount = [(2 * self._elementsCount[0], 2 * self._elementsCount[1]), + (2 * self._elementsCount[2], 2 * self._elementsCount[1]), + (2 * self._elementsCount[0], 2 * self._elementsCount[2])] + + # set derivatives in up and right direction for each ellipse in square region and the circle around. + # We need to do this so the octant derivatives comes out as we want. + squareDerivatives = [[1, 3], [2, 3], [1, -2]] + circleDerivatives = [[1, 3], [2, 3], [-2, 3]] + + for i in range(3): + majorAxis = ellipseAxes[i][0] + minorAxis = ellipseAxes[i][1] + alongAxis = ellipseAxes[i][2] + elementsCountAcrossMajor = elementsCount[i][0] + elementsCountAcrossMinor = elementsCount[i][1] + coreMajorRadius = coreRadius[i][0] + coreMinorRadius = coreRadius[i][1] + ellipse = Ellipse2D(centre, majorAxis, minorAxis, + elementsCountAcrossMajor, elementsCountAcrossMinor, elementsCountAcrossShell, + elementsCountAcrossTransition, shellProportion, coreMajorRadius, coreMinorRadius, + ellipseShape=EllipseShape.Ellipse_SHAPE_FULL) + + self.copy_ellipses_nodes_to_shield_nodes(ellipse, i+1, squareDerivatives[i], circleDerivatives[i]) + + def copy_ellipses_nodes_to_shield_nodes(self, ellipse, ellipse_number, squareDerivatives, circleDerivatives): + """ + Copy coordinates and derivatives of ellipse to shield. + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + shield = ellipse.getShield() + + # Modify the shield to get only the quarter that you want. make others None so you can generate the nodes. + shield.remap_derivatives(squareDerivatives, circleMapping=circleDerivatives) + for n3 in range(self._elementsCount[2] + 1): + for n2 in range(self._elementsCount[0] + 1): + for n1 in range(self._elementsCount[1] + 1): + n3s, n2s, n1s = n3, n2, n1 + n3e, n2e, n1e = 0, n2, n1 + if ellipse_number == 1: + if n3 > 0: + continue + if n2 == 0 and n1 == self._elementsCount[1] - 1: + n1s = n1 + 1 + n1e = n1 + self._elementsCount[1] + elif ellipse_number == 2: + if n2 < self._elementsCount[0]: + continue + if n3 == self._elementsCount[2] and n1 == self._elementsCount[1] - 1: + n1s = n1 + 1 + n2e = n3 + self._elementsCount[2] + n1e = n1 + self._elementsCount[1] + if n3 == 0: + if n1 == self._elementsCount[1]: + btd2[n3][n2s][n1s] = shield.pd2[0][n2e][n1e] + else: + btd2[n3][n2s][n1s] = shield.pd2[0][n2e][n1e] + continue + elif ellipse_number == 3: + if n1 > 0: + continue + if n2 == 0 and n3 == self._elementsCount[2] - 1: + n3s = self._elementsCount[2] + n1e = self._elementsCount[2] - n3 + if n3 == 0: + if n2 == 0: + btd2[n3][n2][n1s] = shield.pd2[0][n2][n1e] + elif 0 < n2 < self._elementsCount[0]: + btd2[n3][n2][n1s] = shield.pd2[0][n2][n1e] + continue + else: + if n2 == self._elementsCount[0]: + if n3 == self._elementsCount[2]: + btd1[n3][n2][n1s] = shield.pd2[0][n2][n1e] + else: + btd1[n3][n2][n1s] = shield.pd1[0][n2][n1e] + continue + + if shield.px[n3e][n2e][n1e]: + btx[n3s][n2s][n1s] = shield.px[n3e][n2e][n1e] + btd1[n3s][n2s][n1s] = shield.pd1[n3e][n2e][n1e] + btd2[n3s][n2s][n1s] = shield.pd2[n3e][n2e][n1e] + btd3[n3s][n2s][n1s] = shield.pd3[n3e][n2e][n1e] + + def create_surface_and_interior_nodes(self): + """ + Create octant exterior surface nodes and interior nodes + :return: + """ + self.calculate_surface_quadruple_point() + self.sample_triple_curves_on_sphere() + self.sample_regular_curves_on_sphere() + self.create_interior_nodes() + self.smooth_derivatives_to_surface() + + def calculate_surface_quadruple_point(self): + """ + Calculate coordinates and derivatives of the quadruple point on the surface, where 3 hex elements merge. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n1y = n1z - 1 + n2z = self._elementsCount[0] + n3z = self._elementsCount[2] + n3y = n3z - 1 + + elementsAroundEllipse12 = self._elementsCount[0] + self._elementsCount[1] - 2 + radiansAroundEllipse12 = math.pi / 2 + radiansPerElementAroundEllipse12 = radiansAroundEllipse12 / elementsAroundEllipse12 + elementsAroundEllipse13 = self._elementsCount[0] + self._elementsCount[2] - 2 + radiansAroundEllipse13 = math.pi / 2 + radiansPerElementAroundEllipse13 = radiansAroundEllipse13 / elementsAroundEllipse13 + + theta_2 = n3y * radiansPerElementAroundEllipse13 + theta_3 = n1y * radiansPerElementAroundEllipse12 + phi_3 = calculate_azimuth(theta_3, theta_2) + # ratio = -0.1 * (min(self._elementsCount) - 2) + 1 if self._elementsCount[0] <= 2 else 0.2 + ratio = 1 + # local_x = intersection_of_two_great_circles_on_sphere(btx[0][0][n1y-1], + # btx[n3z][n2z][n1z], btx[0][2][n1z], btx[n3z][0][0]) + local_x = spherical_to_cartesian(1.0, theta_3, ratio * phi_3 + (1-ratio)*math.pi/2) + + x = local_to_global_coordinates(local_x, self._axes, self._centre) + + a1, a2, a3 = local_orthogonal_unit_vectors(x, self._axes[2], self._centre) + n3r, n2r, n1r = self.get_triple_curves_end_node_parameters(1, index_output=True) + btx[n3r][n2r][n1r] = x + btd1[n3r][n2r][n1r] = a1 # initialise + btd2[n3r][n2r][n1r] = a2 # initialise + btd3[n3r][n2r][n1r] = a3 + + def sample_triple_curves_on_sphere(self): + """ + Sample points on the triple curves of the 'quadruple point' on the sphere surface. + :return: + """ + # sample on curve 1 of the triple curves and smooth the end derivatives. + n3r1, n2r1, n1r1 = self.get_triple_curves_end_node_parameters(1, index_output=True) + n3r2, n2r2, n1r2 = self.get_triple_curves_end_node_parameters(1, cx=1, index_output=True) + self.sample_curves_between_two_nodes_on_sphere([n3r1, n2r1, n1r1], [n3r2, n2r2, n1r2], + self._elementsCount[0] - 1, [1, None], [1], [1]) + # curve 2 + n3r1, n2r1, n1r1 = self.get_triple_curves_end_node_parameters(1, cx=2, index_output=True) + n3r2, n2r2, n1r2 = self.get_triple_curves_end_node_parameters(1, index_output=True) + self.sample_curves_between_two_nodes_on_sphere([n3r1, n2r1, n1r1], [n3r2, n2r2, n1r2], + self._elementsCount[1] - 1, [1], [-2], [None, -2]) + # curve 3. + n3r1, n2r1, n1r1 = self.get_triple_curves_end_node_parameters(1, cx=3, index_output=True) + n3r2, n2r2, n1r2 = self.get_triple_curves_end_node_parameters(1, index_output=True) + self.sample_curves_between_two_nodes_on_sphere([n3r1, n2r1, n1r1], [n3r2, n2r2, n1r2], + self._elementsCount[2] - 1, [2], [2], [None, None]) + + def sample_regular_curves_on_sphere(self): + """ + Create all other nodes on the sphere except the nodes on the triple curves and ellipses. + :return: + """ + n2a = 1 + n1z = self._elementsCount[1] + n2z = self._elementsCount[0] + n3z = self._elementsCount[2] + + # regular curves crossing curve 1 + for n2 in range(2, self._elementsCount[0]): + # bottom right + self.sample_curves_between_two_nodes_on_sphere([0, n2, n1z], [n3z, n2, n1z], self._elementsCount[2] - 1, + [2], [2], [None, None]) + # top + self.sample_curves_between_two_nodes_on_sphere([n3z, n2, 0], [n3z, n2, n1z], self._elementsCount[1] - 1, + [1], [-2], [None, -2]) + + # regular curves crossing curve 2 + for n1 in range(1, self._elementsCount[1] - 1): + # bottom left. Top is done before. + self.sample_curves_between_two_nodes_on_sphere([0, 0, n1], [n3z, 0, n1], self._elementsCount[2] - 1, + [2], [2], [None, None]) + + # smooth regular curves crossing curve 1 + for n2 in range(n2a + 1, n2z): + self.smooth_derivatives_regular_surface_curve(2, n2, [[1], [None, None]], [[-2], [-2]], [[None, -2], [-2]]) + + # smooth regular curves crossing curve 2 + for n1 in range(1, n1z - 1): + self.smooth_derivatives_regular_surface_curve(1, n1, [[2], [None, None]], [[2], [1]], [[None, 1], [1]]) + + # smooth regular curves crossing curve 3 + for n3 in range(1, self._elementsCount[2] - 1): + self.smooth_derivatives_regular_surface_curve(3, n3, [[1], [None, None]], [[1], [1]], [[None, 1], [1]]) + + def create_interior_nodes(self): + """ + Create box nodes. + :return: + """ + self.calculate_interior_quadruple_point() + self.sample_interior_curves() + self.smooth_regular_interior_curves() + + def smooth_derivatives_to_surface(self): + """ + Smooth derivatives leading to quadruple point where 3 hex elements merge. + :return: + """ + n3z = self._elementsCount[2] + n1z = self._elementsCount[1] + + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + for n3 in range(1, self._elementsCount[2] + 1): + for n2 in range(self._elementsCount[0]): + for n1 in range(1, self._elementsCount[1]+1): + if self.on_sphere(n3, n2, n1): + # find indices of the neighbour node inside and contribution of its derivatives. + n3r = n3 - 1 if n3 == n3z else n3 + n2r = n2 + 1 if n2 == 0 else n2 + n1r = n1 - 1 if n1 == n1z else n1 + co = [-1, 1, 1] + co[0] = -1 if n2 == 0 else 0 + co[1] = 1 if n3 == n3z else 0 + co[2] = 1 if n1 == n1z else 0 + + tx = [] + td3 = [] + tx.append(btx[n3r][n2r][n1r]) + td3.append( + [(co[0]*btd1[n3r][n2r][n1r][c] + co[1]*btd2[n3r][n2r][n1r][c] + + co[2]*btd3[n3r][n2r][n1r][c]) for c in range(3)]) + + tx.append(btx[n3][n2][n1]) + td3.append(btd3[n3][n2][n1]) + + td3 = smoothCubicHermiteDerivativesLine(tx, td3, fixStartDirection=True, fixEndDirection=True) + btd3[n3][n2][n1] = td3[1] + + def sample_curves_between_two_nodes_on_sphere(self, id1, id2, elementsOut, dStart, dbetween, dEnd): + """ + Samples curves on the sphere surface between two points given by their indexes. + :param id1, id2: [n3,n2,n1] for the first and second points. + :param dStart, dBetween, dEnd: Specifies the derivatives that are used for this curve at the beginning, end and + in between. e.g. dStart=[2, -1, None] means d2 for the first node, -d1 for the second node and skip the third + one. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + # Find what index is constant + if id1[0] != id2[0]: + varying_index = 3 + elementsCount = self._elementsCount[2] + elif id1[1] != id2[1]: + varying_index = 2 + elementsCount = self._elementsCount[0] + elif id1[2] != id2[2]: + varying_index = 1 + elementsCount = self._elementsCount[1] + + else: + raise ValueError("None of n1, n2, or n3 is constant. Only on the constant curves.") + + btd = {1: btd1, 2: btd2, 3: btd3} + idi = {0: id1[0], 1: id1[1], 2: id1[2]} + + nx, nd1 = sample_curves_on_sphere(btx[id1[0]][id1[1]][id1[2]], btx[id2[0]][id2[1]][id2[2]], self._centre, + elementsOut) + + nit = 0 + for ni in range(elementsCount + 1): + idi[3 - varying_index] = ni + + if ni < len(dStart): + if dStart[ni]: + btd[dStart[ni]][idi[0]][idi[1]][idi[2]] = nd1[nit] if dStart[ni] > 0 else vector.scaleVector( + nd1[nit], -1) + nit += 1 + elif ni > elementsCount - len(dEnd): + nie = ni - elementsCount + len(dEnd) - 1 + if dEnd[nie]: + btd[abs(dEnd[nie])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dEnd[nie] > 0 else vector.scaleVector( + nd1[nit], -1) + nit += 1 + else: + btx[idi[0]][idi[1]][idi[2]] = nx[nit] + + a1, a2, a3 = local_orthogonal_unit_vectors(nx[nit], self._axes[2], self._centre) + btd1[idi[0]][idi[1]][idi[2]] = a1 # initialise + btd2[idi[0]][idi[1]][idi[2]] = a2 # initialise + btd3[idi[0]][idi[1]][idi[2]] = a3 # initialise + + btd[abs(dbetween[0])][idi[0]][idi[1]][idi[2]] = nd1[nit] if dbetween[0] > 0 else vector.scaleVector( + nd1[nit], -1) + nit += 1 + + def smooth_derivatives_regular_surface_curve(self, constant_index, nc, dStart, dBetween, dEnd): + """ + Smooth derivatives for each constant index curve. e.g. n2 = 3 + :param constant_index: Specifies n1, n2 or n3 is constant. + :param nc: Index that is constant across the curve. + :param dStart, dBetween, dEnd: See sample_curves_between_two_nodes_on_sphere. The difference here is + the values are given for two curves that connect one end to the other end of the sphere surface. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n2z = self._elementsCount[0] + n3z = self._elementsCount[2] + + if constant_index == 1: + elementsCount = [self._elementsCount[2], self._elementsCount[0]] + elif constant_index == 2: + elementsCount = [self._elementsCount[1], self._elementsCount[2]] + elif constant_index == 3: + elementsCount = [self._elementsCount[1], self._elementsCount[0]] + + btd = {1: btd1, 2: btd2, 3: btd3} + + tx = [] + td = [] + for se in range(2): + for ni in range(elementsCount[se] + 1): + if constant_index == 1: + ids = [ni, 0, nc] if se == 0 else [n3z, ni, nc] + elif constant_index == 2: + ids = [n3z, nc, ni] if se == 0 else [elementsCount[se] - ni, nc, n1z] + elif constant_index == 3: + ids = [nc, 0, ni] if se == 0 else [nc, ni, n1z] + + if ni < len(dStart[se]): + if dStart[se][ni]: + tx.append(btx[ids[0]][ids[1]][ids[2]]) + if dStart[0][ni] > 0: + td.append(btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]]) + else: + td.append(vector.scaleVector(btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]], -1)) + elif ni > elementsCount[se] - len(dEnd[se]): + nie = ni - elementsCount[se] + len(dEnd[se]) - 1 + if dEnd[se][nie]: + tx.append(btx[ids[0]][ids[1]][ids[2]]) + if dEnd[se][nie] > 0: + td.append(btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]]) + else: + td.append(vector.scaleVector(btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]], -1)) + else: + tx.append(btx[ids[0]][ids[1]][ids[2]]) + if dBetween[se][0] > 0: + td.append(btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]]) + else: + td.append(vector.scaleVector(btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]], -1)) + + td = smoothCubicHermiteDerivativesLine(tx, td, fixStartDirection=True, fixEndDirection=True) + + nit = 0 + for se in range(2): + for ni in range(elementsCount[se] + 1): + if constant_index == 1: + ids = [ni, 0, nc] if se == 0 else [n3z, ni, nc] + elif constant_index == 2: + ids = [n3z, nc, ni] if se == 0 else [elementsCount[se] - ni, nc, n1z] + elif constant_index == 3: + ids = [nc, 0, ni] if se == 0 else [nc, ni, n1z] + + if ni < len(dStart[se]): + if dStart[se][ni]: + if dStart[0][ni] > 0: + btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]] = td[nit] + else: + btd[abs(dStart[se][ni])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + nit += 1 + elif ni > elementsCount[se] - len(dEnd[se]): + nie = ni - elementsCount[se] + len(dEnd[se]) - 1 + if dEnd[se][nie]: + if dEnd[se][nie] > 0: + btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]] = td[nit] + else: + btd[abs(dEnd[se][nie])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + nit += 1 + else: + if dBetween[se][0] > 0: + btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]] = td[nit] + else: + btd[abs(dBetween[se][0])][ids[0]][ids[1]][ids[2]] = vector.scaleVector(td[nit], -1) + nit += 1 + + def calculate_interior_quadruple_point(self): + """ + Calculate coordinates and derivatives of the quadruple point inside, where 3 hex elements merge. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n1y = n1z - 1 + n3z = self._elementsCount[2] + n3y = n3z - 1 + n2z = self._elementsCount[0] + + if self._elementsCount[2] == min(self._elementsCount): + cx = 3 + elif self._elementsCount[1] == min(self._elementsCount): + cx = 2 + else: + cx = 1 + n3r0, n2r0, n1r0 = self.get_triple_curves_end_node_parameters(0, cx=cx, index_output=True) + n3r, n2r, n1r = self.get_triple_curves_end_node_parameters(1, cx=cx, index_output=True) + + ts = vector.magnitude(vector.addVectors([btx[n3r0][n2r0][n1r0], btx[n3r][n2r][n1r]], [1, -1])) + ra = vector.addVectors([btx[n3z][0][n1z], self._centre], [1, -1]) + radius = vector.magnitude(ra) + local_x = vector.scaleVector(ra, (1 - ts/radius)) + x = vector.addVectors([local_x, self._centre], [1, 1]) + n3r0, n2r0, n1r0 = self.get_triple_curves_end_node_parameters(0, index_output=True) + n3r1, n2r1, n1r1 = self.get_triple_curves_end_node_parameters(1, index_output=True) + btx[n3r0][n2r0][n1r0] = x + btd1[n3r0][n2r0][n1r0] = [-(btx[n3r1][n2r1][n1r1][0] - btx[n3r0][n2r0][n1r0][0]), 0.0, 0.0] + btd2[n3r0][n2r0][n1r0] = [0.0, 0.0, (btx[n3r1][n2r1][n1r1][2] - btx[n3r0][n2r0][n1r0][2])] + btd3[n3r0][n2r0][n1r0] = [0.0, (btx[n3r1][n2r1][n1r1][1] - btx[n3r0][n2r0][n1r0][1]), 0.0] + + def sample_interior_curves(self): + """ + Sample box curves. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n1y = n1z - 1 + n3z = self._elementsCount[2] + n3y = n3z - 1 + n2z = self._elementsCount[0] + + tx, td1 = sampleCubicHermiteCurves([btx[n3y][1][n1y], btx[n3y][n2z][n1y]], + [btd1[0][1][n1y], btd1[0][n2z][n1y]], self._elementsCount[0] - 1)[:2] + + for n2 in range(2, self._elementsCount[0]): + btx[n3y][n2][n1y] = tx[n2-1] + btd1[n3y][n2][n1y] = td1[n2-1] + btd2[n3y][n2][n1y] = [0.0, 0.0, (btx[n3z][n2][n1z][2] - btx[n3y][n2][n1y][2])] + btd3[n3y][n2][n1y] = [0.0, (btx[n3z][n2][n1z][1] - btx[n3y][n2][n1y][1]), 0.0] + + # curve 2 and parallel curves + for n2 in range(1, self._elementsCount[0]): + tx, td3 = sampleCubicHermiteCurves([btx[n3y][n2][0], btx[n3y][n2][n1y]], + [btd3[0][n2][0], btd3[0][n2][n1y]], self._elementsCount[1] - 1)[:2] + + for n1 in range(1, self._elementsCount[1] - 1): + btx[n3y][n2][n1] = tx[n1] + btd3[n3y][n2][n1] = td3[n1] + + for n2 in range(1, self._elementsCount[0]): + for n1 in range(1, self._elementsCount[1] - 1): + if n2 == 1: + btd1[n3y][n2][n1] = [btx[n3y][1][n1][0] - btx[n3z][0][n1][0], 0.0, 0.0] + btd2[n3y][n2][n1] = [0.0, 0.0, -btx[n3y][1][n1][2] + btx[n3z][0][n1][2]] + else: + btd1[n3y][n2][n1] = vector.addVectors([btx[n3y][n2][n1], btx[n3y][n2+1][n1]], [-1, 1]) + btd2[n3y][n2][n1] = vector.addVectors([btx[n3y][n2][n1], btx[0][n2][n1]], [1, -1]) + + # sample along curve0_3 + for n2 in range(1, self._elementsCount[0]): + for n1 in range(1, self._elementsCount[1]): + tx, td2 = sampleCubicHermiteCurves([btx[0][n2][n1], btx[n3y][n2][n1]], + [btd2[0][n2][0], btd2[n3y][n2][0]], self._elementsCount[2]-1)[:2] + + for n3 in range(1, self._elementsCount[2] - 1): + btx[n3][n2][n1] = tx[n3] + btd2[n3][n2][n1] = td2[n3] + + for n3 in range(1, self._elementsCount[2] - 1): + for n2 in range(1, self._elementsCount[0]): + for n1 in range(1, self._elementsCount[1]): + if n2 == 1 and n1 == n1y: + btd1[n3][n2][n1] = [btx[n3][n2][n1][0] - btx[n3][n2-1][n1+1][0], 0.0, 0.0] + btd3[n3][n2][n1] = [0.0, btx[n3][n2-1][n1+1][1] - btx[n3][n2][n1][1], 0.0] + else: + btd1[n3][n2][n1] = vector.addVectors([btx[n3][n2+1][n1], btx[n3][n2][n1]], [1, -1]) + btd3[n3][n2][n1] = vector.addVectors([btx[n3][n2][n1+1], btx[n3][n2][n1]], [1, -1]) + + def smooth_regular_interior_curves(self): + """ + Smooth box curves. + :return: + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n1z = self._elementsCount[1] + n1y = n1z - 1 + n3z = self._elementsCount[2] + n3y = n3z - 1 + n2z = self._elementsCount[0] + + # smooth d1 in regular 1 + if self._elementsCount[0] >= 3: + for n3 in range(1, self._elementsCount[2]): + for n1 in range(1, self._elementsCount[1]): + tx = [] + td1 = [] + for n2 in range(1, self._elementsCount[0]+1): + tx.append(btx[n3][n2][n1]) + td1.append(btd1[n3][n2][n1]) + td1 = smoothCubicHermiteDerivativesLine(tx, td1, fixEndDirection=True) + for n2 in range(1, self._elementsCount[0]+1): + btd1[n3][n2][n1] = td1[n2-1] + else: + for n3 in range(1, self._elementsCount[2]): + for n1 in range(1, self._elementsCount[1]): + btd1[n3][1][n1] = vector.addVectors([btx[n3][2][n1], btx[n3][1][n1]], [1, -1]) + btd1[n3][2][n1] = vector.setMagnitude(btd1[n3][2][n1], vector.magnitude(btd1[n3][1][n1])) + + # smooth d3 in regular + if self._elementsCount[1] >= 3: + for n3 in range(1, self._elementsCount[2]): + for n2 in range(1, self._elementsCount[0]): + tx = [] + td3 = [] + for n1 in range(self._elementsCount[1]): + tx.append(btx[n3][n2][n1]) + td3.append(btd3[n3][n2][n1]) + + td3 = smoothCubicHermiteDerivativesLine(tx, td3, fixStartDirection=True) + + for n1 in range(self._elementsCount[1]): + btd3[n3][n2][n1] = td3[n1] + else: + for n3 in range(1, self._elementsCount[2]): + for n2 in range(1, self._elementsCount[0]): + btd3[n3][n2][1] = vector.addVectors([btx[n3][n2][1], btx[n3][n2][0]], [1, -1]) + btd3[n3][n2][0] = vector.setMagnitude(btd3[n3][n2][0], vector.magnitude(btd3[n3][n2][1])) + + # regular curves d2 + for n2 in range(1, self._elementsCount[0]): + for n1 in range(1, self._elementsCount[1]): + if self._elementsCount[2] >= 3: + tx = [] + td2 = [] + for n3 in range(self._elementsCount[2]): + tx.append(btx[n3][n2][n1]) + td2.append(btd2[n3][n2][n1]) + td2 = smoothCubicHermiteDerivativesLine(tx, td2, fixStartDirection=True) + for n3 in range(self._elementsCount[2]): + btd2[n3][n2][n1] = td2[n3] + else: + btd2[1][n2][n1] = vector.addVectors([btx[1][n2][n1], btx[0][n2][n1]], [1, -1]) + btd2[0][n2][n1] = vector.setMagnitude(btd2[0][n2][n1], vector.magnitude(btd2[1][n2][n1])) + + def get_triple_curves_end_node_parameters(self, rx, cx=None, index_output=False): + """ + Find the indexes or node parameters for the 6 end nodes of unique curves of triple curves on the surface and + inside. + if cx is not given, it returns the quadruple points identified by rx. + :param cx: curve index. Curve 1 connects quadruples to ellipse 23. + Similarly, curve 2, connects quadruples to ellipse 13 + :param rx: 1 means on the exterior surface, 0 in one below. + :return: indexes of the point or the node parameters. + """ + btx = self._shield3D.px + btd1 = self._shield3D.pd1 + btd2 = self._shield3D.pd2 + btd3 = self._shield3D.pd3 + + n3z = self._elementsCount[2] + n2z = self._elementsCount[0] + n1z = self._elementsCount[1] + n3y = n3z - 1 + n1y = n1z - 1 + + if not cx: + n3r = n3y + rx + n2r = 1 - rx + n1r = n1y + rx + else: + if cx == 1: + n3r = n3y + rx + n2r = n2z + n1r = n1y + rx + elif cx == 2: + n3r = n3y + rx + n2r = 1 - rx + n1r = 0 + elif cx == 3: + n3r = 0 + n2r = 1 - rx + n1r = n1y + rx + else: + raise ValueError("curve index must be 1,2 or 3.") + + if index_output: + return n3r, n2r, n1r + else: + return btx[n3r][n2r][n1r], btd1[n3r][n2r][n1r], btd2[n3r][n2r][n1r], btd3[n3r][n2r][n1r] + + def on_sphere(self, n3, n2, n1): + """ + Check if the given point is on the sphere. + :param n3, n2, n1: node indexes in data structure [n3][n2][n1] + :return: True, if it is on the sphere. + """ + n3z = self._elementsCount[2] + n1z = self._elementsCount[1] + + btx = self._shield3D.px + + return (n3 == n3z or n2 == 0 or n1 == n1z) and btx[n3][n2][n1] + + def get_shield(self): + return self._shield3D + + +def calculate_azimuth(theta, theta_p): + """ + Given polar angles of a point on the sphere surfaces, calculate the azimuth angle. + :param theta: polar angle. In orthonormal coordinate system (axis1, axis2, axis3) with right-hand rule, + theta is angle between common axis and point projection on plane of theta. In case theta=theta_3 and + theta_p = theta_1, theta is between axis2 and projection + :param theta_p: polar angle wrt other direction. + :return: Azimuth angle. + """ + return math.atan(1/(math.tan(theta_p)*math.cos(theta))) + + +def local_orthogonal_unit_vectors(x, axis3, origin): + """ + Find local orthogonal unit vectors for a point on a sphere + :param x: coordinates of the point. + :param axis3: The third axis in Cartesian coordinate system (axis1, axis2, axis3) + :return: e1, e2, e3. Unit vectors. e3 is normal to the boundary, e2 is in (e3, axis3) plane and e1 normal to them. + """ + r = vector.addVectors([x, origin], [1, -1]) + e3 = vector.normalise(r) + e2 = vector.vectorRejection(axis3, e3) + e2 = vector.normalise(e2) + e1 = vector.crossproduct3(e2, e3) + + return e1, e2, e3 + + +def calculate_arc_length(x1, x2, origin): + """ + Calculate the arc length between points x1 and x2. + :param x1, x2: points coordinates. + :return: arc length + """ + r1 = vector.addVectors([x1, origin], [1, -1]) + r2 = vector.addVectors([x2, origin], [1, -1]) + radius = vector.magnitude(r1) + angle = vector.angleBetweenVectors(r1, r2) + return radius * angle + + +def sample_curves_on_sphere(x1, x2, origin, elementsOut): + """ + + :param x1, x2: points coordinates. + :param elementsOut: + :return: + """ + r1 = vector.addVectors([x1, origin], [1, -1]) + r2 = vector.addVectors([x2, origin], [1, -1]) + deltax = vector.addVectors([r1, r2], [-1, 1]) + normal = vector.crossproduct3(r1, deltax) + angle = vector.angleBetweenVectors(r1, r2) + anglePerElement = angle/elementsOut + arcLengthPerElement = calculate_arc_length(x1, x2, origin)/elementsOut + + nx = [] + nd1 = [] + for n1 in range(elementsOut + 1): + radiansAcross = n1 * anglePerElement + r = vector.rotateVectorAroundVector(r1, normal, radiansAcross) + x = vector.addVectors([r, origin], [1, 1]) + d1 = vector.setMagnitude(vector.crossproduct3(normal, r), arcLengthPerElement) + nx.append(x) + nd1.append(d1) + + return nx, nd1 + + +def spherical_to_cartesian(r, theta, phi): + """ + :param r: Radius. + :param theta: in radians. + :param phi: azimuth angle in radians + :return: x=[x1, x2, x3] coordinates. + """ + return [r*math.sin(phi)*math.cos(theta), r*math.sin(phi)*math.sin(theta), r*math.cos(phi)] + + +def cartesian_to_spherical(x): + """ + :return: [r, theta, phi]. + """ + r = vector.magnitude(x) + theta = math.atan2(x[1], x[0]) + phi = math.acos(x[2]/r) + return r, theta, phi + + +def local_to_global_coordinates(local_x, local_axes, local_origin=None): + """ + Get global coordinates of a point with local coordinates x = [x1, x2, x3] and axes of local coordinate system. + :param local_x: Coordinates in local coordinates system as a list of 3 components. + :param local_origin: Origin of local coordinates system specified as a list of 3 components wrt global coordinates + system. + :param local_axes: Axes of local coordinates system, specified as a list of list 3X3 with respect to global + coordinates system. + :return: Global coordinate system. + """ + if local_origin is None: + local_origin = [0.0, 0.0, 0.0] + normalised_axes = [vector.normalise(v) for v in local_axes] + return vector.addVectors([vector.addVectors(normalised_axes, local_x), local_origin]) + + +def intersection_of_two_great_circles_on_sphere(p1, q1, p2, q2): + """ + Find the intersection between arcs P1Q1 and P2Q2 on sphere. + :param p1, q1, p2, q2: arcs extremities coordinates. + :return: Point Sx, intersection between the arcs. + """ + normal_to_plane_OP1Q1 = vector.crossproduct3(p1, q1) + normal_to_plane_OP2Q2 = vector.crossproduct3(p2, q2) + + planes_intersection_vector = vector.crossproduct3(normal_to_plane_OP1Q1, normal_to_plane_OP2Q2) + if vector.magnitude(planes_intersection_vector) == 0: + sx = None + else: + sx = vector.setMagnitude(planes_intersection_vector, vector.magnitude(p1)) + p1q1_angle = vector.angleBetweenVectors(p1, q1) + p1s_angle = vector.angleBetweenVectors(p1, sx) + p2s_angle = vector.angleBetweenVectors(p2, sx) + if p1s_angle > p1q1_angle or p2s_angle > p1q1_angle: + sx = vector.scaleVector(sx, -1) + + return sx + + +def point_projection_on_sphere(p1, radius): + """ + Find closest point to p1 on the sphere. + :param p1: point. + :return: + """ + return vector.setMagnitude(p1, radius) diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py index 340753a1..510e3b58 100644 --- a/src/scaffoldmaker/utils/vector.py +++ b/src/scaffoldmaker/utils/vector.py @@ -36,12 +36,20 @@ def setMagnitude(v, mag): scale = mag/math.sqrt(sum(c*c for c in v)) return [ c*scale for c in v ] -def addVectors(v1,v2,s1=1.0,s2=1.0): +def addVectors(vectors, scalars = None): ''' - returns s1*v1+s2*v2 where s1 and s2 are scalars. - :return: Vector s1*v1+s2*v2 + returns s1*v1+s2*v2+... where scalars = [s1, s2, ...] and vectors=[v1, v2, ...]. + :return: Resultant vector ''' - return [(s1 * v1[c] + s2 * v2[c]) for c in range(len(v1))] + if not scalars: + scalars = [1]*len(vectors) + else: + assert len(vectors) == len(scalars) + + resultant = [0, 0, 0] + for i in range(len(vectors)): + resultant = [resultant[c] + scalars[i] * vectors[i][c] for c in range(3)] + return resultant def scalarProjection(v1, v2): @@ -57,7 +65,7 @@ def vectorProjection(v1, v2): :return: A projection vector. """ s1 = scalarProjection(v1, v2) - return scalarProduct(s1, normalise(v2)) + return scaleVector(normalise(v2), s1) def vectorRejection(v1, v2): @@ -66,17 +74,17 @@ def vectorRejection(v1, v2): :return: A rejection vector. """ v1p = vectorProjection(v1, v2) - return addVectors(v1, v1p, 1.0, -1.0) + return addVectors([v1, v1p], [1.0, -1.0]) -def scalarProduct(s, v): +def scaleVector(v, s): """ Calculate s * v - :param s: Scalar. :param v: Vector. + :param s: Scalar. :return: """ - return [s * v[c] for c in range(len(v))] + return [s * c for c in v] def parallelVectors(v1, v2): @@ -88,3 +96,21 @@ def parallelVectors(v1, v2): if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2)): return True return False + + +def angleBetweenVectors(v1, v2): + """ + :return: Angle between vectors v1 and v2 in radians + """ + return math.acos(dotproduct(normalise(v1), normalise(v2))) + + +def rotateVectorAroundVector(v, k, a): + """ + Rotate vector v, by an angle a (right-hand rule) in radians around vector k. + :return: rotated vector. + """ + k = normalise(k) + vperp = addVectors([v, crossproduct3(k, v)], [math.cos(a), math.sin(a)]) + vparal = scaleVector(k, dotproduct(k, v)*(1 - math.cos(a))) + return addVectors([vperp, vparal]) diff --git a/tests/test_sphere.py b/tests/test_sphere.py new file mode 100644 index 00000000..7609026f --- /dev/null +++ b/tests/test_sphere.py @@ -0,0 +1,197 @@ +import unittest +from testutils import assertAlmostEqualList + +from opencmiss.zinc.context import Context +from opencmiss.zinc.field import Field +from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetRange +from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.result import RESULT_OK +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, getAnnotationGroupForTerm +from scaffoldmaker.meshtypes.meshtype_3d_solidsphere2 import MeshType_3d_solidsphere2 +from scaffoldmaker.utils.meshrefinement import MeshRefinement + + +class SphereScaffoldTestCase(unittest.TestCase): + + def test_sphere1(self): + """ + Test creation of Sphere scaffold. + """ + scaffold = MeshType_3d_solidsphere2 + parameterSetNames = scaffold.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default"]) + options = scaffold.getDefaultOptions("Default") + self.assertEqual(16, len(options)) + self.assertEqual(4, options.get("Number of elements across axis 1")) + self.assertEqual(4, options.get("Number of elements across axis 2")) + self.assertEqual(4, options.get("Number of elements across axis 3")) + self.assertEqual(0, options.get("Number of elements across shell")) + self.assertEqual(1, options.get("Number of elements across transition")) + self.assertEqual(1.0, options.get("Radius1")) + self.assertEqual(1.0, options.get("Radius2")) + self.assertEqual(1.0, options.get("Radius3")) + self.assertEqual(1.0, options.get("Shell element thickness proportion")) + self.assertEqual([0, 4], options.get("Range of elements required in direction 1")) + self.assertEqual([0, 4], options.get("Range of elements required in direction 2")) + self.assertEqual([0, 4], options.get("Range of elements required in direction 3")) + self.assertEqual([1, 2, 3], options.get("Box derivatives")) + + context = Context("Test") + region = context.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold.generateMesh(region, options) + self.assertEqual(2, len(annotationGroups)) + + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(32, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(108, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(128, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(53, nodes.getSize()) + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # check coordinates range, sphere volume + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-1.0, -1.0, -1.0], 1.0E-6) + assertAlmostEqualList(self, maximums, [1.0, 1.0, 1.0], 1.0E-6) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + surfaceGroup = AnnotationGroup(region, ("sphere surface", "")) + is_exterior = fieldmodule.createFieldIsExterior() + surfaceMeshGroup = surfaceGroup.getMeshGroup(mesh2d) + surfaceMeshGroup.addElementsConditional(is_exterior) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 12.460033954564986, delta=2.0E-1) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 4.132033912594377, delta=1.0E-3) + + # check some annotationGroups: + expectedSizes3d = { + "box group": 8, + "transition group": 24, + } + for name in expectedSizes3d: + group = getAnnotationGroupForTerm(annotationGroups, (name, '')) + size = group.getMeshGroup(mesh3d).getSize() + self.assertEqual(expectedSizes3d[name], size, name) + + # refine 8x8x8 and check result + refineRegion = region.createRegion() + refineFieldmodule = refineRegion.getFieldmodule() + options['Refine number of elements'] = 2 + meshrefinement = MeshRefinement(region, refineRegion, []) + scaffold.refineMesh(meshrefinement, options) + + refineCoordinates = refineFieldmodule.findFieldByName("coordinates").castFiniteElement() + + refineFieldmodule.defineAllFaces() + mesh3d = refineFieldmodule.findMeshByDimension(3) + self.assertEqual(256, mesh3d.getSize()) + mesh2d = refineFieldmodule.findMeshByDimension(2) + self.assertEqual(816, mesh2d.getSize()) + mesh1d = refineFieldmodule.findMeshByDimension(1) + self.assertEqual(880, mesh1d.getSize()) + nodes = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(321, nodes.getSize()) + datapoints = refineFieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # obtain surface area and volume again. If they are the same as previous, mesh is conform. + with ChangeManager(refineFieldmodule): + one = refineFieldmodule.createFieldConstant(1.0) + surfaceGroup = AnnotationGroup(refineRegion, ("sphere surface", "")) + is_exterior = refineFieldmodule.createFieldIsExterior() + surfaceMeshGroup = surfaceGroup.getMeshGroup(mesh2d) + surfaceMeshGroup.addElementsConditional(is_exterior) + + surfaceAreaField = refineFieldmodule.createFieldMeshIntegral(one, refineCoordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + volumeField = refineFieldmodule.createFieldMeshIntegral(one, refineCoordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + refineFieldcache = refineFieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(refineFieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 12.460033954564986, delta=5.0E-1) + result, volume = volumeField.evaluateReal(refineFieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 4.132033912594377, delta=3.0E-1) + + # check ellipsoid. + scaffold1 = MeshType_3d_solidsphere2 + parameterSetNames = scaffold1.getParameterSetNames() + self.assertEqual(parameterSetNames, ["Default"]) + options = scaffold1.getDefaultOptions("Default") + options['Number of elements across axis 1'] = 4 + options['Number of elements across axis 2'] = 6 + options['Number of elements across axis 3'] = 8 + options['Range of elements required in direction 1'] = [0, 4] + options['Range of elements required in direction 2'] = [0, 6] + options['Range of elements required in direction 3'] = [0, 8] + + options['Radius1'] = 0.5 + options['Radius2'] = 0.8 + options['Radius3'] = 1.0 + + context2 = Context("Test2") + region = context2.getDefaultRegion() + self.assertTrue(region.isValid()) + annotationGroups = scaffold1.generateMesh(region, options) + self.assertEqual(2, len(annotationGroups)) + fieldmodule = region.getFieldmodule() + self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) + mesh3d = fieldmodule.findMeshByDimension(3) + self.assertEqual(136, mesh3d.getSize()) + mesh2d = fieldmodule.findMeshByDimension(2) + self.assertEqual(452, mesh2d.getSize()) + mesh1d = fieldmodule.findMeshByDimension(1) + self.assertEqual(510, mesh1d.getSize()) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self.assertEqual(195, nodes.getSize()) + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self.assertEqual(0, datapoints.getSize()) + + # check coordinates range, sphere volume + coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() + self.assertTrue(coordinates.isValid()) + minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) + assertAlmostEqualList(self, minimums, [-0.5, -0.8, -1.0], 1.0E-6) + assertAlmostEqualList(self, maximums, [0.5, 0.8, 1.0], 1.0E-6) + + with ChangeManager(fieldmodule): + one = fieldmodule.createFieldConstant(1.0) + surfaceGroup = AnnotationGroup(region, ("sphere surface", "")) + is_exterior = fieldmodule.createFieldIsExterior() + surfaceMeshGroup = surfaceGroup.getMeshGroup(mesh2d) + surfaceMeshGroup.addElementsConditional(is_exterior) + + surfaceAreaField = fieldmodule.createFieldMeshIntegral(one, coordinates, surfaceMeshGroup) + surfaceAreaField.setNumbersOfPoints(4) + volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d) + volumeField.setNumbersOfPoints(3) + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 7.3015173697377245, delta=2.0E-1) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 1.6741674010981926, delta=1.0E-3) + + +if __name__ == "__main__": + unittest.main()