diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_solidcylinder1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_solidcylinder1.py index ed7c9422..0c2fe4e8 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_solidcylinder1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_solidcylinder1.py @@ -1,6 +1,6 @@ """ Generates a solid cylinder using a ShieldMesh of all cube elements, - with variable numbers of elements in major, minor and length directions. + with variable numbers of elements in major, minor, shell and axial directions. """ from __future__ import division @@ -19,7 +19,7 @@ class MeshType_3d_solidcylinder1(Scaffold_base): """ Generates a solid cylinder using a ShieldMesh of all cube elements, -with variable numbers of elements in major, minor and length directions. +with variable numbers of elements in major, minor, shell and axial directions. """ centralPathDefaultScaffoldPackages = { 'Cylinder 1': ScaffoldPackage(MeshType_1d_path1, { @@ -52,6 +52,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Central path': copy.deepcopy(centralPathOption), 'Number of elements across major': 4, 'Number of elements across minor': 4, + 'Number of elements across shell': 0, 'Number of elements along': 1, 'Lower half': False, 'Use cross derivatives': False, @@ -67,6 +68,7 @@ def getOrderedOptionNames(): 'Central path', 'Number of elements across major', 'Number of elements across minor', + 'Number of elements across shell', 'Number of elements along', 'Lower half', 'Refine', @@ -122,6 +124,10 @@ def checkOptions(cls, options): options['Number of elements across minor'] += 1 if options['Number of elements along'] < 1: options['Number of elements along'] = 1 + Rcrit = min(options['Number of elements across major']-4, options['Number of elements across minor']-4)//2 + if options['Number of elements across shell'] > Rcrit: + dependentChanges = True + options['Number of elements across shell'] = Rcrit return dependentChanges @@ -140,6 +146,7 @@ def generateBaseMesh(region, options): if not full: elementsCountAcrossMajor //= 2 elementsCountAcrossMinor = options['Number of elements across minor'] + elementsCountAcrossShell = options['Number of elements across shell'] elementsCountAlong = options['Number of elements along'] useCrossDerivatives = options['Use cross derivatives'] @@ -150,7 +157,7 @@ def generateBaseMesh(region, options): cylinderShape = CylinderShape.CYLINDER_SHAPE_FULL if full else CylinderShape.CYLINDER_SHAPE_LOWER_HALF - base = CylinderEnds(elementsCountAcrossMajor, elementsCountAcrossMinor, [0.0, 0.0, 0.0], + base = CylinderEnds(elementsCountAcrossMajor, elementsCountAcrossMinor, elementsCountAcrossShell, [0.0, 0.0, 0.0], cylinderCentralPath.alongAxis[0], cylinderCentralPath.majorAxis[0], cylinderCentralPath.minorRadii[0]) cylinder1 = CylinderMesh(fm, coordinates, elementsCountAlong, base, diff --git a/src/scaffoldmaker/utils/cylindermesh.py b/src/scaffoldmaker/utils/cylindermesh.py index c549a014..580b2268 100644 --- a/src/scaffoldmaker/utils/cylindermesh.py +++ b/src/scaffoldmaker/utils/cylindermesh.py @@ -42,7 +42,7 @@ class CylinderEnds: Stores base ellipse parameters. """ - def __init__(self, elementsCountAcrossMajor, elementsCountAcrossMinor, + def __init__(self, elementsCountAcrossMajor, elementsCountAcrossMinor, elementsCountAcrossShell=0, centre=None, alongAxis=None, majorAxis=None, minorRadius=None): """ :param elementsCountAcrossMajor: Number of elements across major axis. Must be at least 2 + elementsCountRim for @@ -61,6 +61,7 @@ def __init__(self, elementsCountAcrossMajor, elementsCountAcrossMinor, self._minorAxis = vector.setMagnitude(vector.crossproduct3(alongAxis, majorAxis), minorRadius) self._elementsCountAcrossMinor = elementsCountAcrossMinor self._elementsCountAcrossMajor = elementsCountAcrossMajor + self._elementsCountAcrossShell = elementsCountAcrossShell self._majorRadius = vector.magnitude(majorAxis) self.px = None self.pd1 = None @@ -150,6 +151,7 @@ def __init__(self, fieldModule, coordinates, elementsCountAlong, base=None, end= self._elementsCountAcrossMajor = base._elementsCountAcrossMajor self._elementsCountUp = base._elementsCountAcrossMajor // 2 \ if cylinderShape == CylinderShape.CYLINDER_SHAPE_FULL else base._elementsCountAcrossMajor + self._elementsCountAcrossShell = base._elementsCountAcrossShell self._elementsCountAlong = elementsCountAlong self._elementsCountAround = 2 * (self._elementsCountUp - 2) + self._elementsCountAcrossMinor self._startNodeIdentifier = 1 @@ -165,7 +167,8 @@ def __init__(self, fieldModule, coordinates, elementsCountAlong, base=None, end= self._cylinderCentralPath = cylinderCentralPath if cylinderCentralPath: self.calculateEllipseParams(cylinderCentralPath=self._cylinderCentralPath) - self._base = CylinderEnds(base._elementsCountAcrossMajor, base._elementsCountAcrossMinor, self._centres[0], + self._base = CylinderEnds(base._elementsCountAcrossMajor, base._elementsCountAcrossMinor, + base._elementsCountAcrossShell, self._centres[0], None, self._majorAxis[0], self._minorRadii[0]) else: self._length = vector.magnitude(base._alongAxis) @@ -195,7 +198,7 @@ def createCylinderMesh3d(self, fieldModule, coordinates): nodes = fieldModule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) mesh = fieldModule.findMeshByDimension(3) - elementsCountRim = 0 + elementsCountRim = self._elementsCountAcrossShell shieldMode = ShieldShape.SHIELD_SHAPE_FULL if self._cylinderShape is self._cylinderShape.CYLINDER_SHAPE_FULL \ else ShieldShape.SHIELD_SHAPE_LOWER_HALF @@ -210,7 +213,7 @@ def createCylinderMesh3d(self, fieldModule, coordinates): self._ellipses = [] for n3 in range(n3Count + 1): ellipse = Ellipse2D(self._centres[n3], self._majorAxis[n3], self._minorAxis[n3], - self._elementsCountAcrossMajor, self._elementsCountAcrossMinor, + self._elementsCountAcrossMajor, self._elementsCountAcrossMinor, self._elementsCountAcrossShell, ellipseShape=ellipseShape) self._ellipses.append(ellipse) self.copyEllipsesNodesToShieldNodes(n3) @@ -240,6 +243,7 @@ def createCylinderMesh3d(self, fieldModule, coordinates): if self._end is None: self._end = CylinderEnds(self._elementsCountAcrossMajor, self._elementsCountAcrossMinor, + self._elementsCountAcrossShell, self._centres[-1], self._shield.pd2[-1][0][1], vector.setMagnitude(self._base._majorAxis, self._majorRadii[-1]), self._minorRadii[-1]) @@ -372,7 +376,7 @@ class Ellipse2D: Generate a 2D ellipse. """ - def __init__(self, centre, majorAxis, minorAxis, elementsCountAcrossMajor, elementsCountAcrossMinor, + def __init__(self, centre, majorAxis, minorAxis, elementsCountAcrossMajor, elementsCountAcrossMinor, elementsCountAcrossShell, ellipseShape=EllipseShape.Ellipse_SHAPE_FULL): """ :param centre: Ellipse centre. @@ -387,7 +391,7 @@ def __init__(self, centre, majorAxis, minorAxis, elementsCountAcrossMajor, eleme self.minorAxis = minorAxis self.majorRadius = vector.magnitude(majorAxis) self.minorRadius = vector.magnitude(minorAxis) - elementsCountRim = 0 + elementsCountRim = elementsCountAcrossShell shieldMode = ShieldShape.SHIELD_SHAPE_FULL if ellipseShape is EllipseShape.Ellipse_SHAPE_FULL\ else ShieldShape.SHIELD_SHAPE_LOWER_HALF shield = ShieldMesh(elementsCountAcrossMinor, elementsCountAcrossMajor, elementsCountRim, @@ -397,6 +401,7 @@ def __init__(self, centre, majorAxis, minorAxis, elementsCountAcrossMajor, eleme self.elementsCountAround = shield.elementsCountAroundFull self.elementsCountUp = shield.elementsCountUp self.elementsCountAcrossMinor = elementsCountAcrossMinor + self.elementsCountAcrossShell = elementsCountAcrossShell self.nodeId = shield.nodeId self.px = shield.px[0] self.pd1 = shield.pd1[0] @@ -416,13 +421,9 @@ def generate2DEllipseMesh(self): self.createRegularRowCurves(rscx, rscd1, rscd3) self.createRegularColumnCurves() self.__shield.getTriplePoints(0) - n1b = 1 - m1a = self.elementsCountAcrossMinor - m1b = m1a - 1 - m1c = m1a - 2 - n2b = 1 - self.smoothTriplePointsCurves(n2b, n1b, m1a) + self.smoothTriplePointsCurves() self.__shield.smoothDerivativesToTriplePoints(0, fixAllDirections=True) + self.smoothDerivativesAroundShell() if self.ellipseShape == EllipseShape.Ellipse_SHAPE_FULL: self.generateNodesForUpperHalf() @@ -452,14 +453,14 @@ def setRimNodes(self, nx, nd1, nd2, nd3): btd2 = self.pd2 btd3 = self.pd3 - elementsCountRim = 0 + elementsCountRim = self.elementsCountAcrossShell for n in range(self.elementsCountAround + 1): n1, n2 = self.__shield.convertRimIndex(n) btx[n2][n1] = nx[n] if n2 > elementsCountRim: # regular rows btd1[n2][n1] = nd1[n] btd3[n2][n1] = nd3[n] - if n2 >= 2: + if n2 >= 2 + elementsCountRim: btd3[n2][n1] = vector.setMagnitude(self.minorAxis, vector.dotproduct(nd3[n], self.minorAxis)) else: # around rim btd1[n2][n1] = nd1[n] @@ -513,19 +514,31 @@ def createRegularRowCurves(self, rscx, rscd1, rscd3): btd2 = self.pd2 btd3 = self.pd3 - elementsCountRim = 0 - for n2 in range(elementsCountRim + 2, self.elementsCountUp + 1): + elementsCountRim = self.elementsCountAcrossShell + n2c = elementsCountRim + 2 + n2m = self.elementsCountUp + n1a = elementsCountRim + for n2 in range(n2c, n2m + 1): btx[n2], btd3[n2], pe, pxi, psf = sampleCubicHermiteCurves( [btx[n2][0], rscx[n2], btx[n2][-1]], [vector.setMagnitude(btd3[n2][0], -1.0), rscd3[n2], btd3[n2][-1]], self.elementsCountAcrossMinor, lengthFractionStart=1, lengthFractionEnd=1, arcLengthDerivatives=True) btd1[n2] = interpolateSampleCubicHermite([[-btd1[n2][0][c] for c in range(3)], rscd1[n2], btd1[n2][-1]], [[0.0, 0.0, 0.0]] * 3, pe, pxi, psf)[0] - if n2 == self.elementsCountUp: + + if n2 < n2m: + for n1 in range(n1a+1): + btd1[n2][n1] = [-btd1[n2][n1][c] for c in range(3)] + btd3[n2][n1] = [-btd3[n2][n1][c] for c in range(3)] + elif n2 == n2m: + btd1[n2][0] = [-btd1[n2][0][c] for c in range(3)] + btd3[n2][0] = [-btd3[n2][0][c] for c in range(3)] for n1 in range(1, self.elementsCountAcrossMinor): - btd1[n2][n1] = vector.setMagnitude(btd1[self.elementsCountUp][-1], 1.0) - btd3[n2][0] = [-btd3[n2][0][c] for c in range(3)] - btd1[n2][0] = [-btd1[n2][0][c] for c in range(3)] + if n1 > n1a: + btd1[n2][n1] = vector.setMagnitude(btd1[n2m][-1], 1.0) + else: + btd1[n2][n1] = vector.setMagnitude(btd1[n2m][-1], -1.0) + btd3[n2][n1] = [-btd3[n2][n1][c] for c in range(3)] def createRegularColumnCurves(self): """ @@ -536,30 +549,32 @@ def createRegularColumnCurves(self): btd2 = self.pd2 btd3 = self.pd3 - for n1 in range(2, self.elementsCountAcrossMinor - 1): + elementsCountRim = self.elementsCountAcrossShell + n1c = elementsCountRim + 2 + n1y = self.elementsCountAcrossMinor - elementsCountRim - 1 + n2a = elementsCountRim + n2c = 2 + n2a + n2m = self.elementsCountUp + for n1 in range(n1c, n1y): tx, td1, pe, pxi, psf = sampleCubicHermiteCurves( - [btx[0][n1], btx[2][n1]], [[-btd3[0][n1][c] for c in range(3)], btd1[2][n1]], 2, + [btx[0][n1], btx[n2c][n1]], [[-btd3[0][n1][c] for c in range(3)], btd1[n2c][n1]], 2+elementsCountRim, lengthFractionStart=1, arcLengthDerivatives=True) - for n2 in range(3, self.elementsCountUp + 1): + for n2 in range(n2c + 1, n2m + 1): tx.append(btx[n2][n1]) td1.append(btd1[n2][n1]) td1 = smoothCubicHermiteDerivativesLine(tx, td1, fixStartDirection=True, fixEndDirection=True) - td3 = \ - interpolateSampleCubicHermite([btd1[0][n1], btd3[2][n1]], [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[ - 0] - for n2 in range(self.elementsCountUp + 1): - if n2 < 2: + td3 = interpolateSampleCubicHermite([btd1[0][n1], btd3[n2c][n1]], [[0.0, 0.0, 0.0]] * 2, pe, pxi, psf)[0] + for n2 in range(n2m + 1): + btd1[n2][n1] = td1[n2] + if n2 < n2c: btx[n2][n1] = tx[n2] - if n2 == 0: - btd3[n2][n1] = [-td1[0][c] for c in range(3)] + if n2 <= n2a: + btd1[n2][n1] = td3[n2] + btd3[n2][n1] = [-td1[n2][c] for c in range(3)] else: btd3[n2][n1] = td3[n2] - if n2 == 0: - btd1[n2][n1] = td3[n2] - else: - btd1[n2][n1] = td1[n2] - def smoothTriplePointsCurves(self, n2b, n1b, m1a): + def smoothTriplePointsCurves(self): """ Smooth row and column curves passing triple points (i.e., row 1 and columns 1 and -2). """ @@ -568,19 +583,30 @@ def smoothTriplePointsCurves(self, n2b, n1b, m1a): btd2 = self.pd2 btd3 = self.pd3 - # smooth shield row 1 - btd3[n2b][n1b:m1a] = smoothCubicHermiteDerivativesLine(btx[n2b][n1b:m1a], btd3[n2b][n1b:m1a]) + n1b = 1 + self.elementsCountAcrossShell + n1z = self.elementsCountAcrossMinor - self.elementsCountAcrossShell + n1y = n1z - 1 + n2b = 1 + self.elementsCountAcrossShell + n2m = self.elementsCountUp + + # smooth shield row n2b + btd3[n2b][n1b:n1z] = smoothCubicHermiteDerivativesLine(btx[n2b][n1b:n1z], btd3[n2b][n1b:n1z]) - # smooth Shield columns 1, -2 - for n1 in [1, -2]: + # smooth Shield columns n1b, n1y + for n1 in [n1b, n1y]: tx = [] td1 = [] - for n2 in range(1, self.elementsCountUp + 1): + for n2 in range(n2b, n2m + 1): tx.append(btx[n2][n1]) td1.append(btd1[n2][n1]) - td1 = smoothCubicHermiteDerivativesLine(tx, td1, fixEndDirection=True) - for n in range(self.elementsCountUp): - btd1[n + 1][n1] = td1[n] + td1 = smoothCubicHermiteDerivativesLine(tx, td1, fixEndDirection=True, fixStartDerivative=True) + for n in range(n2m-n2b+1): + btd1[n + n2b][n1] = td1[n] + + def smoothDerivativesAroundShell(self): + """ Smooth curves around shell layers""" + for rx in range(1, self.elementsCountAcrossShell+1): + self.__shield.smoothDerivativesAroundRim(0, n3d=None, rx=rx) def generateNodesForUpperHalf(self): """ diff --git a/src/scaffoldmaker/utils/shieldmesh.py b/src/scaffoldmaker/utils/shieldmesh.py index fc916ac8..0d9e4acc 100644 --- a/src/scaffoldmaker/utils/shieldmesh.py +++ b/src/scaffoldmaker/utils/shieldmesh.py @@ -13,7 +13,8 @@ from scaffoldmaker.utils.mirror import Mirror from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds -from scaffoldmaker.utils.interpolation import DerivativeScalingMode, sampleCubicHermiteCurves, smoothCubicHermiteDerivativesLine +from scaffoldmaker.utils.interpolation import DerivativeScalingMode, sampleCubicHermiteCurves, \ + smoothCubicHermiteDerivativesLine, interpolateSampleCubicHermite from scaffoldmaker.utils.tracksurface import TrackSurface, TrackSurfacePosition, calculate_surface_axes from enum import Enum @@ -31,7 +32,7 @@ class ShieldMesh: ''' def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, trackSurface : TrackSurface=None, - elementsCountAlong=1, shieldMode = ShieldShape.SHIELD_SHAPE_LOWER_HALF, shieldType = ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR): + elementsCountAlong=1, shieldMode=ShieldShape.SHIELD_SHAPE_LOWER_HALF, shieldType=ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR): ''' Data structure for defining a shield-shaped mesh which is flat on the top and rounded around the bottom and/or the same mirror mirrored on top. @@ -53,8 +54,8 @@ def __init__(self, elementsCountAcross, elementsCountUpFull, elementsCountRim, t triple points, and across the shorter bottom row. Extra elements up add regular rows of nodes/elements on top, and extra non-rim elements across add regular columns of nodes/elements up the centre. - :param elementsCountAcross: Number of elements across top of shield. Must be at least 4 + elementsCountRim. - :param elementsCountUpFull: Number of elements up central axis of shield. Must be at least 2 + elementsCountRim if half and 4 + elementsCountRim if full. + :param elementsCountAcross: Number of elements across top of shield. Must be at least 4 + 2*elementsCountRim. + :param elementsCountUpFull: Number of elements up central axis of shield. Must be at least 2 + elementsCountRim if half and 4 + 2*elementsCountRim if full. :param elementsCountAlong: Number of elements through wall for ventricle case (only 1 element is supported) and along cylinder axis in cylinder case. :param elementsCountRim: Number of elements around bottom rim (not top) outside of 'triple points'. :param trackSurface: Optional trackSurface to store or restrict points to. @@ -129,8 +130,8 @@ def getTriplePoints(self, n3): tx, td1 = sampleCubicHermiteCurves( [self.px[n3][n2a][n1c], self.px[n3][n2c][n1b]],[[(-self.pd1[n3][n2a][n1c][c] - self.pd3[n3][n2a][n1c][c]) for c in range(3)], self.pd1[n3][n2c][n1b]],2, arcLengthDerivatives=True)[0:2] ltx.append(tx[1]) - tx, td1 = sampleCubicHermiteCurves( - [ self.px[n3][n2a][n1b], self.px[n3][n2c][n1c] ], [ [-self.pd3[n3][n2a][n1b][c] for c in range(3)], [ (self.pd1[n3][n2c][n1c][c] + self.pd3[n3][n2c][n1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] + # tx, td1 = sampleCubicHermiteCurves( + # [ self.px[n3][n2a][n1b], self.px[n3][n2c][n1c] ], [ [-self.pd3[n3][n2a][n1b][c] for c in range(3)], [ (self.pd1[n3][n2c][n1c][c] + self.pd3[n3][n2c][n1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] ltx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( [ self.px[n3][n2c][n1a], self.px[n3][n2b][n1c] ], [ [ (self.pd1[n3][n2c][n1a][c] - self.pd3[n3][n2c][n1a][c]) for c in range(3) ], self.pd3[n3][n2b][n1c] ], 2, arcLengthDerivatives = True)[0:2] @@ -171,8 +172,8 @@ def getTriplePoints(self, n3): tx, td1 = sampleCubicHermiteCurves( [ self.px[n3][n2a][m1c], self.px[n3][n2c][m1b] ], [ [ (self.pd1[n3][n2a][m1c][c] - self.pd3[n3][n2a][m1c][c]) for c in range(3) ], self.pd1[n3][n2c][m1b] ], 2, arcLengthDerivatives = True)[0:2] rtx.append(tx[1]) - tx, td1 = sampleCubicHermiteCurves( - [ self.px[n3][n2a][m1b], self.px[n3][n2c][m1c] ], [ [-self.pd3[n3][n2a][m1b][c] for c in range(3)], [ (-self.pd3[n3][n2c][m1c][c] + self.pd1[n3][n2c][m1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] + # tx, td1 = sampleCubicHermiteCurves( + # [ self.px[n3][n2a][m1b], self.px[n3][n2c][m1c] ], [ [-self.pd3[n3][n2a][m1b][c] for c in range(3)], [ (-self.pd3[n3][n2c][m1c][c] + self.pd1[n3][n2c][m1c][c]) for c in range(3) ] ], 2, arcLengthDerivatives = True)[0:2] rtx.append(tx[1]) tx, td1 = sampleCubicHermiteCurves( [ self.px[n3][n2c][m1a], self.px[n3][n2b][m1c] ], [ [ (-self.pd1[n3][n2c][m1a][c] - self.pd3[n3][n2c][m1a][c]) for c in range(3) ], [ -d for d in self.pd3[n3][n2b][m1c] ] ], 2, arcLengthDerivatives = True)[0:2] @@ -216,6 +217,8 @@ def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): ''' n1a = self.elementsCountRim n1b = n1a + 1 + n1z = self.elementsCountAcross - self.elementsCountRim + n1y = n1z - 1 m1a = self.elementsCountAcross - self.elementsCountRim m1b = m1a - 1 n2a = self.elementsCountRim @@ -223,23 +226,33 @@ def smoothDerivativesToTriplePoints(self, n3, fixAllDirections=False): n2c = n2a + 2 if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: # left - tx = [] - td3 = [] - for n2 in range(0, n2c): - tx .append(self.px [n3][n2][n1b]) - td3.append([-d for d in self.pd3[n3][n2][n1b]] if (n2 < n2b) else [ (self.pd1[n3][n2][n1b][c] + self.pd3[n3][n2][n1b][c]) for c in range(3) ]) - td3 = smoothCubicHermiteDerivativesLine(tx, td3, fixAllDirections=fixAllDirections, fixEndDerivative=True, magnitudeScalingMode = DerivativeScalingMode.HARMONIC_MEAN) + tx, td3, pe, pxi, psf = sampleCubicHermiteCurves([self.px[n3][0][n1b], self.px[n3][n2b][n1b]], + [[-d for d in self.pd3[n3][0][n1b]], + [self.pd1[n3][n2b][n1b][c]+self.pd3[n3][n2b][n1b][c] for c in range(3)]], + self.elementsCountRim+1, lengthFractionStart=1, + arcLengthDerivatives=True) + td1 = interpolateSampleCubicHermite([self.pd1[n3][0][n1b], self.pd3[n3][n2b][n1b]], [[0.0, 0.0, 0.0]] * 2, + pe, pxi, psf)[0] for n2 in range(0, n2b): + if n2 > 0: + self.px[n3][n2][n1b] = tx[n2] + self.pd1[n3][n2][n1b] = td1[n2] self.pd3[n3][n2][n1b] = [-d for d in td3[n2]] + # right - tx = [] - td3 = [] - for n2 in range(0, n2c): - tx .append(self.px [n3][n2][m1b]) - td3.append([-d for d in self.pd3[n3][n2][m1b]] if (n2 < n2b) else [ (-self.pd3[n3][n2][m1b][c] + self.pd1[n3][n2][m1b][c]) for c in range(3) ]) - td3 = smoothCubicHermiteDerivativesLine(tx, td3, fixAllDirections=fixAllDirections, fixEndDerivative=True, magnitudeScalingMode = DerivativeScalingMode.HARMONIC_MEAN) + tx, td3, pe, pxi, psf = sampleCubicHermiteCurves([self.px[n3][0][n1y], self.px[n3][n2b][n1y]], + [[-d for d in self.pd3[n3][0][n1y]], + [self.pd1[n3][n2b][n1y][c]-self.pd3[n3][n2b][n1y][c] for c in range(3)]], + self.elementsCountRim+1, lengthFractionStart=1, + arcLengthDerivatives=True) + td1 = interpolateSampleCubicHermite([self.pd1[n3][0][n1y], self.pd3[n3][n2b][n1y]], [[0.0, 0.0, 0.0]] * 2, + pe, pxi, psf)[0] for n2 in range(0, n2b): - self.pd3[n3][n2][m1b] = [-d for d in td3[n2]] + if n2 > 0: + self.px[n3][n2][n1y] = tx[n2] + self.pd1[n3][n2][n1y] = td1[n2] + self.pd3[n3][n2][n1y] = [-d for d in td3[n2]] + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: # left tx = [] @@ -275,23 +288,34 @@ def smoothDerivativesAroundRim(self, n3, n3d=None, rx=0): for ix in range(self.elementsCountAroundFull + 1): n1, n2 = self.convertRimIndex(ix, rx) tx.append(self.px[n3][n2][n1]) - if n2 > self.elementsCountRim: # regular rows - if n1 <= self.elementsCountRim: - td1.append([ -d for d in self.pd2[n3d][n2][n1] ]) - else: - td1.append(self.pd2[n3d][n2][n1]) - else: + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: td1.append(self.pd1[n3d][n2][n1]) - td1 = smoothCubicHermiteDerivativesLine(tx, td1) + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: + if n2 > self.elementsCountRim: # regular rows + if n1 <= self.elementsCountRim: + td1.append([ -d for d in self.pd2[n3d][n2][n1] ]) + else: + td1.append(self.pd2[n3d][n2][n1]) + else: + td1.append(self.pd1[n3d][n2][n1]) + + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: + td1 = smoothCubicHermiteDerivativesLine(tx, td1, fixStartDirection=True, fixEndDirection=True) + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: + td1 = smoothCubicHermiteDerivativesLine(tx, td1) + for ix in range(self.elementsCountAroundFull + 1): n1, n2 = self.convertRimIndex(ix, rx) - if n2 > self.elementsCountRim: # regular rows - if n1 <= self.elementsCountRim: - self.pd2[n3][n2][n1] = [ -d for d in td1[ix] ] - else: - self.pd2[n3][n2][n1] = td1[ix] - else: + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: self.pd1[n3][n2][n1] = td1[ix] + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: + if n2 > self.elementsCountRim: # regular rows + if n1 <= self.elementsCountRim: + self.pd2[n3][n2][n1] = [ -d for d in td1[ix] ] + else: + self.pd2[n3][n2][n1] = td1[ix] + else: + self.pd1[n3][n2][n1] = td1[ix] def generateNodesForOtherHalf(self, mirrorPlane): """ @@ -383,7 +407,9 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes e2a = self.elementsCountRim e2b = self.elementsCountRim + 1 e2c = self.elementsCountRim + 2 - e2d = 2*self.elementsCountUp-1 + e2z = 2*self.elementsCountUp-1-self.elementsCountRim + e2y = e2z - 1 + e2x = e2z - 2 for e3 in range(self.elementsCountAlong): for e2 in range(self.elementsCountUpFull): @@ -396,10 +422,18 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: nids = [self.nodeId[0][e2][e1], self.nodeId[0][e2][e1 + 1], self.nodeId[0][e2 + 1][e1],self.nodeId[0][e2 + 1][e1 + 1], self.nodeId[1][e2][e1], self.nodeId[1][e2][e1 + 1], self.nodeId[1][e2 + 1][e1], self.nodeId[1][e2 + 1][e1 + 1]] - if (e2 < e2b) or (e2 == e2d): + if (e2 < e2b) or (e2 > e2y): if (e1 < e1b) or (e1 > e1y): continue # no element due to triple point closure - if (e2 == e2a) or (e2 == e2d): + if (e2 < e2a) or (e2 > e2z): + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: + if e2 < e2a: + nids = [self.nodeId[e3][e2+1][e1], self.nodeId[e3][e2+1][e1+1], self.nodeId[e3+1][e2+1][e1], self.nodeId[e3+1][e2+1][e1+1], + self.nodeId[e3][e2][e1], self.nodeId[e3][e2][e1+1], self.nodeId[e3+1][e2][e1], self.nodeId[e3+1][e2][e1+1]] + elif e2 > e2z: + nids = [self.nodeId[e3][e2][e1+1], self.nodeId[e3][e2][e1], self.nodeId[e3+1][e2][e1+1], self.nodeId[e3+1][e2][e1], + self.nodeId[e3][e2+1][e1+1], self.nodeId[e3][e2+1][e1], self.nodeId[e3+1][e2+1][e1+1], self.nodeId[e3+1][e2+1][e1]] + elif (e2 == e2a) or (e2 == e2z): # bottom and top row elements if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e2 == e2a: @@ -414,7 +448,7 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes remapEftNodeValueLabel(eft1, [ 2, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) else: remapEftNodeValueLabel(eft1, [ 6, 8 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [] ), ( Node.VALUE_LABEL_D_DS3, [1] ) ]) - elif e2 == e2d: + elif e2 == e2z: eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) scalefactors = [-1.0] @@ -437,7 +471,7 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes else: remapEftNodeValueLabel(eft1, [4, 8], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, [1]),(Node.VALUE_LABEL_D_DS2, [])]) - elif (e2 == e2b) or (e2 == e2d-e2b): + elif (e2 == e2b) or (e2 == e2y): if (e1 <= e1a) or (e1 >= e1z): # map top 2 triple point elements eft1 = tricubichermite.createEftNoCrossDerivatives() @@ -445,12 +479,21 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes scalefactors = [ -1.0 ] if e1 < e1a: e2r = e1 - nids[0] = self.nodeId[0][e2r ][e1b] - nids[1] = self.nodeId[0][e2r + 1][e1b] - nids[4] = self.nodeId[1][e2r ][e1b] - nids[5] = self.nodeId[1][e2r + 1][e1b] - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: + if e2 == e2b: + nids = [self.nodeId[e3][e2c][e1+1], self.nodeId[e3][e2r+1][e1b], self.nodeId[e3+1][e2c][e1+1], self.nodeId[e3+1][e2r+1][e1b], + self.nodeId[e3][e2c][e1], self.nodeId[e3][e2r][e1b], self.nodeId[e3+1][e2c][e1], self.nodeId[e3+1][e2r][e1b]] + if e2 == e2y: + e2r = 2*self.elementsCountUp - e1-1 + nids = [self.nodeId[e3][e2r][e1b], self.nodeId[e3][e2y][e1+1], self.nodeId[e3+1][e2r][e1b], self.nodeId[e3+1][e2y][e1+1], + self.nodeId[e3][e2r+1][e1b], self.nodeId[e3][e2y][e1], self.nodeId[e3+1][e2r+1][e1b], self.nodeId[e3+1][e2y][e1]] + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: + nids[0] = self.nodeId[0][e2r ][e1b] + nids[1] = self.nodeId[0][e2r + 1][e1b] + nids[4] = self.nodeId[1][e2r ][e1b] + nids[5] = self.nodeId[1][e2r + 1][e1b] + remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [] ) ]) elif e1 == e1a: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: if e2 == e2b: @@ -458,9 +501,9 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes nids[2] = self.nodeId[e3+1][e2a][e1b] tripleN = [5, 7] remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS3, [])]) - elif e2 == e2d-e2b: - nids[1] = self.nodeId[e3][e2d+1][e1b] - nids[3] = self.nodeId[e3+1][e2d+1][e1b] + elif e2 == e2y: + nids[1] = self.nodeId[e3][e2z+1][e1b] + nids[3] = self.nodeId[e3+1][e2z+1][e1b] tripleN = [6, 8] remapEftNodeValueLabel(eft1, tripleN, Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, [1]), (Node.VALUE_LABEL_D_DS3, [])]) remapEftNodeValueLabel(eft1, [ 1, 2, 3, 4 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS1, [1] ) ]) @@ -477,9 +520,9 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes nids[4] = self.nodeId[e3][e2a][e1z] nids[6] = self.nodeId[e3+1][e2a][e1z] remapEftNodeValueLabel(eft1, [ 1, 3 ], Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D_DS1, [1] ), ( Node.VALUE_LABEL_D_DS3, [] ) ]) - elif e2 == e2d-e2b: - nids[5] = self.nodeId[e3][e2d+1][e1z] - nids[7] = self.nodeId[e3+1][e2d+1][e1z] + elif e2 == e2y: + nids[5] = self.nodeId[e3][e2z+1][e1z] + nids[7] = self.nodeId[e3+1][e2z+1][e1z] remapEftNodeValueLabel(eft1, [2, 4], Node.VALUE_LABEL_D_DS3,[(Node.VALUE_LABEL_D_DS1, []), (Node.VALUE_LABEL_D_DS3, [])]) elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: nids[1] = self.nodeId[0][e2a][e1z] @@ -489,15 +532,29 @@ def generateElements(self, fieldmodule, coordinates, startElementIdentifier, mes remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS2,[(Node.VALUE_LABEL_D_DS1, [])]) elif e1 > e1z: e2r = self.elementsCountAcross - e1 - nids[0] = self.nodeId[0][e2r ][e1z] - nids[1] = self.nodeId[0][e2r - 1][e1z] - nids[4] = self.nodeId[1][e2r ][e1z] - nids[5] = self.nodeId[1][e2r - 1][e1z] - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) - remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ) ]) + if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: + if e2 == e2b: + nids = [self.nodeId[e3][e2r][e1z], self.nodeId[e3][e2c][e1], self.nodeId[e3+1][e2r][e1z], self.nodeId[e3+1][e2c][e1], + self.nodeId[e3][e2r-1][e1z], self.nodeId[e3][e2c][e1+1], self.nodeId[e3+1][e2r-1][e1z], self.nodeId[e3+1][e2c][e1+1]] + elif e2 == e2y: + e2r = e2z+e1-e1z + nids[1] = self.nodeId[e3][e2r][e1z] + nids[3] = self.nodeId[e3+1][e2r][e1z] + nids[5] = self.nodeId[e3][e2r+1][e1z] + nids[7] = self.nodeId[e3+1][e2r+1][e1z] + elif self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_REGULAR: + nids[0] = self.nodeId[0][e2r ][e1z] + nids[1] = self.nodeId[0][e2r - 1][e1z] + nids[4] = self.nodeId[1][e2r ][e1z] + nids[5] = self.nodeId[1][e2r - 1][e1z] + remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D_DS2, [1] ) ]) + remapEftNodeValueLabel(eft1, [ 1, 2, 5, 6 ], Node.VALUE_LABEL_D_DS2, [ ( Node.VALUE_LABEL_D_DS1, [] ) ]) else: if self._type == ShieldRimDerivativeMode.SHIELD_RIM_DERIVATIVE_MODE_AROUND: - if (e1 <= e1a): + if e1 < e1a: + nids = [ self.nodeId[e3][e2 + 1][e1 + 1], self.nodeId[e3][e2][e1 + 1], self.nodeId[e3+1][e2 + 1][e1 + 1], self.nodeId[e3+1][e2][e1 + 1], + self.nodeId[e3][e2 + 1][e1], self.nodeId[e3][e2][e1], self.nodeId[e3+1][e2 + 1][e1], self.nodeId[e3+1][e2][e1]] + elif e1 == e1a: # map left column elements eft1 = tricubichermite.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) diff --git a/tests/test_cylinder.py b/tests/test_cylinder.py index 73879fb0..b9079f2b 100644 --- a/tests/test_cylinder.py +++ b/tests/test_cylinder.py @@ -19,9 +19,10 @@ def test_cylinder1(self): parameterSetNames = scaffold.getParameterSetNames() self.assertEqual(parameterSetNames, ["Default"]) options = scaffold.getDefaultOptions("Default") - self.assertEqual(9, len(options)) + self.assertEqual(10, len(options)) self.assertEqual(4, options.get("Number of elements across major")) self.assertEqual(4, options.get("Number of elements across minor")) + self.assertEqual(0, options.get("Number of elements across shell")) self.assertEqual(1, options.get("Number of elements along")) context = Context("Test") region = context.getDefaultRegion()