From eda7f10e127da246332b9de4536feecf4b4a32e1 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 13 Nov 2020 17:50:01 +1300 Subject: [PATCH 01/23] Added lungs meshtypes. --- src/scaffoldmaker/scaffolds.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index b523cde6..dce08938 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -29,6 +29,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase1 import MeshType_3d_heartventriclesbase1 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase2 import MeshType_3d_heartventriclesbase2 from scaffoldmaker.meshtypes.meshtype_3d_lens1 import MeshType_3d_lens1 +from scaffoldmaker.meshtypes.meshtype_3d_lungs1 import MeshType_3d_lungs1 from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1 from scaffoldmaker.meshtypes.meshtype_3d_smallintestine1 import MeshType_3d_smallintestine1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 @@ -72,6 +73,7 @@ def __init__(self): MeshType_3d_heartventriclesbase1, MeshType_3d_heartventriclesbase2, MeshType_3d_lens1, + MeshType_3d_lungs1, MeshType_3d_ostium1, MeshType_3d_smallintestine1, MeshType_3d_solidsphere1, From 4fc2e2f67fb9d33af60652560606ab17324e6ec6 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 13 Nov 2020 17:58:26 +1300 Subject: [PATCH 02/23] Created regular elements for the left lung. --- .../meshtypes/meshtype_3d_lungs1.py | 335 ++++++++++++++++++ 1 file changed, 335 insertions(+) create mode 100644 src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py new file mode 100644 index 00000000..acd749b2 --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py @@ -0,0 +1,335 @@ +''' +Generates 3D lung surface mesh, with variable numbers of elements around, along. +''' + +import copy +import math +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm +from scaffoldmaker.annotation.bladder_terms import get_bladder_term +from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.scaffoldpackage import ScaffoldPackage +from scaffoldmaker.utils import interpolation as interp +from scaffoldmaker.utils import matrix +from scaffoldmaker.utils import tubemesh +from scaffoldmaker.utils import vector +from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d +from scaffoldmaker.utils.geometry import createEllipsePoints, createEllipsoidPoints +from scaffoldmaker.utils.interpolation import smoothCubicHermiteDerivativesLine +from scaffoldmaker.utils.tracksurface import TrackSurface, TrackSurfacePosition, calculate_surface_axes +from scaffoldmaker.utils.zinc_utils import exnodeStringFromNodeValues, mesh_destroy_elements_and_nodes_by_identifiers +from opencmiss.zinc.element import Element, Elementbasis +from opencmiss.zinc.field import Field +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates +from opencmiss.zinc.node import Node + + +class MeshType_3d_lungs1(Scaffold_base): + ''' + 3D lung scaffold. + ''' + + @staticmethod + def getName(): + return '3D Lungs 1' + + @staticmethod + def getParameterSetNames(): + return [ + 'Default', + 'Mouse 1', + 'Human 1'] + + @classmethod + def getDefaultOptions(cls, parameterSetName='Default'): + options = { + 'Number of elements around': 8, + 'Number of elements along': 4, + 'Number of elements through wall': 1, + 'Major diameter': 30.0, + 'Minor diameter': 15.0, + 'Height': 70, + 'Use cross derivatives': False, + 'Refine': False, + 'Refine number of elements along': 4, + 'Refine number of elements around': 4 + } + if 'Mouse' in parameterSetName: + options['Major diameter'] = 20.0 + options['Minor diameter'] = 10.0 + return options + + @staticmethod + def getOrderedOptionNames(): + optionNames = [ + 'Number of elements around', + 'Number of elements along', + 'Major diameter', + 'Minor diameter', + 'Height', + 'Use cross derivatives', + 'Refine', + 'Refine number of elements around', + 'Refine number of elements along'] + return optionNames + + + # @classmethod + # def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType): + # if optionName == 'Central path LUT': + # return list(cls.centralPathDefaultScaffoldPackages_LUT.keys()) + # assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), cls.__name__ + '.getOptionScaffoldTypeParameterSetNames. ' + \ + # 'Invalid option \'' + optionName + '\' scaffold type ' + scaffoldType.getName() + # return scaffoldType.getParameterSetNames() + + @classmethod + def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): + ''' + :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. + :return: ScaffoldPackage. + ''' + if parameterSetName: + assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ + 'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName() + assert False, cls.__name__ + '.getOptionScaffoldPackage: Option ' + optionName + ' is not a scaffold' + + @classmethod + def checkOptions(cls, options): + for key in [ + 'Number of elements around', + 'Number of elements along', + 'Refine number of elements along', + 'Refine number of elements around']: + if options[key] < 1: + options[key] = 1 + + @classmethod + def generateBaseMesh(cls, region, options): + ''' + Generate the base tricubic Hermite mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: annotationGroups + ''' + + useCrossDerivatives = options['Use cross derivatives'] + + fm = region.getFieldmodule() + fm.beginChange() + coordinates = findOrCreateFieldCoordinates(fm) + + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplateApex = nodes.createNodetemplate() + nodetemplateApex.defineField(coordinates) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + if useCrossDerivatives: + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + else: + nodetemplate = nodetemplateApex + + + cache = fm.createFieldcache() + + + # Create nodes + nodesList = [[[-1.0, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + + [[-1.0, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + + [[-1.0, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] + # + # [[0.0, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # + # [[0.0, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[0.0, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-0.5, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + # [[-1.0, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] + + # [[0.0, 0.0, 4.7], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.5, 0.0]]] + + print('len(nodesList)', len(nodesList)) + + cx = [] + cd1 = [] + cd2 = [] + cd3 = [] + nodeIdentifier = 1 + for n1 in range(len(nodesList)): + x = nodesList[n1][0] + dx_ds1 = nodesList[n1][1] + dx_ds2 = nodesList[n1][2] + dx_ds3 = nodesList[n1][3] + cx.append(x) + cd1.append(dx_ds1) + cd2.append(dx_ds2) + cd3.append(dx_ds3) + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cx[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, cd1[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, cd2[n1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, cd3[n1]) + nodeIdentifier += 1 + + print('len(cx)', len(cx)) + print('len(cd1)', len(cd1)) + print('len(cd2)', len(cd2)) + + + # Create elements + from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite + + mesh = fm.findMeshByDimension(3) + eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives) + eft = eftfactory.createEftBasic() + + elementtemplate = mesh.createElementtemplate() + elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplate.defineField(coordinates, -1, eft) + + # elementtemplate2 = mesh.createElementtemplate() + # elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) + # eft2 = eftfactory.createEftWedgeXi1Zero() + # elementtemplate2.defineField(coordinates, -1, eft2) + + elementIdentifier = 1 + elementsCountAlongMajorDiameter = 4 + elementsCountAlongMinorDiameter = 2 + elementsCountAlong = 2 + # + nel = (elementsCountAlongMinorDiameter + 1) * (elementsCountAlongMajorDiameter + 1) + for e3 in range(elementsCountAlong): + for e2 in range(elementsCountAlongMajorDiameter): + # First row of elements #4 + for e1 in range(elementsCountAlongMinorDiameter): + if (e2 == 0 or e2 == elementsCountAlongMajorDiameter - 1) and e1 == 0: + # nodeIdentifiers = [10, 11, 14, 25, 26, 29] + # element = mesh.createElement(elementIdentifier, elementtemplate2) + # element.setNodesByIdentifier(eft2, nodeIdentifiers) + pass + else: + element = mesh.createElement(elementIdentifier, elementtemplate) + bni1 = e3 * nel + e2 * (elementsCountAlongMinorDiameter + 1) + e1 + 1 + bni2 = e3 * nel + e2 * (elementsCountAlongMinorDiameter + 1) + (e1 + 1) % (elementsCountAlongMinorDiameter + 1) + 1 + bni3 = bni1 + elementsCountAlongMinorDiameter + 1 + bni4 = bni3 + 1 + bni5 = bni1 + nel + bni6 = bni5 + 1 + bni7 = bni5 + elementsCountAlongMinorDiameter + 1 + bni8 = bni7 + 1 + nodeIdentifiers = [bni1, bni2, bni3, bni4, + bni5, bni6, bni7, bni8] + result = element.setNodesByIdentifier(eft, nodeIdentifiers) + elementIdentifier += 1 + + annotationGroups = [] + fm.endChange() + 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(). + ''' + refineElementsCountAlong = options['Refine number of elements along'] + refineElementsCountAround = options['Refine number of elements around'] + + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong) + return + + # def createEftWedgeXi1Zero(self): + # ''' + # Create a basic tricubic hermite element template for elements + # along boundary of tenia coli where nodes on xi1 = 0 are collapsed. + # :return: Element field template + # ''' + # eft = self.createEftBasic() + # + # # remap parameters on xi1 = 1 before collapsing nodes + # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS2, []) + # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS1DS2, []) + # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS1DS2, []) + # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + # + # ln_map = [1, 2, 1, 3, 4, 5, 4, 6] + # remapEftLocalNodes(eft, 6, ln_map) + # assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1Zero: Failed to validate eft' + # return eft + From 6e8f67a72a2a5e6f9e5d99bce0df31a58a90f229 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Mon, 16 Nov 2020 19:40:40 +1300 Subject: [PATCH 03/23] Created the back and front wedge elements. --- .../meshtypes/meshtype_3d_lungs1.py | 326 +++++++++++------- 1 file changed, 205 insertions(+), 121 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py index acd749b2..f604453b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py @@ -73,15 +73,6 @@ def getOrderedOptionNames(): 'Refine number of elements along'] return optionNames - - # @classmethod - # def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType): - # if optionName == 'Central path LUT': - # return list(cls.centralPathDefaultScaffoldPackages_LUT.keys()) - # assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), cls.__name__ + '.getOptionScaffoldTypeParameterSetNames. ' + \ - # 'Invalid option \'' + optionName + '\' scaffold type ' + scaffoldType.getName() - # return scaffoldType.getParameterSetNames() - @classmethod def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): ''' @@ -138,6 +129,9 @@ def generateBaseMesh(cls, region, options): cache = fm.createFieldcache() + elementsCount1 = 2 + elementsCount2 = 4 + elementsCount3 = 4 # Create nodes nodesList = [[[-1.0, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], @@ -186,69 +180,90 @@ def generateBaseMesh(cls, region, options): [[0.0, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], [[-1.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], [[-0.5, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] - # - # [[0.0, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, -1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, -0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 0.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 0.5, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 1.0, 2.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # - # [[0.0, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, -1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, -0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 0.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 0.5, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[0.0, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-0.5, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - # [[-1.0, 1.0, 3.0], [-0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] - - # [[0.0, 0.0, 4.7], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.5, 0.0]]] - + [[0.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + + [[-1.0, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + + [[-1.0, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-1.0, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[-0.5, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], + [[0.0, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] print('len(nodesList)', len(nodesList)) - cx = [] - cd1 = [] - cd2 = [] - cd3 = [] nodeIdentifier = 1 - for n1 in range(len(nodesList)): - x = nodesList[n1][0] - dx_ds1 = nodesList[n1][1] - dx_ds2 = nodesList[n1][2] - dx_ds3 = nodesList[n1][3] - cx.append(x) - cd1.append(dx_ds1) - cd2.append(dx_ds2) - cd3.append(dx_ds3) - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cx[n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, cd1[n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, cd2[n1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, cd3[n1]) - nodeIdentifier += 1 - - print('len(cx)', len(cx)) - print('len(cd1)', len(cd1)) - print('len(cd2)', len(cd2)) + lNodeIds = [] + nix = 0 + for n3 in range(elementsCount3 + 1): + lNodeIds.append([]) + for n2 in range(elementsCount2 + 1): + lNodeIds[n3].append([]) + for n1 in range(elementsCount1 + 1): + lNodeIds[n3][n2].append(None) + # if n3 < elementsCount3: + # if (n1 == 0) and ((n2 == 0) or (n2 == elementsCount2)): + # continue + # else: + # if (n2 == 0) or (n2 == elementsCount2) or (n1 == 0): + # continue + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodesList[nix][0]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, nodesList[nix][1]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, nodesList[nix][2]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, nodesList[nix][3]) + lNodeIds[n3][n2][n1] = nodeIdentifier + nix += 1 + nodeIdentifier += 1 + + # cx = [] + # cd1 = [] + # cd2 = [] + # cd3 = [] + # for n1 in range(len(nodesList)): + # x = nodesList[n1][0] + # dx_ds1 = nodesList[n1][1] + # dx_ds2 = nodesList[n1][2] + # dx_ds3 = nodesList[n1][3] + # cx.append(x) + # cd1.append(dx_ds1) + # cd2.append(dx_ds2) + # cd3.append(dx_ds3) + # node = nodes.createNode(nodeIdentifier, nodetemplate) + # cache.setNode(node) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cx[n1]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, cd1[n1]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, cd2[n1]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, cd3[n1]) + # nodeIdentifier += 1 + # + # print('len(cx)', len(cx)) + # print('len(cd1)', len(cd1)) + # print('len(cd2)', len(cd2)) # Create elements @@ -256,47 +271,133 @@ def generateBaseMesh(cls, region, options): mesh = fm.findMeshByDimension(3) eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives) - eft = eftfactory.createEftBasic() + eftRegular = eftfactory.createEftBasic() + + elementtemplateRegular = mesh.createElementtemplate() + elementtemplateRegular.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplateRegular.defineField(coordinates, -1, eftRegular) - elementtemplate = mesh.createElementtemplate() - elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplate.defineField(coordinates, -1, eft) + elementtemplateCustom = mesh.createElementtemplate() + elementtemplateCustom.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + eft1 = eftfactory.createEftWedgeCollapseXi1AtXi2Zero() + eft2 = eftfactory.createEftWedgeCollapseXi1AtXi2One() + eft3 = eftfactory.createEftWedgeCollapseXi3AtXi2Zero() + # eft4 = eftfactory.createEftWedgeCollapseXi3AtXi2One() + eft5 = eftfactory.createEftWedgeCollapseXi3AtXi1Zero() - # elementtemplate2 = mesh.createElementtemplate() - # elementtemplate2.setElementShapeType(Element.SHAPE_TYPE_CUBE) - # eft2 = eftfactory.createEftWedgeXi1Zero() - # elementtemplate2.defineField(coordinates, -1, eft2) elementIdentifier = 1 - elementsCountAlongMajorDiameter = 4 - elementsCountAlongMinorDiameter = 2 - elementsCountAlong = 2 - # - nel = (elementsCountAlongMinorDiameter + 1) * (elementsCountAlongMajorDiameter + 1) - for e3 in range(elementsCountAlong): - for e2 in range(elementsCountAlongMajorDiameter): - # First row of elements #4 - for e1 in range(elementsCountAlongMinorDiameter): - if (e2 == 0 or e2 == elementsCountAlongMajorDiameter - 1) and e1 == 0: - # nodeIdentifiers = [10, 11, 14, 25, 26, 29] - # element = mesh.createElement(elementIdentifier, elementtemplate2) - # element.setNodesByIdentifier(eft2, nodeIdentifiers) - pass + for e3 in range(elementsCount3): + for e2 in range(elementsCount2): + for e1 in range(elementsCount1): + eft = eftRegular + nodeIdentifiers = [ + lNodeIds[e3 ][e2][e1], lNodeIds[e3 ][e2][e1 + 1], lNodeIds[e3 ][e2 + 1][e1], lNodeIds[e3 ][e2 + 1][e1 + 1], + lNodeIds[e3 + 1][e2][e1], lNodeIds[e3 + 1][e2][e1 + 1], lNodeIds[e3 + 1][e2 + 1][e1], lNodeIds[e3 + 1][e2 + 1][e1 + 1] ] + scalefactors = None + if (e3 < elementsCount3 - 1): + if (e2 == 0) and (e1 == 0): + # Back wedge elements + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(0) + eft = eft1 + scalefactors = [ -1.0 ] + elif (e2 == (elementsCount2 - 1)) and (e1 == 0): + # Front wedge elements + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(2) + eft = eft2 + + if eft is eftRegular: + element = mesh.createElement(elementIdentifier, elementtemplateRegular) + else: + elementtemplateCustom.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplateCustom) + element.setNodesByIdentifier(eft, nodeIdentifiers) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + elementIdentifier += 1 else: - element = mesh.createElement(elementIdentifier, elementtemplate) - bni1 = e3 * nel + e2 * (elementsCountAlongMinorDiameter + 1) + e1 + 1 - bni2 = e3 * nel + e2 * (elementsCountAlongMinorDiameter + 1) + (e1 + 1) % (elementsCountAlongMinorDiameter + 1) + 1 - bni3 = bni1 + elementsCountAlongMinorDiameter + 1 - bni4 = bni3 + 1 - bni5 = bni1 + nel - bni6 = bni5 + 1 - bni7 = bni5 + elementsCountAlongMinorDiameter + 1 - bni8 = bni7 + 1 - nodeIdentifiers = [bni1, bni2, bni3, bni4, - bni5, bni6, bni7, bni8] - result = element.setNodesByIdentifier(eft, nodeIdentifiers) + if ((e2 == 0) or (e2 == elementsCount2 - 1)) and (e1 == 0): + # Top tetrahedron elements + continue + elif (e2 == 0) and (e1 == 1): + # Top back wedge elements + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(5) + eft = eft3 + print('nodeIdentifiers', nodeIdentifiers) + # elif (e2 == 1) and (e1 == 0): + # # Top middle back wedge element + # nodeIdentifiers.pop(4) + # nodeIdentifiers.pop(6) + # eft = eft5 + + if eft is eftRegular: + element = mesh.createElement(elementIdentifier, elementtemplateRegular) + else: + elementtemplateCustom.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplateCustom) + element.setNodesByIdentifier(eft, nodeIdentifiers) + if scalefactors: + element.setScaleFactors(eft, scalefactors) elementIdentifier += 1 + + + # # Top back wedge element + # nodeIdentifiers = [47, 48, 50, 51, 65, 66] + # eft = eft3 + # if eft is eftRegular: + # element = mesh.createElement(elementIdentifier, elementtemplateRegular) + # else: + # elementtemplateCustom.defineField(coordinates, -1, eft) + # element = mesh.createElement(elementIdentifier, elementtemplateCustom) + # element.setNodesByIdentifier(eft, nodeIdentifiers) + # if scalefactors: + # element.setScaleFactors(eft, scalefactors) + # elementIdentifier += 1 + # + # # Top front wedge element + # # nodeIdentifiers = [56, 57, 59, 60, 71, 72] + # # eft = eft4 + # # scalefactors = [-1.0] + # # elementtemplateCustom.defineField(coordinates, -1, eft) + # # element = mesh.createElement(elementIdentifier, elementtemplateCustom) + # # element.setNodesByIdentifier(eft, nodeIdentifiers) + # # # if scalefactors: + # # element.setScaleFactors(eft, scalefactors) + # # elementIdentifier += 1 + # + # # Top middle back wedge element + # nodeIdentifiers = [49, 50, 52, 53, 65, 68] + # eft = eft5 + # if eft is eftRegular: + # element = mesh.createElement(elementIdentifier, elementtemplateRegular) + # else: + # elementtemplateCustom.defineField(coordinates, -1, eft) + # element = mesh.createElement(elementIdentifier, elementtemplateCustom) + # element.setNodesByIdentifier(eft, nodeIdentifiers) + # if scalefactors: + # element.setScaleFactors(eft, scalefactors) + # elementIdentifier += 1 + # + # # Top middle front wedge element + # nodeIdentifiers = [52, 53, 55, 56, 68, 71] + # eft = eft5 + # if eft is eftRegular: + # element = mesh.createElement(elementIdentifier, elementtemplateRegular) + # else: + # elementtemplateCustom.defineField(coordinates, -1, eft) + # element = mesh.createElement(elementIdentifier, elementtemplateCustom) + # element.setNodesByIdentifier(eft, nodeIdentifiers) + # if scalefactors: + # element.setScaleFactors(eft, scalefactors) + # elementIdentifier += 1 + + + annotationGroups = [] fm.endChange() return annotationGroups @@ -314,22 +415,5 @@ def refineMesh(cls, meshrefinement, options): meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong) return - # def createEftWedgeXi1Zero(self): - # ''' - # Create a basic tricubic hermite element template for elements - # along boundary of tenia coli where nodes on xi1 = 0 are collapsed. - # :return: Element field template - # ''' - # eft = self.createEftBasic() - # - # # remap parameters on xi1 = 1 before collapsing nodes - # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS2, []) - # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS1DS2, []) - # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D2_DS1DS2, []) - # remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # - # ln_map = [1, 2, 1, 3, 4, 5, 4, 6] - # remapEftLocalNodes(eft, 6, ln_map) - # assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1Zero: Failed to validate eft' - # return eft + From 4a357054031babc79b446db7eb5899550c0c5ba8 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Mon, 16 Nov 2020 20:03:37 +1300 Subject: [PATCH 04/23] Added different wedge element templates. --- .../utils/eftfactory_tricubichermite.py | 98 +++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 02e26ba0..4f7366ce 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -819,6 +819,104 @@ def createEftWedgeXi1Zero(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1Zero: Failed to validate eft' return eft + def createEftWedgeCollapseXi1AtXi2Zero(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 0. + :return: Element field template + ''' + eft = self.createEftBasic() + setEftScaleFactorIds(eft, [1], []) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 5], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 1, 2, 3, 4, 4, 5, 6] + remapEftLocalNodes(eft, 6, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi2Zero: Failed to validate eft' + return eft + + def createEftWedgeCollapseXi1AtXi2One(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 1. + :return: Element field template + ''' + eft = self.createEftBasic() + + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 2, 3, 3, 4, 5, 6, 6] + remapEftLocalNodes(eft, 6, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi2One: Failed to validate eft' + return eft + + def createEftWedgeCollapseXi3AtXi2Zero(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi2 = 0. + :return: Element field template + ''' + eft = self.createEftBasic() + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS3, []) + remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 2, 3, 4, 1, 2, 5, 6] + remapEftLocalNodes(eft, 6, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi2Zero: Failed to validate eft' + return eft + + # def createEftWedgeCollapseXi3AtXi2One(self): + # ''' + # Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi2 = 1. + # :return: Element field template + # ''' + # eft = self.createEftBasic() + # + # # remap parameters on xi2 = 1 before collapsing nodes + # remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS3, [1])]) + # + # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D_DS3, []) + # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + # + # ln_map = [1, 2, 3, 4, 5, 6, 3, 4] + # remapEftLocalNodes(eft, 6, ln_map) + # assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi2One: Failed to validate eft' + # return eft + + def createEftWedgeCollapseXi3AtXi1Zero(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi1 = 0. + :return: Element field template + ''' + eft = self.createEftBasic() + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D_DS3, []) + remapEftNodeValueLabel(eft, [5, 7], Node.VALUE_LABEL_D_DS1, [(Node.VALUE_LABEL_D_DS3, [])]) + remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 2, 3, 4, 1, 5, 3, 6] + remapEftLocalNodes(eft, 6, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi1Zero: Failed to validate eft' + return eft + def createEftSplitXi1LeftStraight(self): ''' Create an element field template suitable for the inner elements of the From e78f79417c320a71f3c5b8f8205fd9deacc0eda6 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Tue, 17 Nov 2020 17:11:57 +1300 Subject: [PATCH 05/23] Added new wedge and tetrahedron element templates. --- .../utils/eftfactory_tricubichermite.py | 136 +++++++++++++----- 1 file changed, 97 insertions(+), 39 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 4f7366ce..26755900 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -859,62 +859,120 @@ def createEftWedgeCollapseXi1AtXi2One(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi2One: Failed to validate eft' return eft - def createEftWedgeCollapseXi3AtXi2Zero(self): + def createEftWedgeCollapseXi2RightAtXi3One(self): ''' - Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi2 = 0. + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the right wedge elements. :return: Element field template ''' eft = self.createEftBasic() # remap parameters on xi2 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS3, []) - remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS3, [])]) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) + remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 2, 3, 4, 5, 6, 5, 6] + remapEftLocalNodes(eft, 6, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2RightAtXi3One: Failed to validate eft' + return eft + + def createEftWedgeCollapseXi2LeftAtXi3One(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the right wedge elements. + :return: Element field template + ''' + eft = self.createEftBasic() + setEftScaleFactorIds(eft, [1], []) + + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) + # remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - ln_map = [1, 2, 3, 4, 1, 2, 5, 6] + ln_map = [1, 2, 3, 4, 5, 6, 5, 6] remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi2Zero: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2LeftAtXi3One: Failed to validate eft' return eft - # def createEftWedgeCollapseXi3AtXi2One(self): - # ''' - # Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi2 = 1. - # :return: Element field template - # ''' - # eft = self.createEftBasic() - # - # # remap parameters on xi2 = 1 before collapsing nodes - # remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS3, [1])]) - # - # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D_DS3, []) - # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - # remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # - # ln_map = [1, 2, 3, 4, 5, 6, 3, 4] - # remapEftLocalNodes(eft, 6, ln_map) - # assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi2One: Failed to validate eft' - # return eft - - def createEftWedgeCollapseXi3AtXi1Zero(self): - ''' - Create a tricubic hermite element field for a wedge element, where xi3 collapsed on xi1 = 0. + def createEftWedgeCollapseXi1AtXi3One(self): + ''' + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1. :return: Element field template ''' eft = self.createEftBasic() # remap parameters on xi2 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D_DS3, []) - remapEftNodeValueLabel(eft, [5, 7], Node.VALUE_LABEL_D_DS1, [(Node.VALUE_LABEL_D_DS3, [])]) - remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [1, 3, 5, 7], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [5, 7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - ln_map = [1, 2, 3, 4, 1, 5, 3, 6] + ln_map = [1, 2, 3, 4, 5, 5, 6, 6] remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi3AtXi1Zero: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi3One: Failed to validate eft' + return eft + + def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero(self): + ''' + Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 + , and then, xi1 for right nodes on xi3 = 0 is collapsed on xi2 = 0. + :return: Element field template + ''' + eft = self.createEftBasic() + setEftScaleFactorIds(eft, [1], []) + + # remap parameters on xi3 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 2, 5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) + + remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 1, 2, 3, 4, 4, 4, 4] + remapEftLocalNodes(eft, 4, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero: Failed to validate eft' + return eft + + def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero(self): + ''' + Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 + , and then, xi1 for left nodes on xi3 = 0 is collapsed on xi2 = 0. + :return: Element field template + ''' + eft = self.createEftBasic() + setEftScaleFactorIds(eft, [1], []) + + # remap parameters on xi3 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 4, 5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) + + remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, [5], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + + # zero other cross derivative parameters + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + ln_map = [1, 2, 3, 3, 4, 4, 4, 4] + remapEftLocalNodes(eft, 4, ln_map) + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero: Failed to validate eft' return eft def createEftSplitXi1LeftStraight(self): From e13480b2b7059c021b2d51be5a195cfb9387e9d0 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Tue, 17 Nov 2020 17:12:55 +1300 Subject: [PATCH 06/23] Created tetrahedron elements at the top. --- .../meshtypes/meshtype_3d_lungs1.py | 142 ++++++------------ 1 file changed, 50 insertions(+), 92 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py index f604453b..e9a09be0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py @@ -218,6 +218,7 @@ def generateBaseMesh(cls, region, options): nodeIdentifier = 1 lNodeIds = [] nix = 0 + # d1 = [ 1.0, 0.0, 0.0 ] for n3 in range(elementsCount3 + 1): lNodeIds.append([]) for n2 in range(elementsCount2 + 1): @@ -232,6 +233,7 @@ def generateBaseMesh(cls, region, options): # continue node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) + # x = [ n1*1.0, n2*1.0, n3*1.0 ] coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodesList[nix][0]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, nodesList[nix][1]) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, nodesList[nix][2]) @@ -261,10 +263,6 @@ def generateBaseMesh(cls, region, options): # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, cd3[n1]) # nodeIdentifier += 1 # - # print('len(cx)', len(cx)) - # print('len(cd1)', len(cd1)) - # print('len(cd2)', len(cd2)) - # Create elements from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite @@ -282,10 +280,11 @@ def generateBaseMesh(cls, region, options): eft1 = eftfactory.createEftWedgeCollapseXi1AtXi2Zero() eft2 = eftfactory.createEftWedgeCollapseXi1AtXi2One() - eft3 = eftfactory.createEftWedgeCollapseXi3AtXi2Zero() - # eft4 = eftfactory.createEftWedgeCollapseXi3AtXi2One() - eft5 = eftfactory.createEftWedgeCollapseXi3AtXi1Zero() - + eft3 = eftfactory.createEftWedgeCollapseXi2RightAtXi3One() + eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() + eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() + eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero() + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero() elementIdentifier = 1 for e3 in range(elementsCount3): @@ -308,95 +307,54 @@ def generateBaseMesh(cls, region, options): nodeIdentifiers.pop(6) nodeIdentifiers.pop(2) eft = eft2 - - if eft is eftRegular: - element = mesh.createElement(elementIdentifier, elementtemplateRegular) - else: - elementtemplateCustom.defineField(coordinates, -1, eft) - element = mesh.createElement(elementIdentifier, elementtemplateCustom) - element.setNodesByIdentifier(eft, nodeIdentifiers) - if scalefactors: - element.setScaleFactors(eft, scalefactors) - elementIdentifier += 1 else: - if ((e2 == 0) or (e2 == elementsCount2 - 1)) and (e1 == 0): - # Top tetrahedron elements - continue - elif (e2 == 0) and (e1 == 1): + if (e2 == 0) and (e1 == 1): # Top back wedge elements - nodeIdentifiers.pop(4) nodeIdentifiers.pop(5) + nodeIdentifiers.pop(4) eft = eft3 - print('nodeIdentifiers', nodeIdentifiers) - # elif (e2 == 1) and (e1 == 0): - # # Top middle back wedge element - # nodeIdentifiers.pop(4) - # nodeIdentifiers.pop(6) - # eft = eft5 - - if eft is eftRegular: - element = mesh.createElement(elementIdentifier, elementtemplateRegular) - else: - elementtemplateCustom.defineField(coordinates, -1, eft) - element = mesh.createElement(elementIdentifier, elementtemplateCustom) - element.setNodesByIdentifier(eft, nodeIdentifiers) - if scalefactors: - element.setScaleFactors(eft, scalefactors) - elementIdentifier += 1 - - - - # # Top back wedge element - # nodeIdentifiers = [47, 48, 50, 51, 65, 66] - # eft = eft3 - # if eft is eftRegular: - # element = mesh.createElement(elementIdentifier, elementtemplateRegular) - # else: - # elementtemplateCustom.defineField(coordinates, -1, eft) - # element = mesh.createElement(elementIdentifier, elementtemplateCustom) - # element.setNodesByIdentifier(eft, nodeIdentifiers) - # if scalefactors: - # element.setScaleFactors(eft, scalefactors) - # elementIdentifier += 1 - # - # # Top front wedge element - # # nodeIdentifiers = [56, 57, 59, 60, 71, 72] - # # eft = eft4 - # # scalefactors = [-1.0] - # # elementtemplateCustom.defineField(coordinates, -1, eft) - # # element = mesh.createElement(elementIdentifier, elementtemplateCustom) - # # element.setNodesByIdentifier(eft, nodeIdentifiers) - # # # if scalefactors: - # # element.setScaleFactors(eft, scalefactors) - # # elementIdentifier += 1 - # - # # Top middle back wedge element - # nodeIdentifiers = [49, 50, 52, 53, 65, 68] - # eft = eft5 - # if eft is eftRegular: - # element = mesh.createElement(elementIdentifier, elementtemplateRegular) - # else: - # elementtemplateCustom.defineField(coordinates, -1, eft) - # element = mesh.createElement(elementIdentifier, elementtemplateCustom) - # element.setNodesByIdentifier(eft, nodeIdentifiers) - # if scalefactors: - # element.setScaleFactors(eft, scalefactors) - # elementIdentifier += 1 - # - # # Top middle front wedge element - # nodeIdentifiers = [52, 53, 55, 56, 68, 71] - # eft = eft5 - # if eft is eftRegular: - # element = mesh.createElement(elementIdentifier, elementtemplateRegular) - # else: - # elementtemplateCustom.defineField(coordinates, -1, eft) - # element = mesh.createElement(elementIdentifier, elementtemplateCustom) - # element.setNodesByIdentifier(eft, nodeIdentifiers) - # if scalefactors: - # element.setScaleFactors(eft, scalefactors) - # elementIdentifier += 1 - + elif (e2 == elementsCount2 - 1) and (e1 == 1): + # Top front wedge elements + nodeIdentifiers.pop(7) + nodeIdentifiers.pop(6) + eft = eft4 + scalefactors = [ -1.0 ] + elif (e2 == 1) and (e1 == 0): + # Top middle back wedge element + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + eft = eft5 + elif (e2 == 2) and (e1 == 0): + # Top middle front wedge element + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + eft = eft5 + if (e2 == 0) and (e1 == 0): + # Top back tetrahedron elements + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(5) + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(0) + eft = eft6 + scalefactors = [ -1.0 ] + if (e2 == elementsCount2 - 1) and (e1 == 0): + # Top front tetrahedron elements + nodeIdentifiers.pop(7) + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(2) + eft = eft7 + scalefactors = [-1.0] + if eft is eftRegular: + element = mesh.createElement(elementIdentifier, elementtemplateRegular) + else: + elementtemplateCustom.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplateCustom) + element.setNodesByIdentifier(eft, nodeIdentifiers) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + elementIdentifier += 1 annotationGroups = [] fm.endChange() From e808a18e9034fe6a0c057eda8d850476bd3d1956 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 10:43:32 +1300 Subject: [PATCH 07/23] Changed the lung meshtype name. --- src/scaffoldmaker/scaffolds.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index dce08938..26432774 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -29,7 +29,7 @@ from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase1 import MeshType_3d_heartventriclesbase1 from scaffoldmaker.meshtypes.meshtype_3d_heartventriclesbase2 import MeshType_3d_heartventriclesbase2 from scaffoldmaker.meshtypes.meshtype_3d_lens1 import MeshType_3d_lens1 -from scaffoldmaker.meshtypes.meshtype_3d_lungs1 import MeshType_3d_lungs1 +from scaffoldmaker.meshtypes.meshtype_3d_lung1 import MeshType_3d_lung1 from scaffoldmaker.meshtypes.meshtype_3d_ostium1 import MeshType_3d_ostium1 from scaffoldmaker.meshtypes.meshtype_3d_smallintestine1 import MeshType_3d_smallintestine1 from scaffoldmaker.meshtypes.meshtype_3d_solidsphere1 import MeshType_3d_solidsphere1 @@ -73,7 +73,7 @@ def __init__(self): MeshType_3d_heartventriclesbase1, MeshType_3d_heartventriclesbase2, MeshType_3d_lens1, - MeshType_3d_lungs1, + MeshType_3d_lung1, MeshType_3d_ostium1, MeshType_3d_smallintestine1, MeshType_3d_solidsphere1, From 6bb756c82e8452d1cf3b158c8a3606f8fc1f518e Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 11:19:35 +1300 Subject: [PATCH 08/23] Added lung terms. --- src/scaffoldmaker/annotation/lung_terms.py | 27 ++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 src/scaffoldmaker/annotation/lung_terms.py diff --git a/src/scaffoldmaker/annotation/lung_terms.py b/src/scaffoldmaker/annotation/lung_terms.py new file mode 100644 index 00000000..9b8e7df3 --- /dev/null +++ b/src/scaffoldmaker/annotation/lung_terms.py @@ -0,0 +1,27 @@ +""" +Common resource for lungs annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +lung_terms = [ + ("lung", "UBERON:0002048", "ILX:0726937"), + ("left lung", "UBERON:0002168", "ILX:0733450"), + ("right lung", "UBERON:0002167", "ILX:0729582"), + ("upper lobe of left lung", "UBERON:0008952", "ILX:0735339"), + ("lower lobe of left lung", "UBERON:0008953", "ILX:0735534"), + ("upper lobe of right lung", "UBERON:0002170", "ILX:0728821"), + ("middle lobe of right lung", "UBERON:0002174", "ILX:0733737"), + ("lower lobe of right lung", "UBERON:0002171", "ILX:0725712"), + ("right lung accessory lobe", "UBERON:0004890", "ILX:0728696") +] + +def get_lung_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in lung_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Lung annotation term '" + name + "' not found.") From d531c3d88e283f4e56d6a03cfe6e034985d82fe9 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 11:38:40 +1300 Subject: [PATCH 09/23] Added annotation groups and marker point. --- .../meshtypes/meshtype_3d_lung1.py | 232 ++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py new file mode 100644 index 00000000..0256f4ba --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -0,0 +1,232 @@ +''' +Generates 3D lung surface mesh. +''' + +from scaffoldmaker.annotation.annotationgroup import AnnotationGroup +from scaffoldmaker.annotation.lung_terms import get_lung_term +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, \ + findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString +from opencmiss.zinc.element import Element +from opencmiss.zinc.field import Field +from opencmiss.zinc.node import Node + + +class MeshType_3d_lung1(Scaffold_base): + ''' + 3D lung scaffold. + ''' + + @staticmethod + def getName(): + return '3D Lung 1' + + @staticmethod + def getParameterSetNames(): + return [ + 'Default', + 'Mouse 1'] + + @classmethod + def getDefaultOptions(cls, parameterSetName='Default'): + options = {} + options['Base parameter set'] = parameterSetName + + return options + + @staticmethod + def getOrderedOptionNames(): + optionNames = [] + return optionNames + + @classmethod + def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): + ''' + :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. + :return: ScaffoldPackage. + ''' + if parameterSetName: + assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ + 'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName() + assert False, cls.__name__ + '.getOptionScaffoldPackage: Option ' + optionName + ' is not a scaffold' + + @classmethod + def generateBaseMesh(cls, region, options): + ''' + Generate the base tricubic Hermite mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: annotationGroups + ''' + fm = region.getFieldmodule() + fm.beginChange() + coordinates = findOrCreateFieldCoordinates(fm) + + nodes = fm.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) + + mesh = fm.findMeshByDimension(3) + cache = fm.createFieldcache() + + elementsCount1 = 2 + elementsCount2 = 4 + elementsCount3 = 4 + + # Annotation fiducial point + markerGroup = findOrCreateFieldGroup(fm, "marker") + markerName = findOrCreateFieldStoredString(fm, name="marker_name") + markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") + + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup() + markerTemplateInternal = nodes.createNodetemplate() + markerTemplateInternal.defineField(markerName) + markerTemplateInternal.defineField(markerLocation) + + # Create nodes + nodeIdentifier = 1 + lNodeIds = [] + d1 = [0.5, 0.0, 0.0] + d2 = [0.0, 0.5, 0.0] + d3 = [0.0, 0.0, 1.0] + for n3 in range(elementsCount3 + 1): + lNodeIds.append([]) + for n2 in range(elementsCount2 + 1): + lNodeIds[n3].append([]) + for n1 in range(elementsCount1 + 1): + lNodeIds[n3][n2].append([]) + if n3 < elementsCount3: + if (n1 == 0) and ((n2 == 0) or (n2 == elementsCount2)): + continue + else: + if (n2 == 0) or (n2 == elementsCount2) or (n1 == 0): + continue + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + x = [0.5 * (n1 - 1), 0.5 * (n2 - 1), 1.0 * n3] + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, d3) + lNodeIds[n3][n2][n1] = nodeIdentifier + nodeIdentifier += 1 + + # Create elements + mesh = fm.findMeshByDimension(3) + eftfactory = eftfactory_tricubichermite(mesh, None) + eftRegular = eftfactory.createEftBasic() + + elementtemplateRegular = mesh.createElementtemplate() + elementtemplateRegular.setElementShapeType(Element.SHAPE_TYPE_CUBE) + elementtemplateRegular.defineField(coordinates, -1, eftRegular) + + elementtemplateCustom = mesh.createElementtemplate() + elementtemplateCustom.setElementShapeType(Element.SHAPE_TYPE_CUBE) + + lungGroup = AnnotationGroup(region, get_lung_term("lung")) + leftLungGroup = AnnotationGroup(region, get_lung_term("left lung")) + rightLungGroup = AnnotationGroup(region, get_lung_term("right lung")) + annotationGroups = [leftLungGroup, rightLungGroup, lungGroup] + + lungMeshGroup = lungGroup.getMeshGroup(mesh) + leftLungMeshGroup = leftLungGroup.getMeshGroup(mesh) + rightLungMeshGroup = rightLungGroup.getMeshGroup(mesh) + + eft1 = eftfactory.createEftWedgeCollapseXi1AtXi2Zero() + eft2 = eftfactory.createEftWedgeCollapseXi1AtXi2One() + eft3 = eftfactory.createEftWedgeCollapseXi2RightAtXi3One() + eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() + eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() + eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero() + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero() + + elementIdentifier = 1 + for e3 in range(elementsCount3): + for e2 in range(elementsCount2): + for e1 in range(elementsCount1): + eft = eftRegular + nodeIdentifiers = [ + lNodeIds[e3 ][e2][e1], lNodeIds[e3 ][e2][e1 + 1], lNodeIds[e3 ][e2 + 1][e1], lNodeIds[e3 ][e2 + 1][e1 + 1], + lNodeIds[e3 + 1][e2][e1], lNodeIds[e3 + 1][e2][e1 + 1], lNodeIds[e3 + 1][e2 + 1][e1], lNodeIds[e3 + 1][e2 + 1][e1 + 1]] + scalefactors = None + if (e3 < elementsCount3 - 1): + if (e2 == 0) and (e1 == 0): + # Back wedge elements + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(0) + eft = eft1 + scalefactors = [-1.0] + elif (e2 == elementsCount2 - 1) and (e1 == 0): + # Front wedge elements + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(2) + eft = eft2 + else: + if (e2 == 0) and (e1 == 1): + # Top back wedge elements + nodeIdentifiers.pop(5) + nodeIdentifiers.pop(4) + eft = eft3 + elif (e2 == elementsCount2 - 1) and (e1 == 1): + # Top front wedge elements + nodeIdentifiers.pop(7) + nodeIdentifiers.pop(6) + eft = eft4 + scalefactors = [-1.0] + elif (e2 == 1) and (e1 == 0): + # Top middle back wedge element + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + eft = eft5 + elif (e2 == 2) and (e1 == 0): + # Top middle front wedge element + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + eft = eft5 + if (e2 == 0) and (e1 == 0): + # Top back tetrahedron element + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(5) + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(0) + eft = eft6 + scalefactors = [-1.0] + if (e2 == elementsCount2 - 1) and (e1 == 0): + # Top front tetrahedron element + nodeIdentifiers.pop(7) + nodeIdentifiers.pop(6) + nodeIdentifiers.pop(4) + nodeIdentifiers.pop(2) + eft = eft7 + scalefactors = [-1.0] + + if eft is eftRegular: + element = mesh.createElement(elementIdentifier, elementtemplateRegular) + else: + elementtemplateCustom.defineField(coordinates, -1, eft) + element = mesh.createElement(elementIdentifier, elementtemplateCustom) + element.setNodesByIdentifier(eft, nodeIdentifiers) + if scalefactors: + element.setScaleFactors(eft, scalefactors) + elementIdentifier += 1 + leftLungMeshGroup.addElement(element) + lungMeshGroup.addElement(element) + + # Apex annotation point + idx = elementsCount1 * elementsCount2 * (elementsCount3 - 1) + elementsCount1 * (elementsCount2 // 2) + element1 = mesh.findElementByIdentifier(idx) + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, 'APEX') + markerLocation.assignMeshLocation(cache, element1, [1.0, 1.0, 1.0]) + + fm.endChange() + return annotationGroups + From 5a6d4733b3dd23f266a28e8ea8aa201e7f0f475f Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 12:08:18 +1300 Subject: [PATCH 10/23] Changes two functions names. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index 0256f4ba..e5e8b794 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -143,8 +143,8 @@ def generateBaseMesh(cls, region, options): eft3 = eftfactory.createEftWedgeCollapseXi2RightAtXi3One() eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() - eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero() - eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero() + eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() elementIdentifier = 1 for e3 in range(elementsCount3): From 3cb481e42c0aa79b565ae3cd93f59f24365a6d52 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 12:17:39 +1300 Subject: [PATCH 11/23] Added definition for the new element field templates. --- .../utils/eftfactory_tricubichermite.py | 60 ++++++++++++------- 1 file changed, 39 insertions(+), 21 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 26755900..714201a8 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -829,10 +829,12 @@ def createEftWedgeCollapseXi1AtXi2Zero(self): # remap parameters on xi2 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [1, 5], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS1, []) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln_map = [1, 1, 2, 3, 4, 4, 5, 6] @@ -850,8 +852,11 @@ def createEftWedgeCollapseXi1AtXi2One(self): # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln_map = [1, 2, 3, 3, 4, 5, 6, 6] @@ -861,16 +866,20 @@ def createEftWedgeCollapseXi1AtXi2One(self): def createEftWedgeCollapseXi2RightAtXi3One(self): ''' - Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the right wedge elements. + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the + wedge elements. Top right Nodes are collapsed. :return: Element field template ''' eft = self.createEftBasic() - # remap parameters on xi2 = 0 before collapsing nodes + # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln_map = [1, 2, 3, 4, 5, 6, 5, 6] @@ -880,19 +889,21 @@ def createEftWedgeCollapseXi2RightAtXi3One(self): def createEftWedgeCollapseXi2LeftAtXi3One(self): ''' - Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the right wedge elements. + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the + wedge elements. Top left nodes are collapsed. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) - # remap parameters on xi2 = 0 before collapsing nodes + # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - # remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln_map = [1, 2, 3, 4, 5, 6, 5, 6] @@ -902,16 +913,19 @@ def createEftWedgeCollapseXi2LeftAtXi3One(self): def createEftWedgeCollapseXi1AtXi3One(self): ''' - Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1. + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi3 = 1. :return: Element field template ''' eft = self.createEftBasic() - # remap parameters on xi2 = 0 before collapsing nodes + # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [5, 7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) ln_map = [1, 2, 3, 4, 5, 5, 6, 6] @@ -919,24 +933,26 @@ def createEftWedgeCollapseXi1AtXi3One(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi3One: Failed to validate eft' return eft - def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero(self): + def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero(self): ''' Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 - , and then, xi1 for right nodes on xi3 = 0 is collapsed on xi2 = 0. + , and then, xi1 for nodes on both xi2 = 0 and xi3 = 0 collapsed. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [1, 2, 5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - - remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) remapEftNodeValueLabel(eft, [7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - # zero other cross derivative parameters + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) @@ -947,24 +963,26 @@ def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero: Failed to validate eft' return eft - def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero(self): + def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One(self): ''' Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 - , and then, xi1 for left nodes on xi3 = 0 is collapsed on xi2 = 0. + , and then, xi1 for nodes on both xi2 = 1 and xi3 = 0 is collapsed. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [3, 4, 5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - - remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) remapEftNodeValueLabel(eft, [5], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - # zero other cross derivative parameters + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + + # zero cross derivative parameters remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) From 3c16be2ee191d64bfc19ec8cfb56921015c8dcac Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 12:42:35 +1300 Subject: [PATCH 12/23] Minor change. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index e5e8b794..d5a93e78 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -144,7 +144,7 @@ def generateBaseMesh(cls, region, options): eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() - eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One() elementIdentifier = 1 for e3 in range(elementsCount3): From a1fbe9a00bb57de6cc98207c37df76c9b1c2230b Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Wed, 18 Nov 2020 12:43:28 +1300 Subject: [PATCH 13/23] Changed names. --- src/scaffoldmaker/utils/eftfactory_tricubichermite.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 714201a8..5c94d746 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -936,7 +936,7 @@ def createEftWedgeCollapseXi1AtXi3One(self): def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero(self): ''' Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 - , and then, xi1 for nodes on both xi2 = 0 and xi3 = 0 collapsed. + , and then, xi1 for nodes on both xi2 = 0 and xi3 = 0 is collapsed. :return: Element field template ''' eft = self.createEftBasic() @@ -960,7 +960,7 @@ def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero(self): ln_map = [1, 1, 2, 3, 4, 4, 4, 4] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero: Failed to validate eft' return eft def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One(self): @@ -990,7 +990,7 @@ def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One(self): ln_map = [1, 2, 3, 3, 4, 4, 4, 4] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One: Failed to validate eft' return eft def createEftSplitXi1LeftStraight(self): From 5a2c4c60b8c3b1c6cf0d8ff9faf4d573bcd6ed3d Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 14:08:10 +1300 Subject: [PATCH 14/23] delet meshtype_3d_lungs1 --- .../meshtypes/meshtype_3d_lungs1.py | 377 ------------------ 1 file changed, 377 deletions(-) delete mode 100644 src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py deleted file mode 100644 index e9a09be0..00000000 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lungs1.py +++ /dev/null @@ -1,377 +0,0 @@ -''' -Generates 3D lung surface mesh, with variable numbers of elements around, along. -''' - -import copy -import math -from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm -from scaffoldmaker.annotation.bladder_terms import get_bladder_term -from scaffoldmaker.meshtypes.meshtype_1d_path1 import MeshType_1d_path1, extractPathParametersFromRegion -from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.scaffoldpackage import ScaffoldPackage -from scaffoldmaker.utils import interpolation as interp -from scaffoldmaker.utils import matrix -from scaffoldmaker.utils import tubemesh -from scaffoldmaker.utils import vector -from scaffoldmaker.utils.annulusmesh import createAnnulusMesh3d -from scaffoldmaker.utils.geometry import createEllipsePoints, createEllipsoidPoints -from scaffoldmaker.utils.interpolation import smoothCubicHermiteDerivativesLine -from scaffoldmaker.utils.tracksurface import TrackSurface, TrackSurfacePosition, calculate_surface_axes -from scaffoldmaker.utils.zinc_utils import exnodeStringFromNodeValues, mesh_destroy_elements_and_nodes_by_identifiers -from opencmiss.zinc.element import Element, Elementbasis -from opencmiss.zinc.field import Field -from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates -from opencmiss.zinc.node import Node - - -class MeshType_3d_lungs1(Scaffold_base): - ''' - 3D lung scaffold. - ''' - - @staticmethod - def getName(): - return '3D Lungs 1' - - @staticmethod - def getParameterSetNames(): - return [ - 'Default', - 'Mouse 1', - 'Human 1'] - - @classmethod - def getDefaultOptions(cls, parameterSetName='Default'): - options = { - 'Number of elements around': 8, - 'Number of elements along': 4, - 'Number of elements through wall': 1, - 'Major diameter': 30.0, - 'Minor diameter': 15.0, - 'Height': 70, - 'Use cross derivatives': False, - 'Refine': False, - 'Refine number of elements along': 4, - 'Refine number of elements around': 4 - } - if 'Mouse' in parameterSetName: - options['Major diameter'] = 20.0 - options['Minor diameter'] = 10.0 - return options - - @staticmethod - def getOrderedOptionNames(): - optionNames = [ - 'Number of elements around', - 'Number of elements along', - 'Major diameter', - 'Minor diameter', - 'Height', - 'Use cross derivatives', - 'Refine', - 'Refine number of elements around', - 'Refine number of elements along'] - return optionNames - - @classmethod - def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): - ''' - :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. - :return: ScaffoldPackage. - ''' - if parameterSetName: - assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ - 'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName() - assert False, cls.__name__ + '.getOptionScaffoldPackage: Option ' + optionName + ' is not a scaffold' - - @classmethod - def checkOptions(cls, options): - for key in [ - 'Number of elements around', - 'Number of elements along', - 'Refine number of elements along', - 'Refine number of elements around']: - if options[key] < 1: - options[key] = 1 - - @classmethod - def generateBaseMesh(cls, region, options): - ''' - Generate the base tricubic Hermite mesh. See also generateMesh(). - :param region: Zinc region to define model in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - :return: annotationGroups - ''' - - useCrossDerivatives = options['Use cross derivatives'] - - fm = region.getFieldmodule() - fm.beginChange() - coordinates = findOrCreateFieldCoordinates(fm) - - nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - nodetemplateApex = nodes.createNodetemplate() - nodetemplateApex.defineField(coordinates) - nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - nodetemplateApex.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) - if useCrossDerivatives: - nodetemplate = nodes.createNodetemplate() - nodetemplate.defineField(coordinates) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - else: - nodetemplate = nodetemplateApex - - - cache = fm.createFieldcache() - - elementsCount1 = 2 - elementsCount2 = 4 - elementsCount3 = 4 - - # Create nodes - nodesList = [[[-1.0, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.5, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 0.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - - [[-1.0, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.5, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 1.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - - [[-1.0, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.5, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 2.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - - [[-1.0, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.5, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 3.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - - [[-1.0, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, -0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 0.5, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-1.0, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[-0.5, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]], - [[0.0, 1.0, 4.0], [0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 1.0]]] - print('len(nodesList)', len(nodesList)) - - nodeIdentifier = 1 - lNodeIds = [] - nix = 0 - # d1 = [ 1.0, 0.0, 0.0 ] - for n3 in range(elementsCount3 + 1): - lNodeIds.append([]) - for n2 in range(elementsCount2 + 1): - lNodeIds[n3].append([]) - for n1 in range(elementsCount1 + 1): - lNodeIds[n3][n2].append(None) - # if n3 < elementsCount3: - # if (n1 == 0) and ((n2 == 0) or (n2 == elementsCount2)): - # continue - # else: - # if (n2 == 0) or (n2 == elementsCount2) or (n1 == 0): - # continue - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - # x = [ n1*1.0, n2*1.0, n3*1.0 ] - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodesList[nix][0]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, nodesList[nix][1]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, nodesList[nix][2]) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, nodesList[nix][3]) - lNodeIds[n3][n2][n1] = nodeIdentifier - nix += 1 - nodeIdentifier += 1 - - # cx = [] - # cd1 = [] - # cd2 = [] - # cd3 = [] - # for n1 in range(len(nodesList)): - # x = nodesList[n1][0] - # dx_ds1 = nodesList[n1][1] - # dx_ds2 = nodesList[n1][2] - # dx_ds3 = nodesList[n1][3] - # cx.append(x) - # cd1.append(dx_ds1) - # cd2.append(dx_ds2) - # cd3.append(dx_ds3) - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, cx[n1]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, cd1[n1]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, cd2[n1]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, cd3[n1]) - # nodeIdentifier += 1 - # - - # Create elements - from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite - - mesh = fm.findMeshByDimension(3) - eftfactory = eftfactory_tricubichermite(mesh, useCrossDerivatives) - eftRegular = eftfactory.createEftBasic() - - elementtemplateRegular = mesh.createElementtemplate() - elementtemplateRegular.setElementShapeType(Element.SHAPE_TYPE_CUBE) - elementtemplateRegular.defineField(coordinates, -1, eftRegular) - - elementtemplateCustom = mesh.createElementtemplate() - elementtemplateCustom.setElementShapeType(Element.SHAPE_TYPE_CUBE) - - eft1 = eftfactory.createEftWedgeCollapseXi1AtXi2Zero() - eft2 = eftfactory.createEftWedgeCollapseXi1AtXi2One() - eft3 = eftfactory.createEftWedgeCollapseXi2RightAtXi3One() - eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() - eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() - eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1RightAtXi2Zero() - eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1LeftAtXi2Zero() - - elementIdentifier = 1 - for e3 in range(elementsCount3): - for e2 in range(elementsCount2): - for e1 in range(elementsCount1): - eft = eftRegular - nodeIdentifiers = [ - lNodeIds[e3 ][e2][e1], lNodeIds[e3 ][e2][e1 + 1], lNodeIds[e3 ][e2 + 1][e1], lNodeIds[e3 ][e2 + 1][e1 + 1], - lNodeIds[e3 + 1][e2][e1], lNodeIds[e3 + 1][e2][e1 + 1], lNodeIds[e3 + 1][e2 + 1][e1], lNodeIds[e3 + 1][e2 + 1][e1 + 1] ] - scalefactors = None - if (e3 < elementsCount3 - 1): - if (e2 == 0) and (e1 == 0): - # Back wedge elements - nodeIdentifiers.pop(4) - nodeIdentifiers.pop(0) - eft = eft1 - scalefactors = [ -1.0 ] - elif (e2 == (elementsCount2 - 1)) and (e1 == 0): - # Front wedge elements - nodeIdentifiers.pop(6) - nodeIdentifiers.pop(2) - eft = eft2 - else: - if (e2 == 0) and (e1 == 1): - # Top back wedge elements - nodeIdentifiers.pop(5) - nodeIdentifiers.pop(4) - eft = eft3 - elif (e2 == elementsCount2 - 1) and (e1 == 1): - # Top front wedge elements - nodeIdentifiers.pop(7) - nodeIdentifiers.pop(6) - eft = eft4 - scalefactors = [ -1.0 ] - elif (e2 == 1) and (e1 == 0): - # Top middle back wedge element - nodeIdentifiers.pop(6) - nodeIdentifiers.pop(4) - eft = eft5 - elif (e2 == 2) and (e1 == 0): - # Top middle front wedge element - nodeIdentifiers.pop(6) - nodeIdentifiers.pop(4) - eft = eft5 - if (e2 == 0) and (e1 == 0): - # Top back tetrahedron elements - nodeIdentifiers.pop(6) - nodeIdentifiers.pop(5) - nodeIdentifiers.pop(4) - nodeIdentifiers.pop(0) - eft = eft6 - scalefactors = [ -1.0 ] - if (e2 == elementsCount2 - 1) and (e1 == 0): - # Top front tetrahedron elements - nodeIdentifiers.pop(7) - nodeIdentifiers.pop(6) - nodeIdentifiers.pop(4) - nodeIdentifiers.pop(2) - eft = eft7 - scalefactors = [-1.0] - - if eft is eftRegular: - element = mesh.createElement(elementIdentifier, elementtemplateRegular) - else: - elementtemplateCustom.defineField(coordinates, -1, eft) - element = mesh.createElement(elementIdentifier, elementtemplateCustom) - element.setNodesByIdentifier(eft, nodeIdentifiers) - if scalefactors: - element.setScaleFactors(eft, scalefactors) - elementIdentifier += 1 - - annotationGroups = [] - fm.endChange() - 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(). - ''' - refineElementsCountAlong = options['Refine number of elements along'] - refineElementsCountAround = options['Refine number of elements around'] - - meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong) - return - - - From dfbf1902428854eab9164a86e50f4c07710d02e3 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 14:15:36 +1300 Subject: [PATCH 15/23] Deleted "getOptionScaffoldPackage". --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index d5a93e78..2cd60438 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -40,17 +40,6 @@ def getOrderedOptionNames(): optionNames = [] return optionNames - @classmethod - def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): - ''' - :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. - :return: ScaffoldPackage. - ''' - if parameterSetName: - assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ - 'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName() - assert False, cls.__name__ + '.getOptionScaffoldPackage: Option ' + optionName + ' is not a scaffold' - @classmethod def generateBaseMesh(cls, region, options): ''' From 49fb6e989de9b07f645ebe1ff287861e838effde Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 14:18:18 +1300 Subject: [PATCH 16/23] Deleted "beginChange" and "endChange". --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index 2cd60438..b85fc1cb 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -49,7 +49,6 @@ def generateBaseMesh(cls, region, options): :return: annotationGroups ''' fm = region.getFieldmodule() - fm.beginChange() coordinates = findOrCreateFieldCoordinates(fm) nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) @@ -216,6 +215,5 @@ def generateBaseMesh(cls, region, options): markerName.assignString(cache, 'APEX') markerLocation.assignMeshLocation(cache, element1, [1.0, 1.0, 1.0]) - fm.endChange() return annotationGroups From 18d0499a6787130aea05f12b0b307681fbd20240 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 14:22:37 +1300 Subject: [PATCH 17/23] Renamed the apex. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index b85fc1cb..99477492 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -212,7 +212,7 @@ def generateBaseMesh(cls, region, options): markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) nodeIdentifier += 1 cache.setNode(markerPoint) - markerName.assignString(cache, 'APEX') + markerName.assignString(cache, 'apex of left lung') markerLocation.assignMeshLocation(cache, element1, [1.0, 1.0, 1.0]) return annotationGroups From 73d3a1cbe117acdf287aded795d5250eaf5525f7 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 18:46:20 +1300 Subject: [PATCH 18/23] Added new wedge and tetrahedron element templates. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index 99477492..5ee22bf1 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -126,11 +126,11 @@ def generateBaseMesh(cls, region, options): leftLungMeshGroup = leftLungGroup.getMeshGroup(mesh) rightLungMeshGroup = rightLungGroup.getMeshGroup(mesh) - eft1 = eftfactory.createEftWedgeCollapseXi1AtXi2Zero() - eft2 = eftfactory.createEftWedgeCollapseXi1AtXi2One() - eft3 = eftfactory.createEftWedgeCollapseXi2RightAtXi3One() - eft4 = eftfactory.createEftWedgeCollapseXi2LeftAtXi3One() - eft5 = eftfactory.createEftWedgeCollapseXi1AtXi3One() + eft1 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi2Plane([1, 5], Xi2=0) + eft2 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi2Plane([3, 7], Xi2=1) + eft3 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([5, 6], Xi3=1) + eft4 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([7, 8], Xi3=1) + eft5 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi3Plane([5, 7], Xi3=1) eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One() From 7982bde012b2535a1cf02d30db664dce231c4b99 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Thu, 19 Nov 2020 19:04:20 +1300 Subject: [PATCH 19/23] Replaced new wedge element templates. --- .../utils/eftfactory_tricubichermite.py | 166 ++++++++++-------- 1 file changed, 89 insertions(+), 77 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 5c94d746..df34d123 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -819,118 +819,130 @@ def createEftWedgeXi1Zero(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1Zero: Failed to validate eft' return eft - def createEftWedgeCollapseXi1AtXi2Zero(self): + def createEftWedgeCollapseXi1QuadrantAtXi2Plane(self, collapseNodes, Xi2=None): ''' - Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 0. + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 0 or Xi2 = 1. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) - # remap parameters on xi2 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [1, 5], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D_DS1, []) + if Xi2 == 0: + nodes = [1, 2, 5, 6] + # remap parameters on xi2 = 0 before collapsing nodes + if collapseNodes == [1, 5]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + elif collapseNodes == [2, 6]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + elif Xi2 == 1: + nodes = [3, 4, 7, 8] + # remap parameters on xi2 = 1 before collapsing nodes + if collapseNodes == [3, 7]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + elif collapseNodes == [4, 8]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) # zero cross derivative parameters - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [1, 2, 5, 6], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - - ln_map = [1, 1, 2, 3, 4, 4, 5, 6] - remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi2Zero: Failed to validate eft' - return eft - - def createEftWedgeCollapseXi1AtXi2One(self): - ''' - Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 1. - :return: Element field template - ''' - eft = self.createEftBasic() - - # remap parameters on xi2 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [3, 7], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # zero cross derivative parameters - remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [3, 4, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + if Xi2 == 0: + ln_map = [1, 1, 2, 3, 4, 4, 5, 6] + elif Xi2 == 1: + ln_map = [1, 2, 3, 3, 4, 5, 6, 6] - ln_map = [1, 2, 3, 3, 4, 5, 6, 6] remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi2One: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1QuadrantAtXi2Plane: Failed to validate eft' return eft - def createEftWedgeCollapseXi2RightAtXi3One(self): + def createEftWedgeCollapseXi1QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): ''' - Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the - wedge elements. Top right Nodes are collapsed. + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi3 = 0 or xi3 = 1. :return: Element field template ''' eft = self.createEftBasic() + setEftScaleFactorIds(eft, [1], []) - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + if Xi3 == 0: + nodes = [1, 2, 3, 4] + # remap parameters on xi3 = 0 before collapsing nodes + if collapseNodes == [1, 3]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + elif collapseNodes == [2, 4]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + elif Xi3 == 1: + nodes = [5, 6, 7, 8] + # remap parameters on xi3 = 1 before collapsing nodes + if collapseNodes == [5, 7]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + elif collapseNodes == [6, 8]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) # zero cross derivative parameters - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) + + if Xi3 == 0: + ln_map = [1, 1, 2, 2, 3, 4, 5, 6] + elif Xi3 == 1: + ln_map = [1, 2, 3, 4, 5, 5, 6, 6] - ln_map = [1, 2, 3, 4, 5, 6, 5, 6] remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2RightAtXi3One: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1QuadrantAtXi3Plane: Failed to validate eft' return eft - def createEftWedgeCollapseXi2LeftAtXi3One(self): + def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): ''' - Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 1 creating the - wedge elements. Top left nodes are collapsed. + Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 0 or Xi3 = 1. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) + if Xi3 == 0: + nodes = [1, 2, 3, 4] + # remap parameters on xi3 = 0 before collapsing nodes + if collapseNodes == [1, 2]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + elif collapseNodes == [3, 4]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + elif Xi3 == 1: + nodes = [5, 6, 7, 8] + # remap parameters on xi3 = 1 before collapsing nodes + if collapseNodes == [5, 6]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + elif collapseNodes == [7, 8]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) # zero cross derivative parameters - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - - ln_map = [1, 2, 3, 4, 5, 6, 5, 6] - remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2LeftAtXi3One: Failed to validate eft' - return eft - - def createEftWedgeCollapseXi1AtXi3One(self): - ''' - Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi3 = 1. - :return: Element field template - ''' - eft = self.createEftBasic() - - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [5, 7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # zero cross derivative parameters - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + if Xi3 == 0: + ln_map = [1, 2, 1, 2, 3, 4, 5, 6] + elif Xi3 == 1: + ln_map = [1, 2, 3, 4, 5, 6, 5, 6] - ln_map = [1, 2, 3, 4, 5, 5, 6, 6] remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1AtXi3One: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2QuadrantAtXi3Plane: Failed to validate eft' return eft def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero(self): From 337f95ac0d888a115d68a33ca77a6221ec9bfbbc Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 20 Nov 2020 08:19:19 +1300 Subject: [PATCH 20/23] Replaced 1 function to create 4 cases of tetrahedron elements that their peak node is on xi3=1. --- .../utils/eftfactory_tricubichermite.py | 89 +++++++++---------- 1 file changed, 44 insertions(+), 45 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index df34d123..85850a94 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -945,64 +945,63 @@ def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2QuadrantAtXi3Plane: Failed to validate eft' return eft - def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero(self): + def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode): ''' Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 - , and then, xi1 for nodes on both xi2 = 0 and xi3 = 0 is collapsed. + , and then two nodes are collapsed on the other end of xi3. + :param peakNode: is the top node of the tetrahedron in this case that xi3 = 1. + :param sideNode: is the node at the other end of xi3 which has another node collapsed onto it. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) + nodes = [5, 6, 7, 8] # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) - remapEftNodeValueLabel(eft, [7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - - # remap parameters on xi2 = 0 before collapsing nodes - remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + if peakNode == 8: + remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + if sideNode == 2: + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + elif peakNode == 7: + remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) + remapEftNodeValueLabel(eft, [8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + if sideNode == 1: + # remap parameters on xi2 = 0 before collapsing nodes + remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [2], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + elif peakNode == 6: + remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, [5], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + if sideNode == 3: + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + elif peakNode == 5: + remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) + remapEftNodeValueLabel(eft, [6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + if sideNode == 3: + # remap parameters on xi2 = 1 before collapsing nodes + remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, [4], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) # zero cross derivative parameters - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) - - ln_map = [1, 1, 2, 3, 4, 4, 4, 4] - remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero: Failed to validate eft' - return eft - - def createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One(self): - ''' - Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 - , and then, xi1 for nodes on both xi2 = 1 and xi3 = 0 is collapsed. - :return: Element field template - ''' - eft = self.createEftBasic() - setEftScaleFactorIds(eft, [1], []) - - # remap parameters on xi3 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D_DS2, []) - remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) - remapEftNodeValueLabel(eft, [5], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - - # remap parameters on xi2 = 1 before collapsing nodes - remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - # zero cross derivative parameters - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, [5, 6, 7, 8], Node.VALUE_LABEL_D3_DS1DS2DS3, []) + if (peakNode == 8 and sideNode == 2) or (peakNode == 7 and sideNode == 1): + ln_map = [1, 1, 2, 3, 4, 4, 4, 4] + elif (peakNode == 6 and sideNode == 3) or (peakNode == 5 and sideNode == 3): + ln_map = [1, 2, 3, 3, 4, 4, 4, 4] - ln_map = [1, 2, 3, 3, 4, 4, 4, 4] remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One: Failed to validate eft' return eft def createEftSplitXi1LeftStraight(self): From 0e03132f468951467d1a31ebd0e04379918ca83b Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 20 Nov 2020 08:20:10 +1300 Subject: [PATCH 21/23] Converted eft6 and eft7 to the new one. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index 5ee22bf1..2992a8e2 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -131,8 +131,8 @@ def generateBaseMesh(cls, region, options): eft3 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([5, 6], Xi3=1) eft4 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([7, 8], Xi3=1) eft5 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi3Plane([5, 7], Xi3=1) - eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2Zero() - eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2AtXi3OneXi1AtXi2One() + eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(8, 2) + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(6, 3) elementIdentifier = 1 for e3 in range(elementsCount3): From b6e6086a89cccadda411be539003f86fcfc581fe Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 20 Nov 2020 12:21:14 +1300 Subject: [PATCH 22/23] Used the renamed functions. --- src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py index 2992a8e2..c8945587 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_lung1.py @@ -126,13 +126,13 @@ def generateBaseMesh(cls, region, options): leftLungMeshGroup = leftLungGroup.getMeshGroup(mesh) rightLungMeshGroup = rightLungGroup.getMeshGroup(mesh) - eft1 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi2Plane([1, 5], Xi2=0) - eft2 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi2Plane([3, 7], Xi2=1) - eft3 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([5, 6], Xi3=1) - eft4 = eftfactory.createEftWedgeCollapseXi2QuadrantAtXi3Plane([7, 8], Xi3=1) - eft5 = eftfactory.createEftWedgeCollapseXi1QuadrantAtXi3Plane([5, 7], Xi3=1) - eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(8, 2) - eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(6, 3) + eft1 = eftfactory.createEftWedgeCollapseXi1Quadrant([1, 5]) + eft2 = eftfactory.createEftWedgeCollapseXi1Quadrant([3, 7]) + eft3 = eftfactory.createEftWedgeCollapseXi2Quadrant([5, 6]) + eft4 = eftfactory.createEftWedgeCollapseXi2Quadrant([7, 8]) + eft5 = eftfactory.createEftWedgeCollapseXi1Quadrant([5, 7]) + eft6 = eftfactory.createEftTetrahedronCollapseXi1Xi2Quadrant(8, 2) + eft7 = eftfactory.createEftTetrahedronCollapseXi1Xi2Quadrant(6, 3) elementIdentifier = 1 for e3 in range(elementsCount3): From 73995ccb1766a2e55427acbe7b3bb41946005845 Mon Sep 17 00:00:00 2001 From: Zohreh Ekhlasi Date: Fri, 20 Nov 2020 14:56:11 +1300 Subject: [PATCH 23/23] Merged different functions and have only 3 for wedge and tetrahedron elements. --- .../utils/eftfactory_tricubichermite.py | 146 ++++++++++-------- 1 file changed, 78 insertions(+), 68 deletions(-) diff --git a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py index 85850a94..91b0442c 100644 --- a/src/scaffoldmaker/utils/eftfactory_tricubichermite.py +++ b/src/scaffoldmaker/utils/eftfactory_tricubichermite.py @@ -819,15 +819,40 @@ def createEftWedgeXi1Zero(self): assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeXi1Zero: Failed to validate eft' return eft - def createEftWedgeCollapseXi1QuadrantAtXi2Plane(self, collapseNodes, Xi2=None): + def createEftWedgeCollapseXi1Quadrant(self, collapseNodes): ''' - Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi2 = 0 or Xi2 = 1. + Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi3 = 0 or xi3 = 1. :return: Element field template ''' eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) - if Xi2 == 0: + valid = True + if collapseNodes in [[1, 3], [2, 4]]: + nodes = [1, 2, 3, 4] + # remap parameters on xi3 = 0 before collapsing nodes + if collapseNodes == [1, 3]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + elif collapseNodes == [2, 4]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + else: + valid = False + ln_map = [1, 1, 2, 2, 3, 4, 5, 6] + elif collapseNodes in [[5, 7], [6, 8]]: + nodes = [5, 6, 7, 8] + # remap parameters on xi3 = 1 before collapsing nodes + if collapseNodes == [5, 7]: + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) + elif collapseNodes == [6, 8]: + remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) + remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + else: + valid = False + ln_map = [1, 2, 3, 4, 5, 5, 6, 6] + elif collapseNodes in [[1, 5], [2, 6]]: nodes = [1, 2, 5, 6] # remap parameters on xi2 = 0 before collapsing nodes if collapseNodes == [1, 5]: @@ -836,7 +861,10 @@ def createEftWedgeCollapseXi1QuadrantAtXi2Plane(self, collapseNodes, Xi2=None): elif collapseNodes == [2, 6]: remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) - elif Xi2 == 1: + else: + valid = False + ln_map = [1, 1, 2, 3, 4, 4, 5, 6] + elif collapseNodes in [[3, 7], [4, 8]]: nodes = [3, 4, 7, 8] # remap parameters on xi2 = 1 before collapsing nodes if collapseNodes == [3, 7]: @@ -845,48 +873,14 @@ def createEftWedgeCollapseXi1QuadrantAtXi2Plane(self, collapseNodes, Xi2=None): elif collapseNodes == [4, 8]: remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) - - # zero cross derivative parameters - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS3, []) - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - - if Xi2 == 0: - ln_map = [1, 1, 2, 3, 4, 4, 5, 6] - elif Xi2 == 1: + else: + valid = False ln_map = [1, 2, 3, 3, 4, 5, 6, 6] + else: + valid = False - remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1QuadrantAtXi2Plane: Failed to validate eft' - return eft - - def createEftWedgeCollapseXi1QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): - ''' - Create a tricubic hermite element field for a wedge element, where xi1 collapsed on xi3 = 0 or xi3 = 1. - :return: Element field template - ''' - eft = self.createEftBasic() - setEftScaleFactorIds(eft, [1], []) - - if Xi3 == 0: - nodes = [1, 2, 3, 4] - # remap parameters on xi3 = 0 before collapsing nodes - if collapseNodes == [1, 3]: - remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) - elif collapseNodes == [2, 4]: - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - elif Xi3 == 1: - nodes = [5, 6, 7, 8] - # remap parameters on xi3 = 1 before collapsing nodes - if collapseNodes == [5, 7]: - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) - remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) - elif collapseNodes == [6, 8]: - remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) - remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) + if not valid: + assert False, "createEftWedgeCollapseXi1Quadrant. Not implemented for collapse nodes " + str(collapseNodes) # zero cross derivative parameters remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) @@ -894,16 +888,11 @@ def createEftWedgeCollapseXi1QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - if Xi3 == 0: - ln_map = [1, 1, 2, 2, 3, 4, 5, 6] - elif Xi3 == 1: - ln_map = [1, 2, 3, 4, 5, 5, 6, 6] - remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1QuadrantAtXi3Plane: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi1Quadrant: Failed to validate eft' return eft - def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): + def createEftWedgeCollapseXi2Quadrant(self, collapseNodes): ''' Create a tricubic hermite element field for a wedge element, where xi2 collapsed on xi3 = 0 or Xi3 = 1. :return: Element field template @@ -911,7 +900,8 @@ def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): eft = self.createEftBasic() setEftScaleFactorIds(eft, [1], []) - if Xi3 == 0: + valid = True + if collapseNodes in [[1, 2], [3, 4]]: nodes = [1, 2, 3, 4] # remap parameters on xi3 = 0 before collapsing nodes if collapseNodes == [1, 2]: @@ -920,7 +910,10 @@ def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): elif collapseNodes == [3, 4]: remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) - elif Xi3 == 1: + else: + valid = False + ln_map = [1, 2, 1, 2, 3, 4, 5, 6] + elif collapseNodes in [[5, 6], [7, 8]]: nodes = [5, 6, 7, 8] # remap parameters on xi3 = 1 before collapsing nodes if collapseNodes == [5, 6]: @@ -929,6 +922,14 @@ def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): elif collapseNodes == [7, 8]: remapEftNodeValueLabel(eft, collapseNodes, Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + else: + valid = False + ln_map = [1, 2, 3, 4, 5, 6, 5, 6] + else: + valid = False + + if not valid: + assert False, "createEftWedgeCollapseXi2Quadrant. Not implemented for collapse nodes " + str(collapseNodes) # zero cross derivative parameters remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) @@ -936,20 +937,15 @@ def createEftWedgeCollapseXi2QuadrantAtXi3Plane(self, collapseNodes, Xi3=None): remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - if Xi3 == 0: - ln_map = [1, 2, 1, 2, 3, 4, 5, 6] - elif Xi3 == 1: - ln_map = [1, 2, 3, 4, 5, 6, 5, 6] - remapEftLocalNodes(eft, 6, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2QuadrantAtXi3Plane: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftWedgeCollapseXi2Quadrant: Failed to validate eft' return eft - def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode): + def createEftTetrahedronCollapseXi1Xi2Quadrant(self, peakNode, sideNode): ''' - Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 = 1 + Create a tricubic hermite element field for a tetrahedron element, where xi1 and xi2 are collapsed on xi3 , and then two nodes are collapsed on the other end of xi3. - :param peakNode: is the top node of the tetrahedron in this case that xi3 = 1. + :param peakNode: is the top node of the tetrahedron. :param sideNode: is the node at the other end of xi3 which has another node collapsed onto it. :return: Element field template ''' @@ -960,6 +956,8 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) # remap parameters on xi3 = 1 before collapsing nodes remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D_DS2, []) + + valid = True if peakNode == 8: remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) remapEftNodeValueLabel(eft, [7], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) @@ -967,6 +965,9 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) # remap parameters on xi2 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [1], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + else: + valid = False + ln_map = [1, 1, 2, 3, 4, 4, 4, 4] elif peakNode == 7: remapEftNodeValueLabel(eft, [5, 6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [])]) remapEftNodeValueLabel(eft, [8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) @@ -974,6 +975,9 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) # remap parameters on xi2 = 0 before collapsing nodes remapEftNodeValueLabel(eft, [1, 2], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [2], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + else: + valid = False + ln_map = [1, 1, 2, 3, 4, 4, 4, 4] elif peakNode == 6: remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) remapEftNodeValueLabel(eft, [5], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [])]) @@ -981,6 +985,9 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [3], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [])]) + else: + valid = False + ln_map = [1, 2, 3, 3, 4, 4, 4, 4] elif peakNode == 5: remapEftNodeValueLabel(eft, [7, 8], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS2, [1])]) remapEftNodeValueLabel(eft, [6], Node.VALUE_LABEL_D_DS3, [(Node.VALUE_LABEL_D_DS1, [1])]) @@ -988,6 +995,14 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) # remap parameters on xi2 = 1 before collapsing nodes remapEftNodeValueLabel(eft, [3, 4], Node.VALUE_LABEL_D_DS1, []) remapEftNodeValueLabel(eft, [4], Node.VALUE_LABEL_D_DS2, [(Node.VALUE_LABEL_D_DS1, [1])]) + else: + valid = False + ln_map = [1, 2, 3, 3, 4, 4, 4, 4] + else: + valid = False + + if not valid: + assert False, "createEftTetrahedronCollapseXi1Xi2Quadrant. Not implemented for peak node " + str(peakNode) + " and side node " + str(sideNode) # zero cross derivative parameters remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS1DS2, []) @@ -995,13 +1010,8 @@ def createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One(self, peakNode, sideNode) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D2_DS2DS3, []) remapEftNodeValueLabel(eft, nodes, Node.VALUE_LABEL_D3_DS1DS2DS3, []) - if (peakNode == 8 and sideNode == 2) or (peakNode == 7 and sideNode == 1): - ln_map = [1, 1, 2, 3, 4, 4, 4, 4] - elif (peakNode == 6 and sideNode == 3) or (peakNode == 5 and sideNode == 3): - ln_map = [1, 2, 3, 3, 4, 4, 4, 4] - remapEftLocalNodes(eft, 4, ln_map) - assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2QuadrantAtXi3One: Failed to validate eft' + assert eft.validate(), 'eftfactory_tricubichermite.createEftTetrahedronCollapseXi1Xi2Quadrant: Failed to validate eft' return eft def createEftSplitXi1LeftStraight(self):