diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py index 5d59cf1c..c8ccb322 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_solidsphere2.py @@ -8,9 +8,9 @@ from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates from scaffoldmaker.annotation.annotationgroup import AnnotationGroup from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils import vector from scaffoldmaker.utils.meshrefinement import MeshRefinement from scaffoldmaker.utils.spheremesh import SphereMesh, SphereShape -from scaffoldmaker.utils import vector class MeshType_3d_solidsphere2(Scaffold_base): @@ -30,11 +30,11 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Number of elements across axis 2': 4, 'Number of elements across axis 3': 4, 'Number of elements across shell': 0, + 'Shell element thickness proportion': 1.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], @@ -51,12 +51,12 @@ def getOrderedOptionNames(): '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 shell', + 'Shell element thickness proportion', # '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', @@ -96,14 +96,16 @@ def checkOptions(cls, options): 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: + if options['Range of elements required in direction {}'.format(i)][1] >= \ + options['Number of elements across axis {}'.format(i)] - 1\ + - options['Number of elements across shell']: 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: + if options['Range of elements required in direction {}'.format(i)][0] <= 1 + \ + options['Number of elements across shell']: options['Range of elements required in direction {}'.format(i)][0] = 0 maxelems = [options['Number of elements across axis 1'], @@ -114,25 +116,29 @@ def checkOptions(cls, options): 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]): + dependentChanges = True 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: + dependentChanges = True 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 + elementsCount_min = min(options['Number of elements across axis 1'], + options['Number of elements across axis 2'], + options['Number of elements across axis 3'],) + + max_shell_elements = elementsCount_min//2 - 2 + + if options['Number of elements across shell'] > max_shell_elements: + dependentChanges = True + options['Number of elements across shell'] = max_shell_elements + + if options['Shell element thickness proportion'] < 0.15: + options['Shell element thickness proportion'] = 1.0 return dependentChanges diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index 350c312d..10b343f1 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -83,9 +83,9 @@ def __init__(self): MeshType_3d_lung1, MeshType_3d_ostium1, MeshType_3d_smallintestine1, + MeshType_3d_solidcylinder1, MeshType_3d_solidsphere1, MeshType_3d_solidsphere2, - MeshType_3d_solidcylinder1, MeshType_3d_sphereshell1, MeshType_3d_sphereshellseptum1, MeshType_3d_stellate1, diff --git a/src/scaffoldmaker/utils/shieldmesh.py b/src/scaffoldmaker/utils/shieldmesh.py index b202f833..9b75dd74 100644 --- a/src/scaffoldmaker/utils/shieldmesh.py +++ b/src/scaffoldmaker/utils/shieldmesh.py @@ -788,7 +788,7 @@ def __init__(self, elementsCountAcross, elementsCountRim, box_derivatives=None): self.elementsCountAcross = elementsCountAcross self.elementsCountRim = elementsCountRim # self.elementsCountUpRegular = elementsCountUp - 2 - elementsCountRim - # elementsCountAcrossNonRim = self.elementsCountAcross - 2*elementsCountRim + self.elementsCountAcrossNonRim = [ne - 2*elementsCountRim for ne in elementsCountAcross] # self.elementsCountAroundFull = 2*self.elementsCountUpRegular + elementsCountAcrossNonRim self._boxDerivatives = box_derivatives self._boxMapping = None @@ -879,7 +879,7 @@ def generateNodes(self, fieldmodule, coordinates, startNodeIdentifier, rangeOfRe :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. + :param rangeOfRequiredElements: Only the elements and nodes for the given range is generated. :return: next nodeIdentifier. """ nodeIdentifier = startNodeIdentifier @@ -1123,9 +1123,13 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, ran scalefactors = None self._element_needs_scale_factor = False + def shell_element(e3o, e2o, e1o, e3zo, e1zo): + return e3o > e3zo or e2o < self.elementsCountRim or e1o > e1zo + 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 + element_shell = shell_element(e3o, e2o, e1o, e3zo, e1zo) + e3yo, e2bo, e1yo = e3zo - 1, 1 + self.elementsCountRim, e1zo - 1 eft1 = eft if element_type == self.ELEMENT_REGULAR else\ tricubichermite.createEftNoCrossDerivatives() @@ -1143,87 +1147,168 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, ran 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.QUADRUPLE_DOWN_LEFT) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_12_LEFT, + derivatives=triple12leftderivs) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE_3_LEFT) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_13_DOWN) + else: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE0_13_DOWN) + if e3yo == 0: + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_1, + derivatives=corner1derivs) 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_13_DOWN) else: - self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE0_2_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE_2_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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) else: + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.SURFACE_REGULAR_DOWN_LEFT) 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) + if not element_shell: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE_2_DOWN) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_1, + derivatives=corner1derivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_13_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.SURFACE_REGULAR_DOWN_LEFT) else: self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_2_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) else: + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2], lnm[6]], + self.SURFACE_REGULAR_DOWN_LEFT) 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.QUADRUPLE_RIGHT) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_12_RIGHT, + derivatives=triple12rightderivs) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_3_RIGHT) + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_23_DOWN) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_2, + derivatives=corner2derivs) else: self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_23_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_23_DOWN) else: - self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE0_1_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_CURVE_1_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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) else: + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.SURFACE_REGULAR_DOWN_RIGHT) 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_1_DOWN) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.TRIPLE_23_DOWN) + else: + 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: @@ -1231,157 +1316,301 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, ran derivatives=boundary12rightderivs) self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, derivatives=corner2derivs) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_23_DOWN) else: - self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.TRIPLE_CURVE0_1_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4]], self.TRIPLE_CURVE_1_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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) else: self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.SURFACE_REGULAR_DOWN_RIGHT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2]], + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.QUADRUPLE_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.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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_23_UP) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_3, + derivatives=corner3derivs) + self.remap_eft_node_value_label(eft1, [lnm[1]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_2_UP) else: - self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE0_1_UP) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE_1_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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_1_UP) + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_23_UP) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_3, + derivatives=corner3derivs) + self.remap_eft_node_value_label(eft1, [lnm[1]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_23_UP) + self.remap_eft_node_value_label(eft1, [lnm[1]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.TRIPLE_CURVE_1_UP) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.TRIPLE_CURVE_1_UP) + else: + 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) + if not element_shell: + 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) + if not element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE0_2_UP) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.TRIPLE_CURVE_2_UP) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_2_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], self.SURFACE_REGULAR_UP) self.remap_eft_node_value_label(eft1, [lnm[8]], self.SURFACE_REGULAR_UP) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[5]], self.SURFACE_REGULAR_UP) + self.remap_eft_node_value_label(eft1, [lnm[6]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2]], self.BOUNDARY_13_UP) + self.remap_eft_node_value_label(eft1, [lnm[5], lnm[6]], self.SURFACE_REGULAR_UP) else: self.remap_eft_node_value_label(eft1, [lnm[3], lnm[4], lnm[7], lnm[8]], self.SURFACE_REGULAR_UP) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1], lnm[2], lnm[5], lnm[6]], + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[8]], self.TRIPLE_CURVE_3_LEFT) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_12_LEFT, + derivatives=triple12leftderivs) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.TRIPLE_CURVE_3_LEFT) + else: + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_13_DOWN) if e3o == 0: self.remap_eft_node_value_label(eft1, [lnm[1]], self.CORNER_1, derivatives=corner1derivs) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_1, + derivatives=corner1derivs) else: self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_13_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_13_DOWN) else: self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_LEFT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.SURFACE_REGULAR_DOWN_LEFT) if e3o == 0: self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_LEFT, derivatives=boundary12leftderivs) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_LEFT, + derivatives=boundary12leftderivs) else: self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_DOWN_LEFT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.TRIPLE_CURVE_3_RIGHT) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_12_RIGHT, + derivatives=triple12rightderivs) + else: + 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[1]], self.TRIPLE_CURVE_3_RIGHT) + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.BOUNDARY_23_DOWN) if e3o == 0: self.remap_eft_node_value_label(eft1, [lnm[6]], self.CORNER_2, derivatives=corner2derivs) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_2, + derivatives=corner2derivs) else: self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_23_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_23_DOWN) else: self.remap_eft_node_value_label(eft1, [lnm[8]], self.SURFACE_REGULAR_DOWN_RIGHT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[4]], self.SURFACE_REGULAR_DOWN_RIGHT) if e3o == 0: self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_RIGHT, derivatives=boundary12rightderivs) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) else: self.remap_eft_node_value_label(eft1, [lnm[6]], self.SURFACE_REGULAR_DOWN_RIGHT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.SURFACE_REGULAR_DOWN_RIGHT) elif element_type == self.ELEMENT_DOWN_RIGHT: if e3o == 0: @@ -1392,6 +1621,13 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, ran 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_2, + derivatives=corner2derivs) + self.remap_eft_node_value_label(eft1, [lnm[4]], 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, @@ -1399,33 +1635,67 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, ran 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[1]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_RIGHT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[4]], 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[3]], self.SURFACE_REGULAR_DOWN_RIGHT) + self.remap_eft_node_value_label(eft1, [lnm[1]], self.SURFACE_REGULAR_DOWN_RIGHT) if e2o == e2zo: self.remap_eft_node_value_label(eft1, [lnm[6], lnm[8]], self.BOUNDARY_23_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2], lnm[4]], self.BOUNDARY_23_DOWN) else: self.remap_eft_node_value_label(eft1, [lnm[6], lnm[8]], self.SURFACE_REGULAR_DOWN_RIGHT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2], lnm[4]], + 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 element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6]], self.BOUNDARY_12_LEFT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[8]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.CORNER_1, + derivatives=corner1derivs) + self.remap_eft_node_value_label(eft1, [lnm[4]], 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) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2]], self.BOUNDARY_12_LEFT, + derivatives=boundary12rightderivs) + self.remap_eft_node_value_label(eft1, [lnm[4]], self.SURFACE_REGULAR_DOWN_LEFT) else: self.remap_eft_node_value_label(eft1, [lnm[5], lnm[7]], self.SURFACE_REGULAR_DOWN_LEFT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[6], lnm[8]], self.SURFACE_REGULAR_DOWN_LEFT) if e1o == 0: self.remap_eft_node_value_label(eft1, [lnm[1], lnm[3]], self.BOUNDARY_13_DOWN) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2], lnm[4]], self.BOUNDARY_13_DOWN) else: self.remap_eft_node_value_label(eft1, [lnm[1], lnm[3]], self.SURFACE_REGULAR_DOWN_LEFT) + if element_shell: + self.remap_eft_node_value_label(eft1, [lnm[2], lnm[4]], + self.SURFACE_REGULAR_DOWN_LEFT) else: continue @@ -1547,36 +1817,36 @@ def get_element_type(self, e3, e2, e1): octant_number = 8 if octant_number in [2, 4, 5, 7]: - e3zo = self.elementsCountAcross[2] // 2 - 1 + e3zo = self.elementsCountAcross[2] // 2 - 1 - self.elementsCountRim e2zo = self.elementsCountAcross[1] // 2 - 1 - e1zo = self.elementsCountAcross[0] // 2 - 1 + e1zo = self.elementsCountAcross[0] // 2 - 1 - self.elementsCountRim else: - e3zo = self.elementsCountAcross[2] // 2 - 1 + e3zo = self.elementsCountAcross[2] // 2 - 1 - self.elementsCountRim e2zo = self.elementsCountAcross[0] // 2 - 1 - e1zo = self.elementsCountAcross[1] // 2 - 1 + e1zo = self.elementsCountAcross[1] // 2 - 1 - self.elementsCountRim - e3yo, e2bo, e1yo = e3zo - 1, 1, e1zo - 1 + e3yo, e2bo, e1yo = e3zo - 1, 1 + self.elementsCountRim, 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: + elif e3o == e3yo and e2o <= self.elementsCountRim and e1o == e1yo: element_type = self.ELEMENT_QUADRUPLE_DOWN_LEFT - elif e3o == e3yo and e2o >= e2bo and e1o == e1zo: + 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: + 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: + elif e3o == e3yo and e2o <= self.elementsCountRim and e1o < e1yo: element_type = self.ELEMENT_QUADRUPLE_DOWN - elif e3o == e3zo and e2o >= e2bo and e1o < e1yo: + elif e3o >= e3zo and e2o >= e2bo and e1o < e1yo: element_type = self.ELEMENT_QUADRUPLE_UP - elif e3o < e3yo and e2o == 0 and e1o == e1yo: + elif e3o < e3yo and e2o <= self.elementsCountRim and e1o == e1yo: element_type = self.ELEMENT_QUADRUPLE_LEFT - elif e3o < e3yo and e2o == e2bo and e1o == e1zo: + elif e3o < e3yo and e2o == e2bo and e1o >= e1zo: element_type = self.ELEMENT_QUADRUPLE_RIGHT - elif e3o < e3yo and e2o > e2bo and e1o == e1zo: + elif e3o < e3yo and e2o > e2bo and e1o >= e1zo: element_type = self.ELEMENT_DOWN_RIGHT - elif e3o < e3yo and e2o == 0 and e1o < e1yo: + elif e3o < e3yo and e2o <= self.elementsCountRim and e1o < e1yo: element_type = self.ELEMENT_DOWN_LEFT else: element_type = 0 @@ -1676,30 +1946,28 @@ def getNodeId(self, octant_number, n3, n2, n1, n3yo, n1yo): 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 + + n3_skip, n2_skip, n1_skip = 0, 0, 0 + + if n1o == n1yo: + if n2o <= self.elementsCountRim: + n1_skip = 1 + self.elementsCountRim - n2o + elif n3o >= n3zo: + n1_skip = 1 + n3o - n3zo + + if n2o == self.elementsCountRim + 1: + if n1o >= n1zo: + n2_skip = -(1 + n1o - n1zo) + elif n3o >= n3zo: + n2_skip = -(1 + n3o - n3zo) + + if n3o == n3yo: + if n2o <= self.elementsCountRim: + n3_skip = 1 + self.elementsCountRim - n2o + elif n1o >= n1zo: + n3_skip = 1 + n1o - n1zo + + n3r, n2r, n1r = self.get_global_node_index(octant_number, n3o + n3_skip, n2o + n2_skip, n1o + n1_skip) return self.nodeId[n3r][n2r][n1r] diff --git a/src/scaffoldmaker/utils/spheremesh.py b/src/scaffoldmaker/utils/spheremesh.py index bd0d4125..7bdd7c29 100644 --- a/src/scaffoldmaker/utils/spheremesh.py +++ b/src/scaffoldmaker/utils/spheremesh.py @@ -60,7 +60,8 @@ def __init__(self, fieldModule, coordinates, centre, axes, elementsCountAcross, self._radius = [vector.magnitude(axis) for axis in axes] self._coreRadius = [] self._shield = None - self._elementsCount = elementsCountAcross + self._elementsCount = [ec - 2 * elementsCountAcrossShell for ec in elementsCountAcross] + self._elementsCount_total = elementsCountAcross self._elementsCountAcrossShell = elementsCountAcrossShell self._elementsCountAcrossTransition = elementsCountAcrossTransition self._elementsCountAcrossRim = self._elementsCountAcrossShell + self._elementsCountAcrossTransition - 1 @@ -82,9 +83,10 @@ def __init__(self, fieldModule, coordinates, centre, axes, elementsCountAcross, self._shield3D = None for i in range(3): - elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - shellProportion) + # Make a sphere of radius 1 first, then add the shell layers on top of it and resize it back to R=1.0 + elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - 1.0) self._coreRadius.append( - (1 - shellProportion * elementsCountAcrossShell / elementsAxis) * self._radius[i]) + (1 - 1.0 * elementsCountAcrossShell / elementsAxis) * self._radius[i]) # generate the mesh self.createSphereMesh3d(fieldModule, coordinates) @@ -125,9 +127,7 @@ def createSphereMesh3d(self, fieldModule, coordinates): 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, @@ -142,15 +142,18 @@ def createSphereMesh3d(self, fieldModule, coordinates): else: OctantVariationsType = OctantVariationsAll + # Make a sphere of radius 1 first, then add the shell layers on top of it and resize it back to R=1.0 for octantType in OctantVariationsType: axes, elementsCountAcross, boxDerivatives = self.get_octant_axes_and_elements_count(octantType) octant = OctantMesh(self._centre, axes, elementsCountAcross, - elementsCountAcrossShell, elementsCountAcrossTransition, shellProportion, + 0, elementsCountAcrossTransition, 1.0, sphereShape=SphereShape.SPHERESHIELD_SHAPE_OCTANT_PPP, useCrossDerivatives=False, boxDerivatives=boxDerivatives) self.copy_octant_nodes_to_sphere_shield(octant, octantType) self.modify_octant_common_nodes() + if self._elementsCountAcrossShell > 0: + self.add_shell_layers() self.sphere_to_spheroid() self.generateNodes(nodes, fieldModule, coordinates) @@ -329,24 +332,100 @@ def modify_octant_common_nodes(self): 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): + def add_shell_layers(self): """ - Using the radius in each direction,transform the sphere to ellipsoid. + Add shell layers on top of the sphere. Calculates shell element thickness using shell proportion and adds + shell elements. Note that Radius without the shell elements is 1.0. R= 1.0 + thickness X number of shell layers. :return: """ + + def on_sphere(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] + n2z = self._elementsCount[0] + n1z = self._elementsCount[1] + + x = self._shield3D.px + + return (n3 == n3z or n3 == 0 or n2 == 0 or n2 == n2z or n1 == 0 or n1 == n1z) and x[n3][n2][n1] + + shield3D_with_shell = ShieldMesh3D(self._elementsCount_total, self._elementsCountAcrossShell, + box_derivatives=self._boxDerivatives) + btx = self._shield3D.px btd1 = self._shield3D.pd1 btd2 = self._shield3D.pd2 btd3 = self._shield3D.pd3 + shell_proportion = self._shellProportion + thickness = [2*1.0/ne * shell_proportion for ne in self._elementsCount] + elementsCountShell = self._elementsCountAcrossShell + for n3 in range(self._elementsCount[2] + 1): for n2 in range(self._elementsCount[0] + 1): for n1 in range(self._elementsCount[1] + 1): + # new numbering due to introducing shell layers + n3n, n2n, n1n = n3 + elementsCountShell, n2 + elementsCountShell, n1 + elementsCountShell + # generate rim/shell nodes. Note, sphere radius without shell is 1.0. To keep it 1.0 after adding + # shell elements, we need to shrink it back. See ratio in sphere_to_spheroid. + if on_sphere(n3, n2, n1): + for n in range(elementsCountShell): + n3s, n2s, n1s = n3n, n2n, n1n + nl = n + 1 + if n3 == 0: + n3s = elementsCountShell - nl + if n2 == 0: + n2s = elementsCountShell - nl + if n1 == 0: + n1s = elementsCountShell - nl + if n3 == self._elementsCount[2]: + n3s = n3n + nl + if n2 == self._elementsCount[0]: + n2s = n2n + nl + if n1 == self._elementsCount[1]: + n1s = n1n + nl + + shield3D_with_shell.px[n3s][n2s][n1s] =\ + [(1 + nl * thickness[c]) * btx[n3][n2][n1][c] for c in range(3)] + shield3D_with_shell.pd1[n3s][n2s][n1s] = \ + [(1 + nl * thickness[c]) * btd1[n3][n2][n1][c] for c in range(3)] + shield3D_with_shell.pd2[n3s][n2s][n1s] =\ + [(1 + nl * thickness[c]) * btd2[n3][n2][n1][c] for c in range(3)] + shield3D_with_shell.pd3[n3s][n2s][n1s] =\ + [shell_proportion * btd3[n3][n2][n1][c] for c in range(3)] + # Add the rest of the nodes. + if btx[n3][n2][n1]: + shield3D_with_shell.px[n3n][n2n][n1n] = btx[n3][n2][n1] + shield3D_with_shell.pd1[n3n][n2n][n1n] = btd1[n3][n2][n1] + shield3D_with_shell.pd2[n3n][n2n][n1n] = btd2[n3][n2][n1] + shield3D_with_shell.pd3[n3n][n2n][n1n] = btd3[n3][n2][n1] + + self._shield3D = shield3D_with_shell + + 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 + + ratio = [1/(1 + 2/ne * self._shellProportion * self._elementsCountAcrossShell) for ne in self._elementsCount] + + for n3 in range(self._elementsCount_total[2] + 1): + for n2 in range(self._elementsCount_total[0] + 1): + for n1 in range(self._elementsCount_total[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)] + btx[n3][n2][n1] = [ratio[c] * self._radius[c] * btx[n3][n2][n1][c] for c in range(3)] + btd1[n3][n2][n1] = [ratio[c] * self._radius[c] * btd1[n3][n2][n1][c] for c in range(3)] + btd2[n3][n2][n1] = [ratio[c] * self._radius[c] * btd2[n3][n2][n1][c] for c in range(3)] + btd3[n3][n2][n1] = [ratio[c] * self._radius[c] * btd3[n3][n2][n1][c] for c in range(3)] def generateNodes(self, nodes, fieldModule, coordinates): """ @@ -414,9 +493,9 @@ def __init__(self, centre, axes, elementsCountAcross, self._shield3D = None for i in range(3): - elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - shellProportion) + elementsAxis = elementsCountAcross[i] - elementsCountAcrossShell * (1 - 1.0) self._coreRadius.append( - (1 - shellProportion * elementsCountAcrossShell / elementsAxis) * self._radius[i]) + (1 - 1.0 * elementsCountAcrossShell / elementsAxis) * self._radius[i]) # generate the mesh self.createOctantMesh3d()