diff --git a/scaffoldmaker/meshtypes/meshtype_3d_centrallinetube1.py b/scaffoldmaker/meshtypes/meshtype_3d_centrallinetube1.py deleted file mode 100644 index 726a37b8..00000000 --- a/scaffoldmaker/meshtypes/meshtype_3d_centrallinetube1.py +++ /dev/null @@ -1,144 +0,0 @@ -""" -Generates a 3-D tubular mesh from a central line -with variable numbers of elements around, along and -through wall, with variable major and minor axis length and -wall thickness. -""" - -from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base -from scaffoldmaker.utils.meshrefinement import MeshRefinement -from scaffoldmaker.utils.tubemesh import * - -class MeshType_3d_centrallinetube1(Scaffold_base): - ''' - Generates a 3-D tubular mesh with variable numbers - of elements around, along the central line, and through wall. - The ellipsoidal tube is created from a central line and - lateral axes data - ''' - @staticmethod - def getName(): - return '3D Central Line Tube 1' - - @staticmethod - def getDefaultOptions(parameterSetName='Default'): - return { - 'Number of elements around' : 8, - 'Number of elements along' : 6, - 'Number of elements through wall' : 1, - 'Tube type' : 1, - 'Use cross derivatives' : False, - 'Use linear through wall' : False, - 'Refine' : False, - 'Refine number of elements around' : 1, - 'Refine number of elements along' : 1, - 'Refine number of elements through wall' : 1 - } - - @staticmethod - def getOrderedOptionNames(): - return [ - 'Number of elements around', - 'Number of elements along', - 'Number of elements through wall', - 'Tube type', - 'Use cross derivatives', - 'Use linear through wall', - 'Refine', - 'Refine number of elements around', - 'Refine number of elements along', - 'Refine number of elements through wall' - ] - - @staticmethod - def checkOptions(options): - for key in [ - 'Number of elements along', - 'Number of elements through wall', - 'Refine number of elements around', - 'Refine number of elements along', - 'Refine number of elements through wall']: - if options[key] < 1: - options[key] = 1 - if (options['Number of elements around'] < 2) : - options['Number of elements around'] = 2 - if (options['Tube type'] < 1 or options['Tube type'] > 3 ) : - options['Tube type'] = 1 - - @staticmethod - def generateBaseMesh(region, options): - """ - Generate the base tricubic Hermite mesh. See also generateMesh(). - :param region: Zinc region to define model in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - :return: None - """ - elementsCountAround = options['Number of elements around'] - elementsCountAlong = options['Number of elements along'] - elementsCountThroughWall = options['Number of elements through wall'] - tubeType = options['Tube type'] - useCrossDerivatives = options['Use cross derivatives'] - useCubicHermiteThroughWall = not(options['Use linear through wall']) - - if tubeType == 1: - # Straight tube - cx = [[1.0, 3.0, 0.0], [ 2.0, 0.0, 4.0 ] ] - cd1 = [[ 1.0, -3.0, 4.0 ], [ 1.0, -3.0, 4.0 ]] - cd2 = [ [ 0.0, 0.3, 0.0 ], [ 0.0, 0.3, 0.0 ]] - cd3 = [ [ 0.0, 0.0, 0.5 ], [ 0.0, 0.0, 0.5 ]] - - # thickness in cd2 and cd3 directions and derivatives (rate of change) - t2 = [ 0.5, 0.5 ] - t2d = [ 0.0, 0.0 ] - t3 = [ 0.5, 0.5 ] - t3d = [ 0.0, 0.0 ] - - elif tubeType == 2: - # Curved tube 1 - cx = [ [ 0.0, 0.0, 0.0], [1.0, 1.0, 0.0], [ 2.0, 0.0, 0.0 ] ] - cd1 = [ [ 1.0, 1.0, 0.0 ], [ 1.0, 0.0, 0.0 ], [1.0, -1.0, 0.0] ] - cd2 = [ [ 0.0, 0.1, 0.0 ], [ 0.0, 0.1, 0.0 ], [ 0.0, 0.1, 0.0 ] ] - cd3 = [ [ 0.0, 0.0, 0.2 ], [ 0.0, 0.0, 0.2 ], [ 0.0, 0.0, 0.2] ] - - # thickness in cd2 and cd3 directions and derivatives (rate of change) - t2 = [ 0.1, 0.1, 0.1 ] - t2d = [ 0.0, 0.0, 0.0 ] - t3 = [ 0.1, 0.1, 0.1 ] - t3d = [ 0.0, 0.0, 0.0 ] - - elif tubeType == 3: - # Curved tube 2 - cx = [ [ 0.0, 0.0, 1.0], [1.0, 1.0, 2.0], [ 2.0, 0.0, 4.0 ] ] - cd1 = [ [ 1.0, 1.0, 1.0 ], [ 1.0, 0.0, 2.0 ], [1.0, -1.0, 2.0] ] - cd2 = [ [ 0.0, 0.2, 0.0 ], [ 0.0, 0.2, 0.0 ], [ 0.0, 0.2, 0.0 ] ] - cd3 = [ [ 0.0, 0.0, 0.2 ], [ 0.0, 0.0, 0.2 ], [ 0.0, 0.0, 0.2] ] - - # thickness in cd2 and cd3 directions and derivatives (rate of change) - t2 = [ 0.1, 0.1, 0.1 ] - t2d = [ 0.0, 0.0, 0.0 ] - t3 = [ 0.1, 0.1, 0.1 ] - t3d = [ 0.0, 0.0, 0.0 ] - - nextNodeIdentifier, nextElementIdentifier = generatetubemesh(region, elementsCountAlong, elementsCountAround, elementsCountThroughWall, - cx, cd1, cd2, cd3, t2, t2d, t3, t3d, useCrossDerivatives, useCubicHermiteThroughWall) - - @classmethod - def generateMesh(cls, region, options): - """ - Generate base or refined mesh. - :param region: Zinc region to create mesh in. Must be empty. - :param options: Dict containing options. See getDefaultOptions(). - """ - if not options['Refine']: - cls.generateBaseMesh(region, options) - return - - refineElementsCountAround = options['Refine number of elements around'] - refineElementsCountAlong = options['Refine number of elements along'] - refineElementsCountThroughWall = options['Refine number of elements through wall'] - - baseRegion = region.createRegion() - cls.generateBaseMesh(baseRegion, options) - - meshrefinement = MeshRefinement(baseRegion, region) - meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) diff --git a/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/scaffoldmaker/meshtypes/meshtype_3d_colon1.py new file mode 100644 index 00000000..160d6556 --- /dev/null +++ b/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -0,0 +1,141 @@ +""" +Generates a 3-D colon mesh along the central line, with variable +numbers of elements around, along and through wall, with +variable radius and thickness along. +""" + +from scaffoldmaker.meshtypes.meshtype_3d_haustra1 import MeshType_3d_haustra1, getColonHaustraSegmentInnerPoints +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils.matrix import * +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from scaffoldmaker.utils.tubemesh import * + +class MeshType_3d_colon1(Scaffold_base): + ''' + Generates a 3-D colon mesh with variable numbers + of elements around, along the central line, and through wall. + The colon is created by a function that generates a haustra + segment and uses tubemesh to map the segment along a central + line profile. + ''' + @staticmethod + def getName(): + return '3D Colon 1' + + @staticmethod + def getParameterSetNames(): + return [ + 'Default', + 'Human 1'] + + @staticmethod + def getDefaultOptions(parameterSetName='Default'): + options = MeshType_3d_haustra1.getDefaultOptions(parameterSetName) + options['Number of elements around'] = 15 + options['Number of elements along haustrum'] = 4 + options['Inner radius'] = 1.0 + options['Haustrum length mid derivative factor'] = 2.0 + options['Wall thickness'] = 0.05 + optionsColon = { + 'Number of haustra segments': 30, + 'Tube type': 2 + } + options.update(optionsColon) + if 'Human 1' in parameterSetName: + options['Tube type'] = 3 + return options + + @staticmethod + def getOrderedOptionNames(): + optionNames = MeshType_3d_haustra1.getOrderedOptionNames() + optionNames.remove('Haustrum length') + for optionName in [ + 'Number of haustra segments', + 'Tube type']: + optionNames.insert(3, optionName) + return optionNames + + def checkOptions(options): + MeshType_3d_haustra1.checkOptions(options) + if options['Number of haustra segments'] < 1: + options['Number of haustra segments'] = 1 + if options['Tube type'] < 1: + options['Tube type'] = 1 + if options['Tube type'] > 3: + options['Tube type'] = 3 + + @staticmethod + def generateBaseMesh(region, options): + """ + Generate the base tricubic Hermite mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: annotationGroups + """ + elementsCountAround = options['Number of elements around'] + elementsCountAlongHaustrum = options['Number of elements along haustrum'] + elementsCountThroughWall = options['Number of elements through wall'] + haustraSegmentCount = options['Number of haustra segments'] + radius = options['Inner radius'] + cornerInnerRadiusFactor = options['Corner inner radius factor'] + haustrumInnerRadiusFactor = options['Haustrum inner radius factor'] + haustrumLengthEndDerivativeFactor = options['Haustrum length end derivative factor'] + haustrumLengthMidDerivativeFactor = options['Haustrum length mid derivative factor'] + wallThickness = options['Wall thickness'] + tubeType = options['Tube type'] + useCrossDerivatives = options['Use cross derivatives'] + useCubicHermiteThroughWall = not(options['Use linear through wall']) + elementsCountAlong = int(elementsCountAlongHaustrum*haustraSegmentCount) + + if tubeType == 1: # Straight tube + cx = [[-4.0, 1.0, 3.0], [ 1.0, 2.0, 0.0 ] ] + cd1 = [[ 5.0, 1.0, -3.0 ], [ 5.0, 1.0, -3.0 ]] + elif tubeType == 2: # Human colon in x-y plane + cx = [ [ 0.0, 0.0, 0.0], [0.0, 10.0, 0.0], [5.0, 9.0, 0.0], [ 10.0, 10.0, 0.0 ], [ 10.0, -2.0, 0.0], [ 7.0, -4.0, 0.0] ] + cd1 = [ [ 0.0, 10.0, 0.0 ], [ 5.0, 5.0, 0.0 ], [5.0, 0.0, 0.0], [ 5.0, -5.0, 0.0 ], [ -3.0, -5.0, 0.0 ], [ -3.0, 0.0, 0.0 ]] + elif tubeType == 3: # Human colon in 3D + cx = [ [ 0.0, 0.0, 0.0], [0.0, 10.0, 3.0], [5.0, 9.0, 0.0], [ 10.0, 10.0, 2.0 ], [15.0, 15.0, 7.0], [ 20.0, -2.0, 0.0], [ 10.0, -4.0, -0.0] ] + cd1 = [ [ 0.0, 10.0, 3.0 ], [ 5.0, 5.0, 0.0 ], [5.0, 0.0, 0.0], [ 10.0, -5.0, 0.0 ], [12.0, 12.0, 0.0], [ 5.0, -12.0, -5.0 ], [ -8.0, 0.0, 0.0 ]] + + # find arclength of colon + length = 0.0 + elementsCountIn = len(cx) - 1 + sd1 = smoothCubicHermiteDerivativesLine(cx, cd1, fixAllDirections = True, + magnitudeScalingMode = DerivativeScalingMode.HARMONIC_MEAN) + for e in range(elementsCountIn): + arcLength = getCubicHermiteArcLength(cx[e], sd1[e], cx[e + 1], sd1[e + 1]) + length += arcLength + haustrumLength = length / haustraSegmentCount + + # Generate inner surface of a haustra segment + xHaustraInner, d1HaustraInner, d2HaustraInner, haustraSegmentAxis = getColonHaustraSegmentInnerPoints(elementsCountAround, elementsCountAlongHaustrum, radius, cornerInnerRadiusFactor, + haustrumInnerRadiusFactor, haustrumLengthEndDerivativeFactor, haustrumLengthMidDerivativeFactor, haustrumLength) + + # Generate tube mesh + annotationGroups, nextNodeIdentifier, nextElementIdentifier = generatetubemesh(region, elementsCountAround, elementsCountAlongHaustrum, elementsCountThroughWall, haustraSegmentCount, + cx, cd1, xHaustraInner, d1HaustraInner, d2HaustraInner, wallThickness, haustraSegmentAxis, haustrumLength, useCrossDerivatives, useCubicHermiteThroughWall) + + return annotationGroups + + @classmethod + def generateMesh(cls, region, options): + """ + Generate base or refined mesh. + :param region: Zinc region to create mesh in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup for mesh. + """ + if not options['Refine']: + cls.generateBaseMesh(region, options) + return + + refineElementsCountAround = options['Refine number of elements around'] + refineElementsCountAlong = options['Refine number of elements along haustrum'] + refineElementsCountThroughWall = options['Refine number of elements through wall'] + + baseRegion = region.createRegion() + baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) + + meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) + return meshrefinement.getAnnotationGroups() diff --git a/scaffoldmaker/meshtypes/meshtype_3d_haustra1.py b/scaffoldmaker/meshtypes/meshtype_3d_haustra1.py new file mode 100644 index 00000000..f3db9857 --- /dev/null +++ b/scaffoldmaker/meshtypes/meshtype_3d_haustra1.py @@ -0,0 +1,345 @@ +""" +Generates a single 3-D haustra segment mesh along a central +line, with variable numbers of elements around, along and +through wall, with variable radius and thickness along. +""" +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.utils.matrix import * +from scaffoldmaker.utils.meshrefinement import MeshRefinement +from scaffoldmaker.utils.tubemesh import * + +class MeshType_3d_haustra1(Scaffold_base): + ''' + Generates a single 3-D haustra segment mesh with variable + numbers of elements around, along the central line, and + through wall. The haustra segment has a triangular profile + with rounded corners at the inter-haustral septa, and a + clover profile in the intra-haustral region. + ''' + @staticmethod + def getName(): + return '3D Haustra 1' + + @staticmethod + def getDefaultOptions(parameterSetName='Default'): + return { + 'Number of elements around' : 9, + 'Number of elements along haustrum' : 3, + 'Number of elements through wall' : 1, + 'Inner radius': 0.5, + 'Corner inner radius factor': 0.5, + 'Haustrum inner radius factor': 0.5, + 'Haustrum length end derivative factor': 0.5, + 'Haustrum length mid derivative factor': 1.0, + 'Wall thickness': 0.01, + 'Haustrum length': 1.0, + 'Use cross derivatives' : False, + 'Use linear through wall' : False, + 'Refine' : False, + 'Refine number of elements around' : 1, + 'Refine number of elements along haustrum' : 1, + 'Refine number of elements through wall' : 1 + } + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Number of elements around', + 'Number of elements along haustrum', + 'Number of elements through wall', + 'Inner radius', + 'Corner inner radius factor', + 'Haustrum inner radius factor', + 'Haustrum length end derivative factor', + 'Haustrum length mid derivative factor', + 'Wall thickness', + 'Haustrum length', + 'Use cross derivatives', + 'Use linear through wall', + 'Refine', + 'Refine number of elements around', + 'Refine number of elements along haustrum', + 'Refine number of elements through wall' + ] + + @staticmethod + def checkOptions(options): + for key in [ + 'Number of elements through wall', + 'Refine number of elements around', + 'Refine number of elements along haustrum', + 'Refine number of elements through wall']: + if options[key] < 1: + options[key] = 1 + if (options['Number of elements around'] < 9) : + options['Number of elements around'] = 9 + if (options['Number of elements around'] % 3 != 0) : + options['Number of elements around'] = (options['Number of elements around']//3)*3 + if (options['Number of elements along haustrum'] < 2) : + options['Number of elements along haustrum'] = 2 + for key in [ + 'Inner radius', + 'Haustrum inner radius factor', + 'Haustrum length end derivative factor', + 'Haustrum length mid derivative factor', + 'Wall thickness', + 'Haustrum length']: + if options[key] < 0.0: + options[key] = 0.0 + if options['Corner inner radius factor'] < 0.1: + options['Corner inner radius factor'] = 0.1 + for key in [ + 'Corner inner radius factor', + 'Haustrum length end derivative factor']: + if options[key] > 1.0: + options[key] = 1.0 + + @staticmethod + def generateBaseMesh(region, options): + """ + Generate the base tricubic Hermite mesh. See also generateMesh(). + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: None + """ + elementsCountAround = options['Number of elements around'] + elementsCountAlongHaustrum = options['Number of elements along haustrum'] + elementsCountThroughWall = options['Number of elements through wall'] + radius = options['Inner radius'] + cornerInnerRadiusFactor = options['Corner inner radius factor'] + haustrumInnerRadiusFactor = options['Haustrum inner radius factor'] + haustrumLengthEndDerivativeFactor = options['Haustrum length end derivative factor'] + haustrumLengthMidDerivativeFactor = options['Haustrum length mid derivative factor'] + wallThickness = options['Wall thickness'] + haustrumLength = options['Haustrum length'] + useCrossDerivatives = options['Use cross derivatives'] + useCubicHermiteThroughWall = not(options['Use linear through wall']) + elementsCountAlong = elementsCountAlongHaustrum + haustraSegmentCount = 1 + + cx = [ [ 0.0, 0.0, 0.0 ], [ haustrumLength, 0.0, 0.0 ] ] + cd1 = [ [ haustrumLength, 0.0, 0.0 ], [ haustrumLength, 0.0, 0.0 ] ] + + # Generate inner surface of a haustra segment + xHaustraInner, d1HaustraInner, d2HaustraInner, haustraSegmentAxis = getColonHaustraSegmentInnerPoints(elementsCountAround, elementsCountAlongHaustrum, radius, cornerInnerRadiusFactor, + haustrumInnerRadiusFactor, haustrumLengthEndDerivativeFactor, haustrumLengthMidDerivativeFactor, haustrumLength) + + # Generate tube mesh + annotationGroups, nextNodeIdentifier, nextElementIdentifier = generatetubemesh(region, elementsCountAround, elementsCountAlongHaustrum, elementsCountThroughWall, haustraSegmentCount, + cx, cd1, xHaustraInner, d1HaustraInner, d2HaustraInner, wallThickness, haustraSegmentAxis, haustrumLength, useCrossDerivatives, useCubicHermiteThroughWall) + + return annotationGroups + + @classmethod + def generateMesh(cls, region, options): + """ + Generate base or refined mesh. + :param region: Zinc region to create mesh in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + """ + if not options['Refine']: + cls.generateBaseMesh(region, options) + return + + refineElementsCountAround = options['Refine number of elements around'] + refineElementsCountAlong = options['Refine number of elements along'] + refineElementsCountThroughWall = options['Refine number of elements through wall'] + + baseRegion = region.createRegion() + baseAnnotationGroups = cls.generateBaseMesh(baseRegion, options) + + meshrefinement = MeshRefinement(baseRegion, region, baseAnnotationGroups) + meshrefinement.refineAllElementsCubeStandard3d(refineElementsCountAround, refineElementsCountAlong, refineElementsCountThroughWall) + return meshrefinement.getAnnotationGroups() + +def getColonHaustraSegmentInnerPoints(elementsCountAround, elementsCountAlongHaustrum, radius, cornerInnerRadiusFactor, + haustrumInnerRadiusFactor, haustrumLengthEndDerivativeFactor, haustrumLengthMidDerivativeFactor, haustrumLength): + """ + Generates a 3-D haustra segment mesh with variable numbers + of elements around, along the central line, and through wall. + The haustra segment has a triangular profile with rounded corners + at the inter-haustral septa, and a clover profile in the intra-haustral + region. + :param elementsCountAround: Number of elements around haustra. + :param elementsCountAlongHaustrum: Number of elements along haustrum. + :param radius: Inner radius defined from center of triangular + profile to vertex of the triangle. + :param cornerInnerRadiusFactor: Roundness of triangular corners of + interhaustral septa. Factor is multiplied by inner radius + to get a radius of curvature at the corners. + :param haustrumInnerRadiusFactor: Factor is multiplied by inner + radius to obtain radius of intersecting circles in the middle cross-section + along a haustra segment. + :param haustrumLengthEndDerivativeFactor: Factor is multiplied by haustrum + length to scale derivative along the end of a haustrum length. + :param haustrumLengthMidDerivativeFactor: Factor is multiplied by haustrum + length to scale derivative along the mid length of the haustrum. + :param haustrumLength: Length of a haustrum. + :return: coordinates, derivatives on inner surface of haustra segment. + """ + + # create nodes + x = [ 0.0, 0.0, 0.0 ] + dx_ds1 = [ 0.0, 0.0, 0.0 ] + dx_ds2 = [ 0.0, 0.0, 0.0 ] + dx_ds3 = [ 0.0, 0.0, 0.0 ] + radiansRangeRC = [7*math.pi/4, 0.0, math.pi/4] + cornerRC = cornerInnerRadiusFactor*radius + unitZ = [0.0, 0.0, 1.0] + haustrumRadius = (haustrumInnerRadiusFactor + 1)*radius + xAround = [] + d1Around = [] + xInner = [] + dx_ds1InnerList = [] + xHaustraSide = [] + xHaustraInner = [] + d1InnerHaustraRaw = [] + xInnerRaw = [] + dx_ds2InnerRaw = [] + xInnerList = [] + dx_ds2InnerList = [] + + # Pre-calculate node locations and derivatives on inner triangle + for n1 in range(3): + radiansAround = n1*2*math.pi / 3 + cosRadiansAround = math.cos(radiansAround) + sinRadiansAround = math.sin(radiansAround) + xc = [(radius - cornerRC) * cosRadiansAround, (radius - cornerRC) * sinRadiansAround, 0.0] + + for n in range(3): + radiansRC = radiansAround + radiansRangeRC[n] + cosRadiansRC = math.cos(radiansRC) + sinRadiansRC = math.sin(radiansRC) + x = [xc[0] + cornerRC*cosRadiansRC, xc[1] + cornerRC*sinRadiansRC, 0.0] + xAround.append(x) + d1 = [ cornerRC*math.pi/4 * -sinRadiansRC, cornerRC*math.pi/4 * cosRadiansRC, 0.0] + d1Around.append(d1) + + xSample = xAround[1:9] + xSample.append(xAround[0]) + xSample.append(xAround[1]) + d1Sample = d1Around[1:9] + d1Sample.append(d1Around[0]) + d1Sample.append(d1Around[1]) + sx, sd1, se, sxi, _= sampleCubicHermiteCurves(xSample, d1Sample, elementsCountAround) + xInner = xInner + sx[0:-1] + d1Inner = smoothCubicHermiteDerivativesLoop(sx[0:-1], sd1[0:-1]) + + # Pre-calculate node locations and derivatives on haustra inner cross-section + elementsCountAroundSide = int(elementsCountAround/3) + Ax = xInner[elementsCountAroundSide][0] + Ay = xInner[elementsCountAroundSide][1] + originRC = (Ax*Ax + Ay*Ay - haustrumRadius*haustrumRadius) / (2*(-Ax - haustrumRadius)) + RC = haustrumRadius - originRC + + if originRC > -Ax: + startTheta = math.asin(Ay/RC) + thetaRC = (math.pi - startTheta)*2 + else: + startTheta = math.pi - math.asin(Ay/RC) + thetaRC = math.asin(Ay/RC)*2 + thetaPerElementAround = thetaRC/(elementsCountAround/3) + + for n in range(elementsCountAroundSide + 1): + theta = startTheta + thetaPerElementAround * n + x = [RC*math.cos(theta) - originRC, + RC*math.sin(theta), + 0.0] + xHaustraSide.append(x) + + ang = [-2/3*math.pi, 0.0, 2/3*math.pi] + + for i in range(3): + rotAng = ang[i] + cosRotAng = math.cos(rotAng) + sinRotAng = math.sin(rotAng) + for n in range(elementsCountAroundSide): + theta = startTheta + thetaPerElementAround * n + cosTheta = math.cos(theta) + sinTheta = math.sin(theta) + x = [ (RC*cosTheta - originRC)*cosRotAng - RC*sinTheta*sinRotAng, + (RC*cosTheta - originRC)*sinRotAng + RC*sinTheta*cosRotAng, + 0.0] + xHaustraInner.append(x) + dx_ds1 = [(-RC*sinTheta*cosRotAng - RC*cosTheta*sinRotAng)*thetaPerElementAround, + (-RC*sinTheta*sinRotAng + RC*cosTheta*cosRotAng)*thetaPerElementAround, + 0.0] + d1InnerHaustraRaw.append(dx_ds1) + d1InnerHaustra = smoothCubicHermiteDerivativesLoop(xHaustraInner, d1InnerHaustraRaw) + + # Sample arclength of haustra segment + for n1 in range(elementsCountAround): + if n1%(elementsCountAround/3) > 0.0: + v1 = [xInner[n1][0], xInner[n1][1], 0.0] + startArcLength = haustrumLengthEndDerivativeFactor * haustrumLength + d1 = [ c*startArcLength for c in unitZ] + v2 = [xHaustraInner[n1][0], xHaustraInner[n1][1], haustrumLength/2] + midArcLength = haustrumLengthMidDerivativeFactor * haustrumLength + d2 = [ c*midArcLength for c in unitZ] + v3 = [xInner[n1][0], xInner[n1][1], haustrumLength] + d3 = [ c*startArcLength for c in unitZ] + else: + v1 = [xInner[n1][0], xInner[n1][1], 0.0] + v2 = [xInner[n1][0], xInner[n1][1], haustrumLength/2] + v3 = [xInner[n1][0], xInner[n1][1], haustrumLength] + d3 = d2 = d1 = [c* haustrumLength/3 for c in unitZ] + nx = [v1, v2, v3] + nd1 = [d1, d2, d3] + sx, sd1, se, sxi, _ = sampleCubicHermiteCurves(nx, nd1, elementsCountAlongHaustrum) + xInnerRaw.append(sx) + dx_ds2InnerRaw.append(sd1) + + # Re-arrange sample order & calculate dx_ds1 and dx_ds3 from dx_ds2 + dx_ds1InnerList = dx_ds1InnerList + d1Inner + for n2 in range(elementsCountAlongHaustrum + 1): + xAround = [] + unitdx_ds1Around = [] + for n1 in range(elementsCountAround): + x = xInnerRaw[n1][n2] + xInnerList.append(x) + dx_ds2 = dx_ds2InnerRaw[n1][n2] + dx_ds2InnerList.append(dx_ds2) + unitTangent = normalise(dx_ds2) + # Inter-haustra + if n2 == 0 or n2 > elementsCountAlongHaustrum - 1: + dx_ds1 = d1Inner[n1] + unitdx_ds3 = crossproduct3(normalise(dx_ds1), unitTangent) + else: + # Intra-Haustra + if elementsCountAlongHaustrum == 2: + unitdx_ds1 = normalise(d1InnerHaustra[n1]) + else: + if n1%(elementsCountAround/3) == 0: # intersection points + unitdx_ds1 = normalise(d1InnerHaustra[n1]) + else: # points on clover + if elementsCountAlongHaustrum > 3: + if n2 < int(elementsCountAlongHaustrum/2): # first half of haustrumLength + axisRot = crossproduct3(unitZ, unitTangent) + elif n2 > int(elementsCountAlongHaustrum/2): # second half of haustrumLength + axisRot = crossproduct3(unitTangent, unitZ) + elif elementsCountAlongHaustrum == 3: # 3 elementsAlongHaustrum + axisRot = crossproduct3(unitTangent, unitZ) + + rotFrame = getRotationMatrixFromAxisAngle(axisRot, math.pi/2) + rotNormal = [rotFrame[j][0]*unitTangent[0] + rotFrame[j][1]*unitTangent[1] + rotFrame[j][2]*unitTangent[2] for j in range(3)] + unitdx_ds3 = normalise(rotNormal) + unitdx_ds1 = crossproduct3(unitTangent, unitdx_ds3) + xAround.append(x) + unitdx_ds1Around.append(unitdx_ds1) + + if n2 > 0 and n2 < elementsCountAlongHaustrum: + dx_ds1InnerAroundList = [] + for n1 in range(elementsCountAround): + v1 = xAround[n1] + d1 = unitdx_ds1Around[n1] + v2 = xAround[(n1+1)%elementsCountAround] + d2 = unitdx_ds1Around[(n1+1)%elementsCountAround] + arcLengthAround = computeCubicHermiteArcLength(v1, d1, v2, d2, True) + dx_ds1 = [c*arcLengthAround for c in d1] + dx_ds1InnerAroundList.append(dx_ds1) + d1Smoothed = smoothCubicHermiteDerivativesLoop(xAround, dx_ds1InnerAroundList) + dx_ds1InnerList = dx_ds1InnerList + d1Smoothed + + dx_ds1InnerList = dx_ds1InnerList + d1Inner + + return xInnerList, dx_ds1InnerList, dx_ds2InnerList, unitZ diff --git a/scaffoldmaker/scaffolds.py b/scaffoldmaker/scaffolds.py index d1004908..4f6ae141 100644 --- a/scaffoldmaker/scaffolds.py +++ b/scaffoldmaker/scaffolds.py @@ -8,7 +8,8 @@ from scaffoldmaker.meshtypes.meshtype_2d_tube1 import MeshType_2d_tube1 from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1 from scaffoldmaker.meshtypes.meshtype_3d_boxhole1 import MeshType_3d_boxhole1 -from scaffoldmaker.meshtypes.meshtype_3d_centrallinetube1 import MeshType_3d_centrallinetube1 +from scaffoldmaker.meshtypes.meshtype_3d_colon1 import MeshType_3d_colon1 +from scaffoldmaker.meshtypes.meshtype_3d_haustra1 import MeshType_3d_haustra1 from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1 from scaffoldmaker.meshtypes.meshtype_3d_heart2 import MeshType_3d_heart2 from scaffoldmaker.meshtypes.meshtype_3d_heartarterialroot1 import MeshType_3d_heartarterialroot1 @@ -37,7 +38,8 @@ def __init__(self): MeshType_2d_tube1, MeshType_3d_box1, MeshType_3d_boxhole1, - MeshType_3d_centrallinetube1, + MeshType_3d_colon1, + MeshType_3d_haustra1, MeshType_3d_heart1, MeshType_3d_heart2, MeshType_3d_heartarterialroot1, diff --git a/scaffoldmaker/utils/matrix.py b/scaffoldmaker/utils/matrix.py new file mode 100644 index 00000000..fdb88ad5 --- /dev/null +++ b/scaffoldmaker/utils/matrix.py @@ -0,0 +1,21 @@ +''' +Utility functions for matrix. +''' + +import math + +def getRotationMatrixFromAxisAngle(rotAxis, theta): + """ + Generate the rotation matrix for rotation about an axis. + :param rotAxis: axis of rotation + :param theta: angle of rotation + :return: rotation matrix + """ + cosTheta = math.cos(theta) + sinTheta = math.sin(theta) + C = 1 - cosTheta + rotMatrix = ([[rotAxis[0]*rotAxis[0]*C + cosTheta, rotAxis[0]*rotAxis[1]*C - rotAxis[2]*sinTheta, rotAxis[0]*rotAxis[2]*C + rotAxis[1]*sinTheta], + [rotAxis[1]*rotAxis[0]*C + rotAxis[2]*sinTheta, rotAxis[1]*rotAxis[1]*C + cosTheta, rotAxis[1]*rotAxis[2]*C - rotAxis[0]*sinTheta], + [rotAxis[2]*rotAxis[0]*C - rotAxis[1]*sinTheta, rotAxis[2]*rotAxis[1]*C + rotAxis[0]*sinTheta, rotAxis[2]*rotAxis[2]*C + cosTheta]]) + + return rotMatrix diff --git a/scaffoldmaker/utils/tubemesh.py b/scaffoldmaker/utils/tubemesh.py index 52eaf8e9..77b1a659 100644 --- a/scaffoldmaker/utils/tubemesh.py +++ b/scaffoldmaker/utils/tubemesh.py @@ -1,26 +1,28 @@ ''' -Utility function for generating tubular mesh from a central line. -Created on Oct 9, 2018 - -@author: Mabelle Lin +Utility function for generating tubular mesh from a central line +using a segment profile. ''' from __future__ import division import math from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear from scaffoldmaker.utils.eftfactory_tricubichermite import eftfactory_tricubichermite -from scaffoldmaker.utils.zinc_utils import * from scaffoldmaker.utils.geometry import * from scaffoldmaker.utils.interpolation import * +from scaffoldmaker.utils.matrix import * from scaffoldmaker.utils.vector import * +from scaffoldmaker.utils.zinc_utils import * from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node def generatetubemesh(region, - elementsCountAlong, elementsCountAround, + elementsCountAlongSegment, elementsCountThroughWall, - cx, cd1, cd2, cd3, t2, t2d, t3, t3d, + segmentCountAlong, + cx, cd1, + xInner, d1Inner, d2Inner, wallThickness, + segmentAxis, segmentLength, useCrossDerivatives, useCubicHermiteThroughWall, # or Zinc Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE etc. nextNodeIdentifier = 1, nextElementIdentifier = 1 @@ -28,34 +30,29 @@ def generatetubemesh(region, ''' Generates a 3-D tubular mesh with variable numbers of elements around, along the central axis, and radially through wall. The - ellipsoidal tubular mesh is created from central line and lateral - axes data - :param elementsCountAlong: number of elements along tube + tubular mesh is created from a segment profile which is mapped onto + the central line and lateral axes data :param elementsCountAround: number of elements around tube + :param elementsCountAlongSegment: number of elements along segment profile :param elementsCountThroughWall: number of elements through wall thickness + :param segmentCountAlong: number of segments along the tube :param cx: coordinates on central line - :param cd1: derivative 1 along central line - :param cd2: derivative 2 to the lateral axis - :param cd3: derivative 3 to the lateral axis perpendicular to direction of cd2 - :param t2: wall thickness in cd2 direction - :param t2d: wall thickness derivative in cd2 direction (rate of change) - :param t3: wall thickness in cd3 direction - :param t3d: wall thickness derivative in cd3 direction (rate of change) + :param cd1: derivative along central line + :param xInner: coordinates on inner surface of segment profile + :param d1Inner: derivatives around inner surface of segment profile + :param d2Inner: derivatives along inner surface of segment profile + :param wallThickness: thickness of wall + :param segmentAxis: axis of segment profile + :param segmentLength: length of segment profile :param useCubicHermiteThroughWall: use linear when false - :param nextNodeIdentifier, nextElementIdentifier: Next identifiers to use and increment. - :return: Final values of nextNodeIdentifier, nextElementIdentifier + :return: annotationGroups, nextNodeIdentifier, nextElementIdentifier ''' - - zero = [ 0.0, 0.0, 0.0 ] + zero = [0.0, 0.0, 0.0] + annotationGroups = [] + elementsCountAlong = elementsCountAlongSegment*segmentCountAlong # Sample central line to get same number of elements as elementsCountAlong sx, sd1, se, sxi, _ = sampleCubicHermiteCurves(cx, cd1, elementsCountAlong) - sd2 = interpolateSampleLinear(cd2, se, sxi) - sd3 = interpolateSampleLinear(cd3, se, sxi) - st2 = interpolateSampleLinear(t2, se, sxi) - st2d = interpolateSampleLinear(t2d, se, sxi) - st3 = interpolateSampleLinear(t3, se, sxi) - st3d = interpolateSampleLinear(t3d, se, sxi) # Find unit normals and binormals at each sample points sNormal = [] @@ -80,7 +77,7 @@ def generatetubemesh(region, if magnitude(cp)> 0.0: axisRot = normalise(cp) thetaRot = math.acos(dotproduct(prevUnitTangent, unitTangent)) - rotFrame = rotationMatrixAboutAxis(axisRot, thetaRot) + rotFrame = getRotationMatrixFromAxisAngle(axisRot, thetaRot) rotNormal = [rotFrame[j][0]*prevUnitNormal[0] + rotFrame[j][1]*prevUnitNormal[1] + rotFrame[j][2]*prevUnitNormal[2] for j in range(3)] unitNormal = normalise(rotNormal) unitBinormal = crossproduct3(unitTangent, unitNormal) @@ -130,90 +127,96 @@ def generatetubemesh(region, dx_ds1 = [ 0.0, 0.0, 0.0 ] dx_ds2 = [ 0.0, 0.0, 0.0 ] dx_ds3 = [ 0.0, 0.0, 0.0 ] + xInnerList = [] + d1InnerList = [] + d2InnerList = [] + d3InnerUnitList = [] + xList = [] + dx_ds1List = [] + dx_ds2List = [] + dx_ds3List = [] - for n3 in range(elementsCountThroughWall + 1): - for n2 in range(elementsCountAlong+1): - aThroughWallElement = sd2[n2][1] + st2[n2]*(n3/elementsCountThroughWall) - bThroughWallElement = sd3[n2][2] + st3[n2]*(n3/elementsCountThroughWall) - perimeterAroundWallElement = getApproximateEllipsePerimeter(aThroughWallElement, bThroughWallElement) - arcLengthPerElementAround = perimeterAroundWallElement / elementsCountAround - prevRadiansAround = updateEllipseAngleByArcLength(aThroughWallElement, bThroughWallElement, 0.0, arcLengthPerElementAround) - st2PerWallElement = st2[n2]/elementsCountThroughWall - st3PerWallElement = st3[n2]/elementsCountThroughWall + # Map each face along segment profile to central line + for nSegment in range(segmentCountAlong): + for nAlongSegment in range(elementsCountAlongSegment+1): + n2 = nSegment*elementsCountAlongSegment + nAlongSegment + if nSegment == 0 or (nSegment > 0 and nAlongSegment > 0): + # Rotate to align segment axis with tangent of central line + segmentMid = [0.0, 0.0, segmentLength/elementsCountAlongSegment* nAlongSegment] + unitTangent = normalise(sd1[n2]) + cp = crossproduct3(segmentAxis, unitTangent) + if magnitude(cp)> 0.0: + axisRot = normalise(cp) + thetaRot = math.acos(dotproduct(segmentAxis, unitTangent)) + rotFrame = getRotationMatrixFromAxisAngle(axisRot, thetaRot) + midRot = [rotFrame[j][0]*segmentMid[0] + rotFrame[j][1]*segmentMid[1] + rotFrame[j][2]*segmentMid[2] for j in range(3)] + translateMatrix = [sx[n2][j] - midRot[j] for j in range(3)] + else: + midRot = segmentMid - if n2 < elementsCountAlong: - aThroughWallElementNext = sd2[n2+1][1] + st2[n2+1]*(n3/elementsCountThroughWall) - bThroughWallElementNext = sd3[n2+1][2] + st3[n2+1]*(n3/elementsCountThroughWall) - perimeterAroundWallElementNext = getApproximateEllipsePerimeter(aThroughWallElementNext, bThroughWallElementNext) - arcLengthPerElementAroundNext = perimeterAroundWallElementNext / elementsCountAround + for n1 in range(elementsCountAround): + n = nAlongSegment*elementsCountAround + n1 + x = xInner[n] + d1 = d1Inner[n] + d2 = d2Inner[n] + if magnitude(cp)> 0.0: + xRot1 = [rotFrame[j][0]*x[0] + rotFrame[j][1]*x[1] + rotFrame[j][2]*x[2] for j in range(3)] + d1Rot1 = [rotFrame[j][0]*d1[0] + rotFrame[j][1]*d1[1] + rotFrame[j][2]*d1[2] for j in range(3)] + d2Rot1 = [rotFrame[j][0]*d2[0] + rotFrame[j][1]*d2[1] + rotFrame[j][2]*d2[2] for j in range(3)] + # Rotate to align first vector on face with binormal axis + if n1 == 0: + firstVector = normalise([xRot1[j] - midRot[j] for j in range(3)]) + thetaRot2 = math.acos(dotproduct(normalise(sBinormal[n2]), firstVector)) + cp2 = crossproduct3(normalise(sBinormal[n2]), firstVector) + if magnitude(cp2) > 0.0: + cp2 = normalise(cp2) + signThetaRot2 = dotproduct(unitTangent, cp2) + axisRot2 = unitTangent + rotFrame2 = getRotationMatrixFromAxisAngle(axisRot2, -signThetaRot2*thetaRot2) + else: + rotFrame2 = [ [1, 0, 0], [0, 1, 0], [0, 0, 1]] + xRot2 = [rotFrame2[j][0]*xRot1[0] + rotFrame2[j][1]*xRot1[1] + rotFrame2[j][2]*xRot1[2] for j in range(3)] + d1Rot2 = [rotFrame2[j][0]*d1Rot1[0] + rotFrame2[j][1]*d1Rot1[1] + rotFrame2[j][2]*d1Rot1[2] for j in range(3)] + d2Rot2 = [rotFrame2[j][0]*d2Rot1[0] + rotFrame2[j][1]*d2Rot1[1] + rotFrame2[j][2]*d2Rot1[2] for j in range(3)] + else: + xRot2 = x + d1Rot2 = d1 + d2Rot2 = d2 - for n1 in range(elementsCountAround): - arcLengthAround = n1*arcLengthPerElementAround - radiansAround = -1*updateEllipseAngleByArcLength(aThroughWallElement, bThroughWallElement, 0.0, arcLengthAround) - cosRadiansAround = math.cos(radiansAround) - sinRadiansAround = math.sin(radiansAround) - x = [sx[n2][j] + aThroughWallElement*cosRadiansAround*sBinormal[n2][j] + bThroughWallElement*sinRadiansAround*sNormal[n2][j] for j in range(3)] - dx_ds1 = [(radiansAround - prevRadiansAround)*(aThroughWallElement*-sinRadiansAround*sBinormal[n2][j] + bThroughWallElement*cosRadiansAround*sNormal[n2][j]) for j in range(3)] + xTranslate = [xRot2[j] + translateMatrix[j] for j in range(3)] - # Calculate curvature to find d1 for node - unitNormal = normalise([aThroughWallElement*cosRadiansAround*sBinormal[n2][j] + bThroughWallElement*sinRadiansAround*sNormal[n2][j] for j in range(3)]) - if n2 == 0: - curvature = getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], unitNormal, 0.0) - elif n2 == elementsCountAlong: - curvature = getCubicHermiteCurvature(sx[n2-1], sd1[n2-1], sx[n2], sd1[n2], unitNormal, 1.0) - else: - curvature = 0.5*( - getCubicHermiteCurvature(sx[n2-1], sd1[n2-1], sx[n2], sd1[n2], unitNormal, 1.0) + - getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], unitNormal, 0.0)) - wallDistance = magnitude([aThroughWallElement*cosRadiansAround*sBinormal[n2][j] + bThroughWallElement*sinRadiansAround*sNormal[n2][j] for j in range(3)]) - factor = 1.0 - curvature*wallDistance - d1Wall = [ factor*c for c in sd1[n2]] + xInnerList.append(xTranslate) + d1InnerList.append(d1Rot2) + d2InnerList.append(d2Rot2) + d3Unit = normalise(crossproduct3(normalise(d1Rot2), normalise(d2Rot2))) + d3InnerUnitList.append(d3Unit) - # Calculate curvature to find d1 for downstream node - if n2 < elementsCountAlong: - arcLengthAroundNext = n1*arcLengthPerElementAroundNext - radiansAroundNext = -1*updateEllipseAngleByArcLength(aThroughWallElementNext, bThroughWallElementNext, 0.0, arcLengthAroundNext) - cosRadiansAroundNext = math.cos(radiansAroundNext) - sinRadiansAroundNext = math.sin(radiansAroundNext) - xNext = [sx[n2+1][j] + aThroughWallElementNext*cosRadiansAroundNext*sBinormal[n2+1][j] + bThroughWallElementNext*sinRadiansAroundNext*sNormal[n2+1][j] for j in range(3)] - unitNormalNext = normalise([aThroughWallElementNext*cosRadiansAroundNext*sBinormal[n2+1][j] + bThroughWallElementNext*sinRadiansAroundNext*sNormal[n2+1][j] for j in range(3)]) - if n2 + 1 == elementsCountAlong: - curvatureNext = getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], unitNormalNext, 1.0) - else: - curvatureNext = 0.5*( - getCubicHermiteCurvature(sx[n2], sd1[n2], sx[n2+1], sd1[n2+1], unitNormalNext, 1.0) + - getCubicHermiteCurvature(sx[n2+1], sd1[n2+1], sx[n2+2], sd1[n2+2], unitNormalNext, 0.0)) - wallDistanceNext = magnitude([aThroughWallElementNext*cosRadiansAroundNext*sBinormal[n2+1][j] + bThroughWallElementNext*sinRadiansAroundNext*sNormal[n2+1][j] for j in range(3)]) - factorNext = 1.0 - curvatureNext*wallDistanceNext - d1WallNext = [ factorNext*c for c in sd1[n2+1] ] - arcLength = computeCubicHermiteArcLength(x, d1Wall, xNext, d1WallNext, True) - dx_ds2 = [arcLength*c for c in normalise(d1Wall)] - if n2 == elementsCountAlong - 1: - secondLastX = x - secondLastd1Wall = d1Wall - lastX = xNext - lastd1Wall = d1WallNext - else: - arcLength = computeCubicHermiteArcLength(secondLastX, secondLastd1Wall, lastX, lastd1Wall, True) - dx_ds2 = [arcLength*c for c in normalise(lastd1Wall)] + # Pre-calculate node locations and derivatives on outer boundary + xOuterList, curvatureInner = getOuterCoordinatesAndCurvatureFromInner(xInnerList, d1InnerList, d3InnerUnitList, wallThickness, elementsCountAlong, elementsCountAround) - dx_ds3 = [ st2PerWallElement*cosRadiansAround*sBinormal[n2][j] + st3PerWallElement*sinRadiansAround*sNormal[n2][j] for j in range(3)] # Modify later to calculate with interpolation + # Interpolate to get nodes through wall + for n3 in range(elementsCountThroughWall + 1): + xi3 = 1/elementsCountThroughWall * n3 + x, dx_ds1, dx_ds2, dx_ds3 = interpolatefromInnerAndOuter( xInnerList, xOuterList, wallThickness, xi3, curvatureInner, d1InnerList, d2InnerList, d3InnerUnitList, elementsCountAround, elementsCountAlong, elementsCountThroughWall) + xList = xList + x + dx_ds1List = dx_ds1List + dx_ds1 + dx_ds2List = dx_ds2List + dx_ds2 + dx_ds3List = dx_ds3List + dx_ds3 - node = nodes.createNode(nodeIdentifier, nodetemplate) - cache.setNode(node) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) - if useCubicHermiteThroughWall: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) - if useCrossDerivatives: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) - if useCubicHermiteThroughWall: - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) - nodeIdentifier = nodeIdentifier + 1 - prevRadiansAround = radiansAround + for n in range(len(xList)): + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, xList[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1List[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2List[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3List[n]) + if useCrossDerivatives: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS2DS3, 1, zero) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1, zero) + # print('NodeIdentifier = ', nodeIdentifier, xList[n]) + nodeIdentifier = nodeIdentifier + 1 # # For debugging - Nodes along central line # for pt in range(len(sx)): @@ -242,19 +245,94 @@ def generatetubemesh(region, fm.endChange() - return nodeIdentifier, elementIdentifier + return annotationGroups, nodeIdentifier, elementIdentifier -def rotationMatrixAboutAxis(rotAxis, theta): +def getOuterCoordinatesAndCurvatureFromInner(xInner, d1Inner, d3Inner, wallThickness, elementsCountAlong, elementsCountAround): + """ + Generates coordinates on outer surface and curvature of inner + surface from coordinates and derivatives of inner surface using + wall thickness and normals. + param xInner: Coordinates on inner surface + param d1Inner: Derivatives on inner surface around tube + param d3Inner: Derivatives on inner surface through wall + param wallThickness: Thickness of wall + param elementsCountAlong: Number of elements along tube + param elementsCountAround: Number of elements around tube + return xOuter: Coordinates on outer surface + return curvatureInner: Curvature of coordinates on inner surface """ - Generate the rotation matrix for rotation about an axis. - :param rotAxis: axis of rotation - :param theta: angle of rotation - :return: rotation matrix + xOuter = [] + curvatureInner = [] + for n2 in range(elementsCountAlong + 1): + for n1 in range(elementsCountAround): + n = n2*elementsCountAround + n1 + x = [xInner[n][i] + d3Inner[n][i]*wallThickness for i in range(3)] + prevIdx = n-1 if (n1 != 0) else (n2+1)*elementsCountAround - 1 + nextIdx = n+1 if (n1 < elementsCountAround-1) else n2*elementsCountAround + norm = d3Inner[n] + curvatureAround = 0.5*( + getCubicHermiteCurvature(xInner[prevIdx], d1Inner[prevIdx], xInner[n], d1Inner[n], norm, 1.0) + + getCubicHermiteCurvature(xInner[n], d1Inner[n], xInner[nextIdx], d1Inner[nextIdx], norm, 0.0)) + xOuter.append(x) + curvatureInner.append(curvatureAround) + + return xOuter, curvatureInner + +def interpolatefromInnerAndOuter( xInner, xOuter, thickness, xi3, curvatureInner, d1Inner, d2Inner, d3InnerUnit, + elementsCountAround, elementsCountAlong, elementsCountThroughWall): """ - cosTheta = math.cos(theta) - sinTheta = math.sin(theta) - C = 1 - cosTheta - rotMatrix = ([[rotAxis[0]*rotAxis[0]*C + cosTheta, rotAxis[0]*rotAxis[1]*C - rotAxis[2]*sinTheta, rotAxis[0]*rotAxis[2]*C + rotAxis[1]*sinTheta], - [rotAxis[1]*rotAxis[0]*C + rotAxis[2]*sinTheta, rotAxis[1]*rotAxis[1]*C + cosTheta, rotAxis[1]*rotAxis[2]*C - rotAxis[0]*sinTheta], - [rotAxis[2]*rotAxis[0]*C - rotAxis[1]*sinTheta, rotAxis[2]*rotAxis[1]*C + rotAxis[0]*sinTheta, rotAxis[2]*rotAxis[2]*C + cosTheta]]) - return rotMatrix + Generate coordinates and derivatives at xi3 by interpolating with + inner and outer coordinates and derivatives. + param xInner: Coordinates on inner surface + param xOuter: Coordinates on outer surface + param thickness: Thickness of wall + param curvatureInner: Curvature of coordinates on inner surface + param d1Inner: Derivatives on inner surface around tube + param d2Inner: Derivatives on inner surface along tube + param d3InnerUnit: Unit derivatives on inner surface through wall + param elementsCountAround: Number of elements around tube + param elementsCountAlong: Number of elements along tube + param elementsCountThroughWall: Number of elements through wall + return xList, dx_ds1List, dx_ds2List, dx_ds3List: Coordinates and derivatives on xi3 + """ + xList = [] + dx_ds1List = [] + dx_ds2List = [] + dx_ds3List =[] + + for n2 in range(elementsCountAlong+1): + for n1 in range(elementsCountAround): + n = n2*elementsCountAround + n1 + norm = d3InnerUnit[n] + # x + innerx = xInner[n] + outerx = xOuter[n] + dWall = [thickness*c for c in norm] + x = interpolateCubicHermite(innerx, dWall, outerx, dWall, xi3) + xList.append(x) + # dx_ds1 + factor = 1.0 - curvatureInner[n]*thickness*xi3 + dx_ds1 = [ factor*c for c in d1Inner[n]] + dx_ds1List.append(dx_ds1) + # dx_ds2 + if n2 > 0 and n2 < elementsCountAlong: + prevIdx = (n2-1)*elementsCountAround + n1 + nextIdx = (n2+1)*elementsCountAround + n1 + curvatureAround = 0.5*( + getCubicHermiteCurvature(xInner[prevIdx], d2Inner[prevIdx], xInner[n], d2Inner[n], norm, 1.0)+ + getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[nextIdx], d2Inner[nextIdx], norm, 0.0)) + elif n2 == 0: + nextIdx = (n2+1)*elementsCountAround + n1 + curvatureAround = getCubicHermiteCurvature(xInner[n], d2Inner[n], xInner[nextIdx], d2Inner[nextIdx], norm, 0.0) + else: + prevIdx = (n2-1)*elementsCountAround + n1 + curvatureAround = getCubicHermiteCurvature(xInner[prevIdx], d2Inner[prevIdx], xInner[n], d2Inner[n], norm, 1.0) + + factor = 1.0 - curvatureAround*thickness*xi3 + dx_ds2 = [ factor*c for c in d2Inner[n]] + dx_ds2List.append(dx_ds2) + #dx_ds3 + dx_ds3 = [c * thickness/elementsCountThroughWall for c in norm] + dx_ds3List.append(dx_ds3) + + return xList, dx_ds1List, dx_ds2List, dx_ds3List