diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index a5d1b61a..fa68b8a0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -103,7 +103,7 @@ class MeshType_3d_bladderurethra1(Scaffold_base): }) } ostiumDefaultScaffoldPackages = { - 'Ostium Cat 1': ScaffoldPackage(MeshType_3d_ostium1, { + 'Ureter Cat 1': ScaffoldPackage(MeshType_3d_ostium1, { 'scaffoldSettings': { 'Number of vessels': 1, 'Number of elements around ostium': 8, # implemented for 8 @@ -127,7 +127,7 @@ class MeshType_3d_bladderurethra1(Scaffold_base): 'Refine number of elements through wall': 1 }, }), - 'Ostium Rat 1': ScaffoldPackage(MeshType_3d_ostium1, { + 'Ureter Rat 1': ScaffoldPackage(MeshType_3d_ostium1, { 'scaffoldSettings': { 'Number of vessels': 1, 'Number of elements across common': 2, @@ -172,9 +172,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): else: centralPathOption_LUT = cls.centralPathDefaultScaffoldPackages_LUT['Cat 1'] if 'Rat 1' in parameterSetName: - ostiumOption = cls.ostiumDefaultScaffoldPackages['Ostium Rat 1'] + ureterOption = cls.ostiumDefaultScaffoldPackages['Ureter Rat 1'] else: - ostiumOption = cls.ostiumDefaultScaffoldPackages['Ostium Cat 1'] + ureterOption = cls.ostiumDefaultScaffoldPackages['Ureter Cat 1'] options = { 'Central path LUT': copy.deepcopy(centralPathOption_LUT), 'Number of elements along bladder': 12, @@ -192,10 +192,10 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Urethra diameter 2': 1.0, 'Urethra wall thickness': 0.5, 'Include ureter': False, - 'Ureter': copy.deepcopy(ostiumOption), - 'Ostium position around': 0.65, # should be on the dorsal part (> 0.5) - 'Ostium position down': 0.75, - 'Number of elements radially on annulus': 1, + 'Ureter': copy.deepcopy(ureterOption), + 'Ureter position around': 0.67, # should be on the dorsal part (> 0.5) + 'Ureter position down': 0.83, + 'Number of elements ureter radial': 2, 'Use cross derivatives': False, 'Use linear through wall': True, 'Refine': False, @@ -209,8 +209,8 @@ def getDefaultOptions(cls, parameterSetName='Default'): options['Wall thickness'] = 0.2 options['Neck diameter 1'] = 3.5 options['Neck diameter 2'] = 2.0 - options['Ostium position around'] = 0.65 # should be on the dorsal part (> 0.5) - options['Ostium position down'] = 0.75 + options['Ureter position around'] = 0.67 # should be on the dorsal part (> 0.5) + options['Ureter position down'] = 0.83 options['Urethra diameter 1'] = 0.75 options['Urethra diameter 2'] = 0.65 options['Urethra wall thickness'] = 0.25 @@ -236,8 +236,9 @@ def getOrderedOptionNames(): 'Urethra wall thickness', 'Include ureter', 'Ureter', - 'Ostium position around', - 'Ostium position down', + 'Ureter position around', + 'Ureter position down', + 'Number of elements ureter radial', 'Use cross derivatives', 'Use linear through wall', 'Refine', @@ -292,6 +293,7 @@ def checkOptions(cls, options): 'Number of elements around', 'Number of elements through wall', 'Number of elements along urethra', + 'Number of elements ureter radial', 'Refine number of elements along', 'Refine number of elements around', 'Refine number of elements through wall']: @@ -300,6 +302,14 @@ def checkOptions(cls, options): if options['Number of elements through wall'] > 1: if options['Include ureter']: options['Number of elements through wall'] = 1 + if options['Ureter position around'] < 0.1: + options['Ureter position around'] = 0.1 + elif options['Ureter position around'] > 0.9: + options['Ureter position around'] = 0.9 + if options['Ureter position down'] < 0.15: + options['Ureter position down'] = 0.15 + elif options['Ureter position down'] > 0.95: + options['Ureter position down'] = 0.95 @classmethod def updateSubScaffoldOptions(cls, options): @@ -307,9 +317,9 @@ def updateSubScaffoldOptions(cls, options): Update ostium sub-scaffold options which depend on parent options. ''' bladderWallThickness = options['Wall thickness'] - ostiumOptions = options['Ureter'] - ostiumDefaultOptions = ostiumOptions.getScaffoldSettings() - ostiumDefaultOptions['Ostium wall thickness'] = bladderWallThickness + ureterOptions = options['Ureter'] + ureterDefaultOptions = ureterOptions.getScaffoldSettings() + ureterDefaultOptions['Ostium wall thickness'] = bladderWallThickness @classmethod def generateBaseMesh(cls, region, options): @@ -340,12 +350,12 @@ def generateBaseMesh(cls, region, options): urethraWallThickness = options['Urethra wall thickness'] includeUreter = options['Include ureter'] - ostiumOptions = options['Ureter'] - ostiumDefaultOptions = ostiumOptions.getScaffoldSettings() - elementsCountAroundOstium = ostiumDefaultOptions['Number of elements around ostium'] - elementsCountAnnulusRadially = options['Number of elements radially on annulus'] - ostiumPositionAround = options['Ostium position around'] - ostiumPositionDown = options['Ostium position down'] + ureterOptions = options['Ureter'] + ureterDefaultOptions = ureterOptions.getScaffoldSettings() + elementsCountAroundUreter = ureterDefaultOptions['Number of elements around ostium'] + elementsCountUreterRadial = options['Number of elements ureter radial'] + ureterPositionAround = options['Ureter position around'] + ureterPositionDown = options['Ureter position down'] firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -638,9 +648,9 @@ def generateBaseMesh(cls, region, options): # Obtain elements count along body and neck of the bladder for defining annotation groups if includeUrethra: # bladderLength = length - urethraLength - bodyLength = ostiumPositionDown * bladderLength + bodyLength = ureterPositionDown * bladderLength else: - bodyLength = ostiumPositionDown * length + bodyLength = ureterPositionDown * length elementsCountAlongBody = int(bodyLength / bladderSegmentLength) elementsCountAlongNeck = elementsCountAlongBladder - elementsCountAlongBody @@ -696,7 +706,7 @@ def generateBaseMesh(cls, region, options): if includeUreter: elementsCount1 = elementsCountAround // 2 elementsCount2 = elementsCountAlongBladder - # Create trackSurface at the outer layer of the bladder for ostium 1 + # Create trackSurface at the outer layer of the bladder for ureter 1 nodesOnTrackSurface_x = [] nodesOnTrackSurface_d1 = [] nodesOnTrackSurface_d2 = [] @@ -705,14 +715,14 @@ def generateBaseMesh(cls, region, options): nodesOnTrackSurface_x.append(outerNodes_x[n2 * elementsCountAround + n1]) nodesOnTrackSurface_d1.append(outerNodes_d1[n2 * elementsCountAround + n1]) nodesOnTrackSurface_d2.append(outerNodes_d2[n2 * elementsCountAround + n1]) - trackSurfaceOstium1 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface_x, + trackSurfaceUreter1 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface_x, nodesOnTrackSurface_d1, nodesOnTrackSurface_d2) - ostium1Position = trackSurfaceOstium1.createPositionProportion(ostiumPositionAround, ostiumPositionDown) - ostium1Position.xi1 = 1.0 - ostium1Position.xi2 = 1.0 - ostiumElementPositionAround = ostium1Position.e1 - ostiumElementPositionDown = ostium1Position.e2 - # Create trackSurface at the outer layer of the bladder for ostium 2 + + ureter1Position = trackSurfaceUreter1.createPositionProportion(ureterPositionAround, ureterPositionDown) + ureterElementPositionAround = ureter1Position.e1 + ureterElementPositionDown = ureter1Position.e2 + + # Create trackSurface at the outer layer of the bladder for ureter 2 nodesOnTrackSurface2_x = [] nodesOnTrackSurface2_d1 = [] nodesOnTrackSurface2_d2 = [] @@ -724,18 +734,19 @@ def generateBaseMesh(cls, region, options): nodesOnTrackSurface2_x.append(outerNodes_x[n2 * elementsCountAround]) nodesOnTrackSurface2_d1.append(outerNodes_d1[n2 * elementsCountAround]) nodesOnTrackSurface2_d2.append(outerNodes_d2[n2 * elementsCountAround]) - trackSurfaceOstium2 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface2_x, + + trackSurfaceUreter2 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface2_x, nodesOnTrackSurface2_d1, nodesOnTrackSurface2_d2) - ostium2Position = TrackSurfacePosition(elementsCountAround - ostiumElementPositionAround, - ostiumElementPositionDown - 1, 0.0, 1.0) + ureter2Position = TrackSurfacePosition(elementsCountAround//2 - ureterElementPositionAround + (-1 if ureter1Position.xi1 > 0 else 0), + ureterElementPositionDown, + (1 - ureter1Position.xi1) if ureter1Position.xi1 > 0 else ureter1Position.xi1, + ureter1Position.xi2) annulusMeshGroups = [neckMeshGroup, urinaryBladderMeshGroup] - generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaultOptions, - elementsCountAround, elementsCountAroundOstium, trackSurfaceOstium1, - ostium1Position, trackSurfaceOstium2, ostium2Position, - ostiumElementPositionDown, ostiumElementPositionAround, - xFinal, d1Final, d2Final, nextNodeIdentifier, - nextElementIdentifier, elementsCountAnnulusRadially, - annulusMeshGroups) + generateUreters(region, nodes, mesh, ureterDefaultOptions, elementsCountAround, elementsCountThroughWall, + elementsCountAroundUreter, trackSurfaceUreter1, ureter1Position, trackSurfaceUreter2, + ureter2Position, ureterElementPositionDown, ureterElementPositionAround, xFinal, d1Final, + d2Final, nextNodeIdentifier, nextElementIdentifier, elementsCountUreterRadial, + annulusMeshGroups) fm.endChange() return annotationGroups @@ -807,88 +818,89 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): lumenOfUrethra = findOrCreateAnnotationGroupForTerm(annotationGroups, region, get_bladder_term("lumen of urethra")) lumenOfUrethra.getMeshGroup(mesh2d).addElementsConditional(is_urethra_lumen) -def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaultOptions,elementsCountAround, - elementsCountAroundOstium, trackSurfaceOstium1, ostium1Position, - trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, - ostiumElementPositionAround, xBladder, d1Bladder, d2Bladder, - nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadially, - annulusMeshGroups = []): +def generateUreters(region, nodes, mesh, ureterDefaultOptions,elementsCountAround, elementsCountThroughWall, + elementsCountAroundUreter, trackSurfaceUreter1, ureter1Position, trackSurfaceUreter2, + ureter2Position, ureterElementPositionDown, ureterElementPositionAround, + xBladder, d1Bladder, d2Bladder, nextNodeIdentifier, nextElementIdentifier, + elementsCountUreterRadial, annulusMeshGroups = []): # Create ureters on the surface # Ureter 1 - centerUreter1_x, centerUreter1_d1, centerUreter1_d2 = trackSurfaceOstium1.evaluateCoordinates(ostium1Position, + centerUreter1_x, centerUreter1_d1, centerUreter1_d2 = trackSurfaceUreter1.evaluateCoordinates(ureter1Position, derivatives=True) td1, td2, td3 = calculate_surface_axes(centerUreter1_d1, centerUreter1_d2, [1.0, 0.0, 0.0]) - m1 = 2 * elementsCountAround * (ostiumElementPositionDown - 1) + ostiumElementPositionAround + 2 - ureter1StartCornerx = xBladder[m1] + endPointStartId1 = elementsCountThroughWall + 1 \ + + (elementsCountThroughWall + 1) * elementsCountAround * (ureterElementPositionDown - (1 if ureter1Position.xi2 > 0.5 else 2)) \ + + ureterElementPositionAround + (1 if ureter1Position.xi1 > 0.5 else 0) + ureter1StartCornerx = xBladder[endPointStartId1 - 1] v1 = [(ureter1StartCornerx[c] - centerUreter1_x[c]) for c in range(3)] - ostium1Direction = vector.crossproduct3(td3, v1) + ureter1Direction = vector.crossproduct3(td3, v1) nodeIdentifier, elementIdentifier, (o1_x, o1_d1, o1_d2, _, o1_NodeId, o1_Positions) = \ - generateOstiumMesh(region, ostiumDefaultOptions, trackSurfaceOstium1, ostium1Position, ostium1Direction, + generateOstiumMesh(region, ureterDefaultOptions, trackSurfaceUreter1, ureter1Position, ureter1Direction, startNodeIdentifier=nextNodeIdentifier, startElementIdentifier=nextElementIdentifier, ostiumMeshGroups=None) + # Ureter 2 - centerUreter2_x, centerUreter2_d1, centerUreter2_d2 = trackSurfaceOstium2.evaluateCoordinates(ostium2Position, + centerUreter2_x, centerUreter2_d1, centerUreter2_d2 = trackSurfaceUreter2.evaluateCoordinates(ureter2Position, derivatives=True) ad1, ad2, ad3 = calculate_surface_axes(centerUreter2_d1, centerUreter2_d2, [1.0, 0.0, 0.0]) - m2 = 2 * elementsCountAround * (ostiumElementPositionDown - 1) + elementsCountAround - ostiumElementPositionAround - ureter2StartCornerx = xBladder[m2] + endPointStartId2 = elementsCountThroughWall + 1 \ + + (elementsCountThroughWall + 1) * elementsCountAround * \ + (ureterElementPositionDown - (1 if ureter1Position.xi2 > 0.5 else 2)) \ + + elementsCountAround - ureterElementPositionAround + (-1 if ureter1Position.xi1 > 0.5 else 0) + ureter2StartCornerx = xBladder[endPointStartId2 - 1] v2 = [(ureter2StartCornerx[c] - centerUreter2_x[c]) for c in range(3)] - ostium2Direction = vector.crossproduct3(ad3, v2) + ureter2Direction = vector.crossproduct3(ad3, v2) nodeIdentifier, elementIdentifier, (o2_x, o2_d1, o2_d2, _, o2_NodeId, o2_Positions) = \ - generateOstiumMesh(region, ostiumDefaultOptions, trackSurfaceOstium2, ostium2Position, ostium2Direction, + generateOstiumMesh(region, ureterDefaultOptions, trackSurfaceUreter2, ureter2Position, ureter2Direction, startNodeIdentifier=nodeIdentifier, startElementIdentifier=elementIdentifier) - # Create annulus mesh around ostiums - endPoints1_x = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endPoints1_d1 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endPoints1_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endNode1_Id = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endDerivativesMap = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endPoints2_x = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endPoints2_d1 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endPoints2_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - endNode2_Id = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - - for n3 in range(2): - n1 = 0 - endNode1_Id[n3][n1] = (2 * ostiumElementPositionDown - 2 + n3) * elementsCountAround + \ - ostiumElementPositionAround + 1 + 2 - endNode1_Id[n3][n1 + 1] = endNode1_Id[n3][n1] + 2 * elementsCountAround - endNode1_Id[n3][n1 + 2] = endNode1_Id[n3][n1 + 1] + 2 * elementsCountAround - endNode1_Id[n3][n1 + 3] = endNode1_Id[n3][n1 + 2] + 1 - endNode1_Id[n3][n1 + 4] = endNode1_Id[n3][n1 + 3] + 1 - endNode1_Id[n3][n1 + 5] = endNode1_Id[n3][n1 + 1] + 2 - endNode1_Id[n3][n1 + 6] = endNode1_Id[n3][n1] + 2 - endNode1_Id[n3][n1 + 7] = endNode1_Id[n3][n1] + 1 - - endNode2_Id[n3][n1] = (2 * ostiumElementPositionDown - 2 + n3) * elementsCountAround + elementsCountAround - \ - ostiumElementPositionAround - 1 + 2 - endNode2_Id[n3][n1 + 1] = endNode2_Id[n3][n1] + 2 * elementsCountAround - endNode2_Id[n3][n1 + 2] = endNode2_Id[n3][n1 + 1] + 2 * elementsCountAround - endNode2_Id[n3][n1 + 3] = endNode2_Id[n3][n1 + 2] + 1 - endNode2_Id[n3][n1 + 7] = endNode2_Id[n3][n1] + 1 - if ostiumElementPositionAround == 0: - endNode2_Id[n3][n1 + 4] = endNode2_Id[n3][n1 + 3] - elementsCountAround + 1 - endNode2_Id[n3][n1 + 5] = endNode2_Id[n3][n1 + 4] - 2 * elementsCountAround - endNode2_Id[n3][n1 + 6] = endNode2_Id[n3][n1 + 5] - 2 * elementsCountAround - else: - endNode2_Id[n3][n1 + 4] = endNode2_Id[n3][n1 + 3] + 1 - endNode2_Id[n3][n1 + 5] = endNode2_Id[n3][n1 + 1] + 2 - endNode2_Id[n3][n1 + 6] = endNode2_Id[n3][n1] + 2 + # Create annulus mesh around ureters + endPoints1_x = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endPoints1_d1 = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endPoints1_d2 = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endNode1_Id = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endDerivativesMap = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endPoints2_x = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endPoints2_d1 = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endPoints2_d2 = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + endNode2_Id = [[None] * elementsCountAroundUreter, [None] * elementsCountAroundUreter] + + count = 0 + for n2 in range(3): + endNode1_Id[0][count] = endPointStartId1 + n2 * (elementsCountThroughWall + 1) * elementsCountAround + endNode2_Id[0][count] = endPointStartId2 + n2 * (elementsCountThroughWall + 1) * elementsCountAround + count += 1 + for n1 in range(2): + endNode1_Id[0][count] = endNode1_Id[0][count - 1] + 1 + endNode2_Id[0][count] = endNode2_Id[0][count - 1] + 1 - \ + (0 if (endNode2_Id[0][count - 1] - 2) % elementsCountAround > 0 + else elementsCountAround) + count += 1 + for n2 in range(2): + endNode1_Id[0][count] = endNode1_Id[0][count - 1] - (elementsCountThroughWall + 1) * elementsCountAround + endNode2_Id[0][count] = endNode2_Id[0][count - 1] - (elementsCountThroughWall + 1) * elementsCountAround + count += 1 + endNode1_Id[0][count] = endNode1_Id[0][count - 1] - 1 + endNode2_Id[0][count] = endNode2_Id[0][count - 1] - 1 + \ + (0 if (endNode2_Id[0][count - 1] - 3) % elementsCountAround > 0 else elementsCountAround) + + for n in range(len(endNode1_Id[0])): + endNode1_Id[1][n] = endNode1_Id[0][n] + elementsCountAround + endNode2_Id[1][n] = endNode2_Id[0][n] + elementsCountAround for n3 in range(2): - for n1 in range(elementsCountAroundOstium): + for n1 in range(elementsCountAroundUreter): nc1 = endNode1_Id[n3][n1] - 1 endPoints1_x[n3][n1] = xBladder[nc1] endPoints1_d1[n3][n1] = d1Bladder[nc1] - endPoints1_d2[n3][n1] = [d2Bladder[nc1][c] for c in range(3)] + endPoints1_d2[n3][n1] = d2Bladder[nc1] nc2 = endNode2_Id[n3][n1] - 1 endPoints2_x[n3][n1] = xBladder[nc2] endPoints2_d1[n3][n1] = d1Bladder[nc2] endPoints2_d2[n3][n1] = d2Bladder[nc2] - for n1 in range(elementsCountAroundOstium): + for n1 in range(elementsCountAroundUreter): if n1 == 0: endDerivativesMap[0][n1] = ((-1, 0, 0), (-1, -1, 0), None, (0, 1, 0)) endDerivativesMap[1][n1] = ((-1, 0, 0), (-1, -1, 0), None, (0, 1, 0)) @@ -914,30 +926,76 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul endDerivativesMap[0][n1] = ((-1, 0, 0), (0, -1, 0), None) endDerivativesMap[1][n1] = ((-1, 0, 0), (0, -1, 0), None) + startProportions1 = [] + for n in range(len(o1_Positions)): + startProportions1.append(trackSurfaceUreter1.getProportion(o1_Positions[n])) + + startProportions2 = [] + for n in range(len(o2_Positions)): + startProportions2.append(trackSurfaceUreter2.getProportion(o2_Positions[n])) + + endProportions1 = [] + elementsAroundTrackSurface1 = trackSurfaceUreter1.elementsCount1 + elementsAlongTrackSurface1 = trackSurfaceUreter1.elementsCount2 + + endProportions2 = [] + elementsAroundTrackSurface2 = trackSurfaceUreter2.elementsCount1 + elementsAlongTrackSurface2 = trackSurfaceUreter2.elementsCount2 + + firstIdxAround1 = ureterElementPositionAround + (0 if ureter1Position.xi1 > 0.5 else -1) + firstIdxAlong = ureterElementPositionDown - (0 if ureter1Position.xi2 > 0.5 else 1) + firstIdxAround2 = elementsCountAround//2 - ureterElementPositionAround - (2 if ureter1Position.xi1 > 0.5 else 1) + + for n in range(3): + endProportions1.append([(firstIdxAround1)/elementsAroundTrackSurface1, + (firstIdxAlong + n)/elementsAlongTrackSurface1]) + endProportions2.append([(firstIdxAround2) / elementsAroundTrackSurface2, + (firstIdxAlong + n) / elementsAlongTrackSurface2]) + for n in range(2): + endProportions1.append([(firstIdxAround1 + n + 1) / elementsAroundTrackSurface1, + (firstIdxAlong + 2) / elementsAlongTrackSurface1]) + endProportions2.append([(firstIdxAround2 + n + 1) / elementsAroundTrackSurface2, + (firstIdxAlong + 2) / elementsAlongTrackSurface2]) + for n in range(2): + endProportions1.append([(firstIdxAround1 + 2) / elementsAroundTrackSurface1, + (firstIdxAlong - n + 1) / elementsAlongTrackSurface1]) + endProportions2.append([(firstIdxAround2 + 2) / elementsAroundTrackSurface2, + (firstIdxAlong - n + 1) / elementsAlongTrackSurface2]) + endProportions1.append([(firstIdxAround1 + 1 )/ elementsAroundTrackSurface1, + firstIdxAlong / elementsAlongTrackSurface1]) + endProportions2.append([(firstIdxAround2 + 1) / elementsAroundTrackSurface2, + firstIdxAlong / elementsAlongTrackSurface2]) + nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, o1_x, o1_d1, o1_d2, None, o1_NodeId, None, endPoints1_x, endPoints1_d1, endPoints1_d2, None, endNode1_Id, endDerivativesMap, - elementsCountRadial=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups) + elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, + tracksurface=trackSurfaceUreter1, startProportions=startProportions1, endProportions=endProportions1, + rescaleEndDerivatives=True) nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, o2_x, o2_d1, o2_d2, None, o2_NodeId, None, endPoints2_x, endPoints2_d1, endPoints2_d2, None, endNode2_Id, endDerivativesMap, - elementsCountRadial=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups) + elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, + tracksurface=trackSurfaceUreter2, startProportions=startProportions2, endProportions=endProportions2, + rescaleEndDerivatives=True) - # Store elements to be deleted later from bladder mesh + # Delete elements under annulus mesh element_identifiers = [] - for n3 in range(2): - elementIdxUnderOstium1 = ostiumElementPositionDown * elementsCountAround + ostiumElementPositionAround + \ - n3 * elementsCountAround + 1 - element_identifiers.append(elementIdxUnderOstium1) - element_identifiers.append(elementIdxUnderOstium1 + 1) - elementIdxUnderOstium2 = ostiumElementPositionDown * elementsCountAround + elementsCountAround - \ - ostiumElementPositionAround + n3 * elementsCountAround - element_identifiers.append(elementIdxUnderOstium2) - element_identifiers.append(elementIdxUnderOstium2 - 1) - - # Delete elements under annulus mesh - mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers) + elementToDeleteStartIdx1 = elementsCountThroughWall * elementsCountAround * (ureterElementPositionDown - (0 if ureter1Position.xi2 > 0.5 else 1)) \ + + ureterElementPositionAround + (1 if ureter1Position.xi1 > 0.5 else 0) + elementToDeleteStartIdx2 = elementsCountThroughWall * elementsCountAround * (ureterElementPositionDown - (0 if ureter1Position.xi2 > 0.5 else 1)) \ + + elementsCountAround - ureterElementPositionAround + (-1 if ureter1Position.xi1 > 0.5 else 0) + elementToDeleteStartIdxList = [elementToDeleteStartIdx1, elementToDeleteStartIdx2] + for i in range(2): + elementToDeleteStart = elementToDeleteStartIdxList[i] + elementsToDelete = [elementToDeleteStart, elementToDeleteStart + 1, + elementToDeleteStart + elementsCountAround, + elementToDeleteStart + 1 + elementsCountAround] + element_identifiers += elementsToDelete + + mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers) + return diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index cef25e71..6118a7df 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -27,6 +27,69 @@ class MeshType_3d_colon1(Scaffold_base): ''' centralPathDefaultScaffoldPackages = { + 'Cattle 1' : ScaffoldPackage(MeshType_1d_path1, { + 'scaffoldSettings' : { + 'Coordinate dimensions' : 3, + 'D2 derivatives': True, + 'Length' : 1.0, + 'Number of elements' : 52 + }, + 'meshEdits' : exnodeStringFromNodeValues( + [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ + [ [ -245.3, 444.6, -49.1 ], [ -267.7, -53.1, -20.2 ], [ 0.0, 0.0, 15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -380.3, 484.8, -45.0 ], [ 24.5, 102.7, 15.7 ], [ 0.0, 0.0, 15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -298.1, 510.4, -36.8 ], [ 73.6, 9.9, -16.4 ], [ 0.0, 0.0, 15.0 ], [ -10.8, -2.7, -6.5 ] ], + [ [ -213.1, 527.9, -22.5 ], [ -1.0, -10.8, 125.6 ], [ -22.6, -5.6, 1.4 ], [ -5.3, -12.9, -22.2 ] ], + [ [ -315.5, 570.2, 18.9 ], [ -107.9, 9.3, 21.9 ], [ -3.9, -18.1, -31.1 ], [ 14.1, 6.3, -14.4 ] ], + [ [ -417.4, 555.0, 14.6 ], [ -83.0, -41.3, -0.8 ], [ 6.1, 1.7, -31.0 ], [ 1.4, 10.6, 7.1 ] ], + [ [ -497.3, 488.9, 13.6 ], [ -44.6, -81.6, 10.0 ], [ -1.7, -2.6, -21.1 ], [ -3.0, -0.8, 8.0 ] ], + [ [ -527.0, 392.5, 2.7 ], [ 47.4, -82.0, -7.9 ], [ 0.0, 0.0, -15.0 ], [ 1.6, 2.9, 1.1 ] ], + [ [ -461.2, 345.9, -0.8 ], [ 56.9, -44.5, 2.4 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -415.6, 293.8, 3.9 ], [ 93.2, -62.6, 3.1 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -232.2, 264.9, 0.2 ], [ 140.1, 58.2, -1.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -168.4, 357.2, 1.3 ], [ 10.1, 78.6, -3.2 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -185.3, 419.1, -0.7 ], [ -45.1, 57.1, -0.9 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -253.2, 466.7, -0.3 ], [ -63.4, 24.7, 0.2 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -323.8, 482.5, 0.1 ], [ -68.2, 2.9, -1.2 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -387.5, 485.4, -0.2 ], [ -44.2, -17.1, -1.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -435.6, 433.5, 3.3 ], [ 3.4, -109.5, 1.4 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -370.6, 376.3, -1.1 ], [ 66.9, -29.2, -0.9 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -313.0, 357.9, -0.1 ], [ 40.0, -33.5, 9.6 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -259.2, 340.7, 2.1 ], [ 48.9, 6.4, 1.4 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -246.5, 380.3, -0.8 ], [ -29.7, 33.6, -0.7 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -297.3, 387.1, 0.6 ], [ -59.7, 12.6, -0.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -340.2, 415.6, -1.0 ], [ -86.2, 28.9, -2.9 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -398.3, 443.1, -0.1 ], [ 10.6, 82.1, -2.6 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -329.8, 449.1, -2.1 ], [ 53.2, 14.0, -0.5 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -251.3, 425.9, -0.3 ], [ 43.9, -19.3, 0.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -209.1, 390.6, 0.0 ], [ 26.0, -38.8, 0.9 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -207.8, 350.8, 1.4 ], [ -9.4, -43.6, 1.8 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -245.8, 299.4, 7.6 ], [ -70.3, -36.0, 1.4 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -345.3, 304.1, 3.1 ], [ -100.2, 27.9, -1.9 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -418.4, 361.1, -0.2 ], [ -57.8, 55.8, -1.7 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -479.2, 415.6, 2.2 ], [ -8.8, 73.1, -1.6 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -439.6, 495.7, -2.1 ], [ 61.1, 57.1, -1.3 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -361.6, 522.6, -3.0 ], [ 78.6, 9.9, 0.2 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -270.1, 506.5, -3.8 ], [ 103.6, -33.3, 1.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -148.9, 441.4, -2.1 ], [ 79.7, -91.5, 2.8 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -130.9, 313.3, 4.0 ], [ -4.0, -107.2, 3.1 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -183.9, 251.0, 3.8 ], [ -65.5, -60.2, 3.6 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -280.3, 213.0, 3.4 ], [ -165.1, -18.6, 0.1 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -400.8, 247.5, 6.8 ], [ -127.1, 36.8, 1.3 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -530.5, 290.7, 5.2 ], [ -89.0, 86.5, 0.3 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -568.8, 392.3, 6.9 ], [ -77.4, 67.7, -5.5 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -511.2, 535.1, 2.5 ], [ 86.2, 111.4, -1.0 ], [ 0.0, 0.0, -15.0 ], [ 0.0, 0.0, 0.0 ] ], + [ [ -405.0, 601.7, 6.4 ], [ 143.6, 52.2, 2.6 ], [ 0.0, 0.0, -15.0 ], [ 3.2, -3.1, -1.4 ] ], + [ [ -238.8, 615.9, 16.6 ], [ 63.3, -9.1, 19.1 ], [ 8.0, -7.7, -18.4 ], [ -11.9, 7.9, 0.2 ] ], + [ [ -146.2, 605.9, 36.5 ], [ 49.3, -9.9, -50.6 ], [ -23.9, -3.7, -21.6 ], [ -9.9, 8.8, 14.5 ] ], + [ [ -218.4, 585.3, -2.0 ], [ -124.0, 0.4, -37.5 ], [ -9.2, 13.6, 21.7 ], [ 2.9, -0.2, 23.1 ] ], + [ [ -376.3, 579.6, -40.8 ], [ -189.2, -50.7, -8.8 ], [ -5.0, 3.3, 23.8 ], [ 9.6, -7.2, 10.6 ] ], + [ [ -557.9, 493.9, -24.9 ], [ -30.3, 24.1, 152.8 ], [ 27.2, 29.2, 3.0 ], [ 3.8, -10.1, -29.6 ] ], + [ [ -484.8, 594.4, 0.7 ], [ 132.7, 97.0, 3.5 ], [ 12.4, -4.5, -32.3 ], [ -12.3, -22.3, -22.5 ] ], + [ [ -318.1, 641.9, -8.5 ], [ 166.7, 17.6, 5.5 ], [ 3.0, -13.0, -39.4 ], [ -8.3, -3.3, -0.9 ] ], + [ [ -158.3, 634.7, -1.9 ], [ 176.5, -14.0, 10.8 ], [ -4.3, -11.5, -34.6 ], [ -3.1, 2.7, 5.3 ] ], + [ [ 32.7, 611.7, 13.6 ], [ 205.5, -32.2, 20.0 ], [ -2.4, -7.3, -28.7 ], [ 6.9, 5.6, 6.4 ] ] ] ) + } ), 'Human 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, @@ -103,47 +166,50 @@ class MeshType_3d_colon1(Scaffold_base): 'Coordinate dimensions' : 3, 'D2 derivatives': True, 'Length' : 1.0, - 'Number of elements' : 36 + 'Number of elements' : 39 }, 'meshEdits' : exnodeStringFromNodeValues( [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2 ], [ - [ [ -7.2, 83.3, -20.7 ], [ -65.2, -8.1, 7.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -68.5, 52.8, -9.6 ], [ -40.1, -36.1, 10.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -97.4, -26.3, 5.7 ], [ 18.0, -93.2, 13.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -56.8, -90.5, 14.1 ], [ 65.5, -41.4, 7.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 48.9, -100.8, 24.0 ], [ 112.2, 40.1, 19.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 114.8, -12.6, 38.7 ], [ 8.2, 96.1, 14.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 60.3, 83.5, 43.7 ], [ -108.7, 54.1, 22.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -41.2, 90.7, 56.3 ], [ -89.0, -32.4, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -107.9, -9.7, 76.6 ], [ 11.1, -94.4, 11.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -57.3, -91.9, 81.3 ], [ 71.2, -31.2, 5.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 51.2, -89.4, 97.2 ], [ 99.1, 55.4, 12.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 91.6, 9.3, 103.6 ], [ 4.7, 51.2, 3.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 61.6, 111.8, 109.6 ], [ -85.2, 46.1, 2.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -54.6, 91.9, 129.4 ], [ -92.7, -55.0, 14.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -109.0, 5.6, 156.9 ], [ 23.6, -108.2, 27.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -59.1, -62.5, 170.8 ], [ 74.0, -20.1, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 23.5, -53.2, 179.7 ], [ 84.6, 47.0, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 62.3, 30.1, 187.5 ], [ -12.8, 58.0, 0.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 22.4, 45.2, 181.1 ], [ -23.6, -34.5, -7.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -1.9, 4.9, 180.5 ], [ -41.3, -30.9, 7.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -45.1, -12.6, 194.4 ], [ -40.5, -4.6, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -71.7, -2.2, 197.2 ], [ -25.2, 35.8, -6.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -65.8, 42.1, 182.3 ], [ 26.6, 37.6, -15.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -14.1, 81.2, 163.5 ], [ 41.0, 10.3, -9.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 61.7, 86.1, 156.4 ], [ 77.9, -40.7, 8.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 92.9, 20.5, 150.3 ], [ 0.0, -73.3, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 48.9, -65.0, 142.8 ], [ -82.8, -80.0, -1.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -54.3, -90.8, 134.0 ], [ -60.1, 26.4, -8.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -89.9, 11.2, 115.0 ], [ 34.9, 125.1, -27.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -17.4, 74.2, 91.1 ], [ 78.8, 19.1, -15.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 43.4, 50.2, 73.3 ], [ 30.2, -36.0, -9.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 62.4, -5.1, 63.5 ], [ 10.9, -54.2, -2.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 32.7, -51.7, 56.1 ], [ -38.6, -29.8, -8.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -38.1, -28.6, 46.8 ], [ -66.0, 83.3, -12.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ -5.4, 32.0, 26.0 ], [ 48.1, 17.6, -21.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 146.6, 101.2, -41.2 ], [ 63.3, 35.3, -31.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], - [ [ 312.3, 199.5, -123.4 ], [ 39.7, 24.3, -20.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ] ] ) + [ [ -7.2, 83.3, -20.7 ], [ -65.2, -8.1, 7.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -68.5, 52.8, -9.6 ], [ -40.1, -36.1, 10.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -97.4, -26.3, 5.7 ], [ 18.0, -93.2, 13.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -56.8, -90.5, 14.1 ], [ 65.5, -41.4, 7.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 48.9, -100.8, 24.0 ], [ 112.2, 40.1, 19.0 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 114.8, -12.6, 38.7 ], [ 8.2, 96.1, 14.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 60.3, 83.5, 43.7 ], [ -108.7, 54.1, 22.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -41.2, 90.7, 56.3 ], [ -89.0, -32.4, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -107.9, -9.7, 76.6 ], [ 11.1, -94.4, 11.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -57.3, -91.9, 81.3 ], [ 71.2, -31.2, 5.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 51.2, -89.4, 97.2 ], [ 99.1, 55.4, 12.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 91.6, 9.3, 103.6 ], [ 4.7, 51.2, 3.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 61.6, 111.8, 109.6 ], [ -85.2, 46.1, 2.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -54.6, 91.9, 129.4 ], [ -92.7, -55.0, 14.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -109.0, 5.6, 156.9 ], [ 23.6, -108.2, 27.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -59.1, -62.5, 170.8 ], [ 74.0, -20.1, 14.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 23.5, -53.2, 179.7 ], [ 84.6, 47.0, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 62.3, 30.1, 187.5 ], [ -12.8, 58.0, 0.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 22.4, 45.2, 181.1 ], [ -23.6, -34.5, -7.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -1.9, 4.9, 180.5 ], [ -41.3, -30.9, 7.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -45.1, -12.6, 194.4 ], [ -40.5, -4.6, 6.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -71.7, -2.2, 197.2 ], [ -25.2, 35.8, -6.8 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -65.8, 42.1, 182.3 ], [ 26.6, 37.6, -15.6 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -14.1, 81.2, 163.5 ], [ 41.0, 10.3, -9.5 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 61.7, 86.1, 156.4 ], [ 77.9, -40.7, 8.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 92.9, 20.5, 150.3 ], [ 0.0, -73.3, -5.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 48.9, -65.0, 142.8 ], [ -82.8, -80.0, -1.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -54.3, -90.8, 134.0 ], [ -60.1, 26.4, -8.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -89.9, 11.2, 115.0 ], [ 34.9, 125.1, -27.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -17.4, 74.2, 91.1 ], [ 78.8, 19.1, -15.4 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 43.4, 50.2, 73.3 ], [ 30.2, -36.0, -9.9 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 62.4, -5.1, 63.5 ], [ 10.9, -54.2, -2.7 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 32.7, -51.7, 56.1 ], [ -38.6, -29.8, -8.1 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -38.1, -28.6, 46.8 ], [ -62.5, 82.6, -19.2 ], [ 4.0, 6.8, 13.2 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 5.7, 40.4, 22.4 ], [ 144.8, 18.6, -20.5 ], [ 4.3, -13.3, 12.3 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 53.0, -14.7, -4.1 ], [ -6.0, -25.7, -46.7 ], [ -13.5, -19.4, 4.0 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 24.8, -0.4, -48.8 ], [ -13.4, 23.9, -30.6 ], [ -17.4, -15.1, 6.4 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -20.9, 15.3, -77.9 ], [ -51.2, -30.6, 21.1 ], [ 7.1, -7.7, 12.2 ], [ 0.0, 0.0, 0.5 ] ], + [ [ -47.6, 33.9, -112.2 ], [ 32.6, 30.7, -27.8 ], [ -12.8, -0.5, -1.6 ], [ 0.0, 0.0, 0.5 ] ], + [ [ 19.6, 96.0, -167.5 ], [ 19.9, 19.1, -18.4 ], [ -12.6, 1.3, -8.1 ], [ 0.0, 0.0, 0.5 ] ] ] ) } ), 'Pig 2': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { @@ -169,6 +235,7 @@ def getName(): def getParameterSetNames(): return [ 'Default', + 'Cattle 1', 'Human 1', 'Human 2', 'Mouse 1', @@ -180,6 +247,8 @@ def getParameterSetNames(): def getDefaultOptions(cls, parameterSetName='Default'): if 'Human 2' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 2'] + elif 'Cattle 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Cattle 1'] elif 'Mouse 1' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 1'] elif 'Mouse 2' in parameterSetName: @@ -190,7 +259,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 2'] else: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 1'] - if 'Mouse' in parameterSetName: + if 'Cattle' in parameterSetName: + segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Cattle 1') + elif 'Mouse' in parameterSetName: segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Mouse 1') elif 'Pig' in parameterSetName: segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Pig 1') @@ -217,7 +288,20 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Refine number of elements along' : 1, 'Refine number of elements through wall' : 1 } - if 'Human 2' in parameterSetName: + if 'Cattle 1' in parameterSetName: + options['Number of segments'] = 40 + options['Proximal length'] = 880.0 + options['Transverse length'] = 3500.0 + options['Distal length'] = 1650.0 + options['Proximal inner radius'] = 30.0 + options['Proximal tenia coli width'] = 3.0 + options['Proximal-transverse inner radius'] = 15.0 + options['Proximal-transverse tenia coli width'] = 3.0 + options['Transverse-distal inner radius'] = 9.0 + options['Transverse-distal tenia coli width'] = 3.0 + options['Distal inner radius'] = 10.5 + options['Distal tenia coli width'] = 3.0 + elif 'Human 2' in parameterSetName: options['Proximal length'] = 180.0 options['Transverse length'] = 620.0 options['Distal length'] = 700.0 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index 9dbae5f9..04c05ebe 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -43,6 +43,7 @@ def getName(): def getParameterSetNames(): return [ 'Default', + 'Cattle 1', 'Human 1', 'Mouse 1', 'Pig 1'] @@ -78,7 +79,19 @@ def getDefaultOptions(parameterSetName='Default'): 'Refine number of elements along segment' : 1, 'Refine number of elements through wall' : 1 } - if 'Mouse' in parameterSetName: + if 'Cattle' in parameterSetName: + options['Start inner radius'] = 10.5 + options['End inner radius'] = 10.5 + options['Corner inner radius factor'] = 0.0 + options['Haustrum inner radius factor'] = 0.0 + options['Segment length end derivative factor'] = 0.0 + options['Segment length mid derivative factor'] = 0.0 + options['Number of tenia coli'] = 1 + options['Start tenia coli width'] = 3.0 + options['End tenia coli width'] = 3.0 + options['Tenia coli thickness'] = 0.0 + options['Wall thickness'] = 3.02 + elif 'Mouse' in parameterSetName: options['Start inner radius'] = 0.94 options['End inner radius'] = 0.94 options['Corner inner radius factor'] = 0.0 diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 360f0ae6..f06814c2 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -14,26 +14,71 @@ from scaffoldmaker.utils import vector -def derivativeSignsToExpressionTerms(valueLabels, signs): +def derivativeSignsToExpressionTerms(valueLabels, signs, scaleFactorIdx = None): ''' - Return remap expression terms for summing derivative[i]*sign[i] + Return remap expression terms for summing derivative[i]*sign[i]*scaleFactor :param valueLabels: List of node value labels to possibly include. :param signs: List of 1 (no scaling), -1 (scale by scale factor 1) or 0 (no term). + :param scaleFactorIdx: Optional index of local scale factor to scale all non-zero terms. Default None means no extra scaling. ''' expressionTerms = [] for i in range(len(valueLabels)): if signs[i] is 1: - expressionTerms.append( ( valueLabels[i], [] ) ) + expressionTerms.append((valueLabels[i], ([ scaleFactorIdx ] if scaleFactorIdx else []))) elif signs[i] is -1: - expressionTerms.append( ( valueLabels[i], [1] ) ) + expressionTerms.append((valueLabels[i], ([1, scaleFactorIdx] if scaleFactorIdx else [1]))) return expressionTerms +def getMappedD1D2(gds, derivativesMaps): + ''' + Get vector combinations of d1In, d2In, d3In indicated by derivativesMap. + :param gds: List of global d1, d2 and optionally d3. + :param derivativesMaps: List over d1, d2, d3, and optionally d1b (for + different d1 exiting global node) of list of 3 weights of gds, + each limited to -1.0, 0.0, or -1.0. + :return: Effective d1, d2. Where d1 is around, d2 is radial. + ''' + dslimit = len(gds) + if not (derivativesMaps and derivativesMaps[0]): + d1 = gds[0] + else: + derivativesMap = derivativesMaps[0] + d1 = [ 0.0, 0.0, 0.0 ] + for ds in range(dslimit): + if derivativesMap[ds] != 0.0: + for c in range(3): + d1[c] += derivativesMap[ds]*gds[ds][c] + if len(derivativesMaps) > 3: + # average with d1 map for other side + derivativesMap = derivativesMaps[3] + d1 = [ 0.5*d for d in d1 ] + if not derivativesMap: + for c in range(3): + d1[c] += 0.5*gds[0][c] + else: + for ds in range(dslimit): + if derivativesMap[ds] != 0.0: + for c in range(3): + d1[c] += 0.5*derivativesMap[ds]*gds[ds][c] + if not (derivativesMaps and derivativesMaps[1]): + d2 = gds[1] + else: + derivativesMap = derivativesMaps[1] + d2 = [ 0.0, 0.0, 0.0 ] + for ds in range(dslimit): + if derivativesMap[ds] != 0.0: + for c in range(3): + d2[c] += derivativesMap[ds]*gds[ds][c] + return d1, d2 + + def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, startPointsx, startPointsd1, startPointsd2, startPointsd3, startNodeId, startDerivativesMap, endPointsx, endPointsd1, endPointsd2, endPointsd3, endNodeId, endDerivativesMap, forceStartLinearXi3 = False, forceMidLinearXi3 = False, forceEndLinearXi3 = False, maxStartThickness = None, maxEndThickness = None, useCrossDerivatives = False, - elementsCountRadial = 1, meshGroups = None, tracksurface = None, startProportions = None, endProportions = None): + elementsCountRadial = 1, meshGroups = None, tracksurface = None, startProportions = None, endProportions = None, + rescaleStartDerivatives = False, rescaleEndDerivatives = False, sampleBlend = 0.0): """ Create an annulus mesh from a loop of start points/nodes with specified derivative mappings to a loop of end points/nodes with specified derivative mappings. @@ -72,6 +117,14 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, around as for startPoints. Values only given for tracksurface for outer layer (xi3 == 1). :param endProportions: Proportion around and along of endPoints on track surface. These vary with nodes around as for endPoints. Values only given for tracksurface for outer layer (xi3 == 1). + :param rescaleStartDerivatives, rescaleEndDerivatives: Optional flags to compute and multiply additional scale factors + on start, end or both radial derivatives to fit arc length, needed if derivatives are of the wrong scale for the radial + distances and the chosen elementsCountRadial. If either is True, derivatives and sampled radial nodes are spaced for a + gradual change of derivative from that at the other end. If both are True, scaling is set to give even sampling and arc + length derivatives. + :param sampleBlend: Real value varying from 0.0 to 1.0 controlling weighting of start and end + derivatives when interpolating extra points in-between, where 0.0 = sample with equal end derivatives, + and 1.0 = proportional to current magnitudes, interpolated in between. :return: Final values of nextNodeIdentifier, nextElementIdentifier """ assert (elementsCountRadial >= 1), 'createAnnulusMesh3d: Invalid number of radial elements' @@ -134,178 +187,147 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, pd2[n3] = [ startPointsd2[n3], endPointsd2[n3] ] pd3[n3] = [ startPointsd3[n3] if (startPointsd3 is not None) else None, \ endPointsd3[n3] if (endPointsd3 is not None) else None ] - if elementsCountRadial > 1: - # add in-between points - startPointsd = [ startPointsd1, startPointsd2, startPointsd3 ] - startPointsdslimit = 2 if (startPointsd3 is None) else 3 - endPointsd = [ endPointsd1, endPointsd2, endPointsd3 ] - endPointsdslimit = 2 if (endPointsd3 is None) else 3 - for n3 in range(2): - for n2 in range(1, elementsCountRadial): - px [n3].insert(n2, [ None ]*nodesCountAround) - pd1[n3].insert(n2, [ None ]*nodesCountAround) - pd2[n3].insert(n2, [ None ]*nodesCountAround) - pd3[n3].insert(n2, None if midLinearXi3 else [ None ]*nodesCountAround) - # compute on outside / n3 = 1, then map to inside using thickness - thicknesses = [] - thicknesses.append([ vector.magnitude([ (startPointsx[1][n1][c] - startPointsx[0][n1][c]) for c in range(3) ]) for n1 in range(nodesCountAround) ]) - if maxStartThickness: - for n1 in range(nodesCountAround): - thicknesses[0][n1] = min(thicknesses[0][n1], maxStartThickness) + + if rescaleStartDerivatives: + scaleFactorMapStart = [ [] for n2 in range(2) ] + if rescaleEndDerivatives: + scaleFactorMapEnd = [ [] for n2 in range(2) ] + + # following code adds in-between points, but also handles rescaling for 1 radial element + for n3 in range(2): for n2 in range(1, elementsCountRadial): - thicknesses.append([ None ]*nodesCountAround) - thicknesses.append([ vector.magnitude([ (endPointsx[1][n1][c] - endPointsx[0][n1][c]) for c in range(3) ]) for n1 in range(nodesCountAround) ]) - if maxEndThickness: - for n1 in range(nodesCountAround): - thicknesses[-1][n1] = min(thicknesses[-1][n1], maxEndThickness) - n3 == 1 + px [n3].insert(n2, [ None ]*nodesCountAround) + pd1[n3].insert(n2, [ None ]*nodesCountAround) + pd2[n3].insert(n2, [ None ]*nodesCountAround) + pd3[n3].insert(n2, None if midLinearXi3 else [ None ]*nodesCountAround) + # compute on outside / n3 = 1, then map to inside using thickness + thicknesses = [] + thicknesses.append([ vector.magnitude([ (startPointsx[1][n1][c] - startPointsx[0][n1][c]) for c in range(3) ]) for n1 in range(nodesCountAround) ]) + if maxStartThickness: for n1 in range(nodesCountAround): - ax = startPointsx [n3][n1] - if (startDerivativesMap is None) or (startDerivativesMap[n3][n1][0] is None): - ad1 = startPointsd1[n3][n1] - else: - derivativesMap = startDerivativesMap[n3][n1][0] - ad1 = [ 0.0, 0.0, 0.0 ] - for ds in range(startPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - ad1[c] += derivativesMap[ds]*startPointsd[ds][n3][n1][c] - if len(startDerivativesMap[n3][n1]) > 3: - # average with d1 map for other side - derivativesMap = startDerivativesMap[n3][n1][3] - ad1 = [ 0.5*d for d in ad1 ] - if not derivativesMap: - for c in range(3): - ad1[c] += 0.5*startPointsd[0][n3][n1][c] - else: - for ds in range(startPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - ad1[c] += 0.5*derivativesMap[ds]*startPointsd[ds][n3][n1][c] - if (startDerivativesMap is None) or (startDerivativesMap[n3][n1][1] is None): - ad2 = startPointsd2[n3][n1] - else: - derivativesMap = startDerivativesMap[n3][n1][1] - ad2 = [ 0.0, 0.0, 0.0 ] - for ds in range(startPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - ad2[c] += derivativesMap[ds]*startPointsd[ds][n3][n1][c] - - bx = endPointsx [n3][n1] - if (endDerivativesMap is None) or (endDerivativesMap[n3][n1][0] is None): - bd1 = endPointsd1[n3][n1] - else: - derivativesMap = endDerivativesMap[n3][n1][0] - bd1 = [ 0.0, 0.0, 0.0 ] - for ds in range(endPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - bd1[c] += derivativesMap[ds]*endPointsd[ds][n3][n1][c] - if len(endDerivativesMap[n3][n1]) > 3: - # average with d1 map for other side - derivativesMap = endDerivativesMap[n3][n1][3] - bd1 = [ 0.5*d for d in bd1 ] - if not derivativesMap: - for c in range(3): - bd1[c] += 0.5*endPointsd[0][n3][n1][c] - else: - for ds in range(endPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - bd1[c] += 0.5*derivativesMap[ds]*endPointsd[ds][n3][n1][c] - if (endDerivativesMap is None) or (endDerivativesMap[n3][n1][1] is None): - bd2 = endPointsd2[n3][n1] - else: - derivativesMap = endDerivativesMap[n3][n1][1] - bd2 = [ 0.0, 0.0, 0.0 ] - for ds in range(endPointsdslimit): - if derivativesMap[ds] != 0.0: - for c in range(3): - bd2[c] += derivativesMap[ds]*endPointsd[ds][n3][n1][c] - - # scaling end derivatives to arc length gives even curvature along the curve - if tracksurface: - mx, md2, md1 = \ - tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], - endProportions[n1][0], endProportions[n1][1], - elementsCountRadial, derivativeStart=ad2, derivativeEnd=bd2)[0:3] - else: - arcLength = interp.computeCubicHermiteArcLength(ax, ad2, bx, bd2, rescaleDerivatives = False) - scaledDerivatives = [ vector.setMagnitude(d2, arcLength) for d2 in [ ad2, bd2 ]] - mx, md2, me, mxi = interp.sampleCubicHermiteCurvesSmooth([ ax, bx ], scaledDerivatives, elementsCountRadial, - derivativeMagnitudeStart = vector.magnitude(ad2), - derivativeMagnitudeEnd = vector.magnitude(bd2))[0:4] - md1 = interp.interpolateSampleLinear([ ad1, bd1 ], me, mxi) - thi = interp.interpolateSampleLinear([ thicknesses[0][n1], thicknesses[-1][n1] ], me, mxi) + thicknesses[0][n1] = min(thicknesses[0][n1], maxStartThickness) + for n2 in range(1, elementsCountRadial): + thicknesses.append([ None ]*nodesCountAround) + thicknesses.append([ vector.magnitude([ (endPointsx[1][n1][c] - endPointsx[0][n1][c]) for c in range(3) ]) for n1 in range(nodesCountAround) ]) + if maxEndThickness: + for n1 in range(nodesCountAround): + thicknesses[-1][n1] = min(thicknesses[-1][n1], maxEndThickness) + n3 == 1 + for n1 in range(nodesCountAround): + ax = startPointsx[n3][n1] + ad1, ad2 = getMappedD1D2([ startPointsd1[n3][n1], startPointsd2[n3][n1] ] + ([ startPointsd3[n3][n1] ] if startPointsd3 else []), + startDerivativesMap[n3][n1] if startDerivativesMap else None) + bx = endPointsx[n3][n1] + bd1, bd2 = getMappedD1D2([ endPointsd1[n3][n1], endPointsd2[n3][n1] ] + ([ endPointsd3[n3][n1] ] if endPointsd3 else []), + endDerivativesMap[n3][n1] if endDerivativesMap else None) - #md2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, fixStartDerivative = True, fixEndDerivative = True) + # sample between start and end points and derivatives + # scaling end derivatives to arc length gives even curvature along the curve + aMag = vector.magnitude(ad2) + bMag = vector.magnitude(bd2) + ad2Scaled = vector.setMagnitude(ad2, 0.5*((1.0 + sampleBlend)*aMag + (1.0 - sampleBlend)*bMag)) + bd2Scaled = vector.setMagnitude(bd2, 0.5*((1.0 + sampleBlend)*bMag + (1.0 - sampleBlend)*aMag)) + scaling = interp.computeCubicHermiteDerivativeScaling(ax, ad2Scaled, bx, bd2Scaled) + ad2Scaled = [ d*scaling for d in ad2Scaled ] + bd2Scaled = [ d*scaling for d in bd2Scaled ] + derivativeMagnitudeStart = None if rescaleStartDerivatives else vector.magnitude(ad2) + derivativeMagnitudeEnd = None if rescaleEndDerivatives else vector.magnitude(bd2) + if tracksurface: + mx, md2, md1, md3, mProportions = \ + tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], + endProportions[n1][0], endProportions[n1][1], + elementsCountRadial, + derivativeStart=[ d/elementsCountRadial for d in ad2Scaled ], + derivativeEnd=[ d/elementsCountRadial for d in bd2Scaled ]) + mx, md2, md1 = tracksurface.resampleHermiteCurvePointsSmooth(mx, md2, md1, md3, mProportions, derivativeMagnitudeStart, derivativeMagnitudeEnd)[0:3] + # interpolate thicknesses using xi calculated from radial arclength distances to points + arcLengthInsideToRadialPoint = [ 0.0 ] + [ interp.getCubicHermiteArcLength(mx[n2], md2[n2], mx[n2 + 1], md2[n2 + 1]) for n2 in range(elementsCountRadial) ] + arclengthInsideToOutside = sum(arcLengthInsideToRadialPoint) + thi = [] + for n2 in range(elementsCountRadial + 1): + xi = arcLengthInsideToRadialPoint[n2 - 1]/arclengthInsideToOutside + thi.append(thicknesses[0][n1]*xi + thicknesses[-1][n1]*(1.0 - xi)) + else: + mx, md2, me, mxi = interp.sampleCubicHermiteCurvesSmooth([ ax, bx ], [ ad2Scaled, bd2Scaled ], elementsCountRadial, + derivativeMagnitudeStart, derivativeMagnitudeEnd)[0:4] + md1 = interp.interpolateSampleLinear([ ad1, bd1 ], me, mxi) + thi = interp.interpolateSampleLinear([ thicknesses[0][n1], thicknesses[-1][n1] ], me, mxi) - for n2 in range(1, elementsCountRadial): - px [n3][n2][n1] = mx [n2] - pd1[n3][n2][n1] = md1[n2] - pd2[n3][n2][n1] = md2[n2] - if not tracksurface: - thicknesses[n2][n1] = thi[n2] + # set scalefactors if rescaling, make same on inside for now + if rescaleStartDerivatives: + scaleFactor = vector.magnitude(md2[0])/vector.magnitude(ad2) + scaleFactorMapStart[n3].append(scaleFactor) + if rescaleEndDerivatives: + scaleFactor = vector.magnitude(md2[-1])/vector.magnitude(bd2) + scaleFactorMapEnd[n3].append(scaleFactor) - # Interpolate thicknesses using xi calculated using arclength distances between points - if tracksurface: - for n1 in range(nodesCountAround): - arclengthInsideToOutside = 0.0 - arcLengthInsideToRadialPoint = [] - for n2 in range(elementsCountRadial): - arclengthInsideToOutside += interp.computeCubicHermiteArcLength(px[n3][n2][n1], pd2[n3][n2][n1], px[n3][n2+1][n1], pd2[n3][n2+1][n1], False) - if n2 < elementsCountRadial: - arcLengthInsideToRadialPoint.append(arclengthInsideToOutside) - for n2 in range(elementsCountRadial - 1): - xi = arcLengthInsideToRadialPoint[n2]/arclengthInsideToOutside - thicknesses[n2+1][n1] = thicknesses[0][n1]*xi + thicknesses[-1][n1]*(1-xi) - - # now get inner positions from normal and thickness, derivatives from curvature for n2 in range(1, elementsCountRadial): - # first smooth derivative 1 around outer loop - pd1[1][n2] = interp.smoothCubicHermiteDerivativesLoop(px[1][n2], pd1[1][n2], magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) - - for n1 in range(nodesCountAround): - normal = vector.normalise(vector.crossproduct3(pd1[1][n2][n1], pd2[1][n2][n1])) - thickness = thicknesses[n2][n1] - d3 = [ d*thickness for d in normal ] - px [0][n2][n1] = [ (px [1][n2][n1][c] - d3[c]) for c in range(3) ] - # calculate inner d1 from curvature around - n1m = n1 - 1 - n1p = (n1 + 1)%nodesCountAround - curvature = 0.5*( - interp.getCubicHermiteCurvature(px[1][n2][n1m], pd1[1][n2][n1m], px[1][n2][n1 ], pd1[1][n2][n1 ], normal, 1.0) + - interp.getCubicHermiteCurvature(px[1][n2][n1 ], pd1[1][n2][n1 ], px[1][n2][n1p], pd1[1][n2][n1p], normal, 0.0)) - factor = 1.0 + curvature*thickness - pd1[0][n2][n1] = [ factor*d for d in pd1[1][n2][n1] ] - # calculate inner d2 from curvature radially - n2m = n2 - 1 - n2p = n2 + 1 - curvature = 0.5*( - interp.getCubicHermiteCurvature(px[1][n2m][n1], pd2[1][n2m][n1], px[1][n2 ][n1], pd2[1][n2 ][n1], normal, 1.0) + - interp.getCubicHermiteCurvature(px[1][n2 ][n1], pd2[1][n2 ][n1], px[1][n2p][n1], pd2[1][n2p][n1], normal, 0.0)) - factor = 1.0 + curvature*thickness - pd2[0][n2][n1] = [ factor*d for d in pd2[1][n2][n1] ] - d2Scaled = [factor*d for d in pd2[1][n2][n1]] - if vector.dotproduct(vector.normalise(pd2[1][n2][n1]), vector.normalise(d2Scaled)) == -1: - pd2[0][n2][n1] = [-factor * d for d in pd2[1][n2][n1]] - if not midLinearXi3: - pd3[0][n2][n1] = pd3[1][n2][n1] = d3 - - # smooth derivative 1 around inner loop - pd1[0][n2] = interp.smoothCubicHermiteDerivativesLoop(px[0][n2], pd1[0][n2], magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) - - for n3 in range(0, 1): # was (0, nodesCountWall) - # smooth derivative 2 radially/along annulus - for n1 in range(nodesCountAround): - sd2 = interp.smoothCubicHermiteDerivativesLine( - [ px [n3][n2][n1] for n2 in range(elementsCountRadial + 1) ], - [ pd2[n3][n2][n1] for n2 in range(elementsCountRadial + 1) ], - fixAllDirections = True, fixStartDerivative = True, fixEndDerivative = True, - magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) - for n2 in range(elementsCountRadial + 1): - pd2[n3][n2][n1] = sd2[n2] + px [n3][n2][n1] = mx [n2] + pd1[n3][n2][n1] = md1[n2] + pd2[n3][n2][n1] = md2[n2] + thicknesses[n2][n1] = thi[n2] + + # now get inner positions from normal and thickness, derivatives from curvature + for n2 in range(1, elementsCountRadial): + # first smooth derivative 1 around outer loop + pd1[1][n2] = interp.smoothCubicHermiteDerivativesLoop(px[1][n2], pd1[1][n2], magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + + for n1 in range(nodesCountAround): + normal = vector.normalise(vector.crossproduct3(pd1[1][n2][n1], pd2[1][n2][n1])) + thickness = thicknesses[n2][n1] + d3 = [ d*thickness for d in normal ] + px [0][n2][n1] = [ (px [1][n2][n1][c] - d3[c]) for c in range(3) ] + # calculate inner d1 from curvature around + n1m = n1 - 1 + n1p = (n1 + 1)%nodesCountAround + curvature = 0.5*( + interp.getCubicHermiteCurvature(px[1][n2][n1m], pd1[1][n2][n1m], px[1][n2][n1 ], pd1[1][n2][n1 ], normal, 1.0) + + interp.getCubicHermiteCurvature(px[1][n2][n1 ], pd1[1][n2][n1 ], px[1][n2][n1p], pd1[1][n2][n1p], normal, 0.0)) + factor = 1.0 + curvature*thickness + pd1[0][n2][n1] = [ factor*d for d in pd1[1][n2][n1] ] + # calculate inner d2 from curvature radially + n2m = n2 - 1 + n2p = n2 + 1 + curvature = 0.5*( + interp.getCubicHermiteCurvature(px[1][n2m][n1], pd2[1][n2m][n1], px[1][n2 ][n1], pd2[1][n2 ][n1], normal, 1.0) + + interp.getCubicHermiteCurvature(px[1][n2 ][n1], pd2[1][n2 ][n1], px[1][n2p][n1], pd2[1][n2p][n1], normal, 0.0)) + factor = 1.0 + curvature*thickness + pd2[0][n2][n1] = [ factor*d for d in pd2[1][n2][n1] ] + d2Scaled = [factor*d for d in pd2[1][n2][n1]] + if vector.dotproduct(vector.normalise(pd2[1][n2][n1]), vector.normalise(d2Scaled)) == -1: + pd2[0][n2][n1] = [-factor * d for d in pd2[1][n2][n1]] + if not midLinearXi3: + pd3[0][n2][n1] = pd3[1][n2][n1] = d3 + + # smooth derivative 1 around inner loop + pd1[0][n2] = interp.smoothCubicHermiteDerivativesLoop(px[0][n2], pd1[0][n2], magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + + for n3 in range(0, 1): # was (0, nodesCountWall) + # smooth derivative 2 radially/along annulus + for n1 in range(nodesCountAround): + mx = [ px [n3][n2][n1] for n2 in range(elementsCountRadial + 1) ] + md2 = [ pd2[n3][n2][n1] for n2 in range(elementsCountRadial + 1) ] + # replace mapped start/end d2 + md2[0] = getMappedD1D2([ startPointsd1[n3][n1], startPointsd2[n3][n1] ] + ([ startPointsd3[n3][n1] ] if startPointsd3 else []), + startDerivativesMap[n3][n1] if startDerivativesMap else None)[1] + md2[-1] = getMappedD1D2([ endPointsd1[n3][n1], endPointsd2[n3][n1] ] + ([ endPointsd3[n3][n1] ] if endPointsd3 else []), + endDerivativesMap[n3][n1] if endDerivativesMap else None)[1] + + sd2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, fixAllDirections = True, + fixStartDerivative = not rescaleStartDerivatives, + fixStartDirection = rescaleStartDerivatives, + fixEndDerivative = not rescaleEndDerivatives, + fixEndDirection = rescaleEndDerivatives, + magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + if rescaleStartDerivatives: + scaleFactor = vector.magnitude(sd2[0]) / vector.magnitude(md2[0]) + scaleFactorMapStart[n3].append(scaleFactor) + if rescaleEndDerivatives: + scaleFactor = vector.magnitude(sd2[-1]) / vector.magnitude(md2[-1]) + scaleFactorMapEnd[n3].append(scaleFactor) + + for n2 in range(1, elementsCountRadial): + pd2[n3][n2][n1] = sd2[n2] ############## # Create nodes @@ -373,43 +395,90 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, eftFactory = tricubichermite if nonlinearXi3 else bicubichermitelinear eftStandard = eftFactory.createEftBasic() elementtemplateStandard.defineField(coordinates, -1, eftStandard) - mapStartDerivatives = (e2 == 0) and (startDerivativesMap is not None) + mapStartDerivatives = (e2 == 0) and (startDerivativesMap or rescaleStartDerivatives) mapStartLinearDerivativeXi3 = nonlinearXi3 and rowLinearXi3[e2] - mapEndDerivatives = (e2 == (elementsCountRadial - 1)) and (endDerivativesMap is not None) + mapEndDerivatives = (e2 == (elementsCountRadial - 1)) and (endDerivativesMap or rescaleEndDerivatives) mapEndLinearDerivativeXi3 = nonlinearXi3 and rowLinearXi3[e2 + 1] mapDerivatives = mapStartDerivatives or mapStartLinearDerivativeXi3 or mapEndDerivatives or mapEndLinearDerivativeXi3 for e1 in range(elementsCountAround): en = (e1 + 1)%elementsCountAround nids = [ nodeId[0][e2][e1], nodeId[0][e2][en], nodeId[0][e2 + 1][e1], nodeId[0][e2 + 1][en], nodeId[1][e2][e1], nodeId[1][e2][en], nodeId[1][e2 + 1][e1], nodeId[1][e2 + 1][en] ] + scaleFactors = [] if mapDerivatives: eft1 = eftFactory.createEftNoCrossDerivatives() - setEftScaleFactorIds(eft1, [1], []) + # work out if scaling by global -1 + scaleMinus1 = mapStartLinearDerivativeXi3 or mapEndLinearDerivativeXi3 + if (not scaleMinus1) and startDerivativesMap: + for n3 in range(2): + for i in range(2): + derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] + for map in derivativesMap: + if map and (-1 in map): + scaleMinus1 = True + break + if (not scaleMinus1) and endDerivativesMap: + for n3 in range(2): + for i in range(2): + derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] + for map in derivativesMap: + if map and (-1 in map): + scaleMinus1 = True + break + # make node scale factors vary fastest by local node varying across lower xi + nodeScaleFactorIds = [] + for n3 in range(2): + if mapStartDerivatives and rescaleStartDerivatives: + for i in range(2): + derivativesMap = (startDerivativesMap[n3][e1][1] if (i == 0) else startDerivativesMap[n3][en][1]) if startDerivativesMap else None + nodeScaleFactorIds.append(getQuadrantID(derivativesMap if derivativesMap else (0, 1, 0))) + if mapEndDerivatives and rescaleEndDerivatives: + for i in range(2): + derivativesMap = (endDerivativesMap[n3][e1][1] if (i == 0) else endDerivativesMap[n3][en][1]) if endDerivativesMap else None + nodeScaleFactorIds.append(getQuadrantID(derivativesMap if derivativesMap else (0, 1, 0))) + setEftScaleFactorIds(eft1, [1] if scaleMinus1 else [], nodeScaleFactorIds) + firstNodeScaleFactorIndex = 2 if scaleMinus1 else 1 + firstStartNodeScaleFactorIndex = firstNodeScaleFactorIndex if (mapStartDerivatives and rescaleStartDerivatives) else None + firstEndNodeScaleFactorIndex = (firstNodeScaleFactorIndex + (2 if firstStartNodeScaleFactorIndex else 0)) if (mapEndDerivatives and rescaleEndDerivatives) else None + layerNodeScaleFactorIndexOffset = 4 if (firstStartNodeScaleFactorIndex and firstEndNodeScaleFactorIndex) else 2 + if scaleMinus1: + scaleFactors.append(-1.0) + for n3 in range(2): + if firstStartNodeScaleFactorIndex: + scaleFactors.append(scaleFactorMapStart[n3][e1]) + scaleFactors.append(scaleFactorMapStart[n3][en]) + if firstEndNodeScaleFactorIndex: + scaleFactors.append(scaleFactorMapEnd[n3][e1]) + scaleFactors.append(scaleFactorMapEnd[n3][en]) + if mapStartLinearDerivativeXi3: eftFactory.setEftLinearDerivative2(eft1, [ 1, 5, 2, 6 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapStartDerivatives: for i in range(2): lns = [ 1, 5 ] if (i == 0) else [ 2, 6 ] for n3 in range(2): - derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] + derivativesMap = (startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en]) if startDerivativesMap else (None, None, None) # handle different d1 on each side of node d1Map = derivativesMap[0] if ((i == 1) or (len(derivativesMap) < 4)) else derivativesMap[3] - d2Map = derivativesMap[1] + d2Map = derivativesMap[1] if derivativesMap[1] else (0, 1, 0) d3Map = derivativesMap[2] # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] - if d1Map is not None: + if d1Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) - if d3Map is not None: + if d3Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) - if d2Map is not None: + if d2Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, \ - derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d2Map)) - if d1Map is not None: + derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3 ), d2Map, + (firstStartNodeScaleFactorIndex + i + n3*layerNodeScaleFactorIndexOffset) if rescaleStartDerivatives else None) ) + if d1Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D2_DS1DS2, \ derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d1Map)) - if d3Map is not None: + if d3Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D2_DS2DS3, \ derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d3Map)) if mapEndLinearDerivativeXi3: @@ -418,26 +487,31 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] for n3 in range(2): - derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] + derivativesMap = (endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en]) if endDerivativesMap else (None, None, None) # handle different d1 on each side of node d1Map = derivativesMap[0] if ((i == 1) or (len(derivativesMap) < 4)) else derivativesMap[3] - d2Map = derivativesMap[1] + d2Map = derivativesMap[1] if derivativesMap[1] else (0, 1, 0) d3Map = derivativesMap[2] + # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] - if d1Map is not None: + if d1Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) - if d3Map is not None: + if d3Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) - if d2Map is not None: + if d2Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, \ - derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d2Map)) - if d1Map is not None: + derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3), d2Map, + (firstEndNodeScaleFactorIndex + i + n3*layerNodeScaleFactorIndexOffset) if rescaleEndDerivatives else None) ) + if d1Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D2_DS1DS2, \ derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d1Map)) - if d3Map is not None: + if d3Map: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D2_DS2DS3, \ derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D_DS3 ), d3Map)) + elementtemplateX.defineField(coordinates, -1, eft1) elementtemplate1 = elementtemplateX else: @@ -446,11 +520,9 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, element = mesh.createElement(elementIdentifier, elementtemplate1) result2 = element.setNodesByIdentifier(eft1, nids) - if mapDerivatives: - result3 = element.setScaleFactors(eft1, [ -1.0 ]) - #else: - # result3 = '-' - #print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + if scaleFactors: + result3 = element.setScaleFactors(eft1, scaleFactors); + #print('create element annulus', element.isValid(), elementIdentifier, eft1.validate(), result2, result3 if scaleFactors else None, nids) elementIdentifier += 1 if rowMeshGroups: @@ -460,3 +532,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, fm.endChange() return nodeIdentifier, elementIdentifier + +def getQuadrantID(d): + """ + Returns a scale factor ID based on direction of the derivative. Index starts + at 1 in direction (1, 0, 0) and progresses in anti-clockwise direction. + :param d: directional derivative + :return: scale factor ID + """ + maps = [(1, 0, 0), (1, 1, 0), (0, 1, 0), (-1, 1, 0), (-1, 0, 0), (-1, -1, 0), (0, -1, 0), (1, -1, 0), # d3 = 0 + (1, 0, 1), (1, 1, 1), (0, 1, 1), (-1, 1, 1), (-1, 0, 1), (-1, -1, 1), (0, -1, 1), (1, -1, 1), (0, 0, 1), # d3 = 1 + (1, 0, -1), (1, 1, -1), (0, 1, -1), (-1, 1, -1), (-1, 0, -1), (-1, -1, -1), (0, -1, -1), (1, -1, -1), (0, 0, -1)] # d3 = -1 + for i in range(len(maps)): + if d == maps[i]: + return i + 1 diff --git a/src/scaffoldmaker/utils/eft_utils.py b/src/scaffoldmaker/utils/eft_utils.py index 14ea709b..bcb69f11 100644 --- a/src/scaffoldmaker/utils/eft_utils.py +++ b/src/scaffoldmaker/utils/eft_utils.py @@ -117,13 +117,13 @@ def scaleEftNodeValueLabels(eft, localNodeIndexes, valueLabels, addScaleFactorIn scaleFactorIndexes += addScaleFactorIndexes eft.setTermScaling(f, t, scaleFactorIndexes) -def setEftScaleFactorIds(eft, generalScaleFactorIds, nodeScaleFactorIds): +def setEftScaleFactorIds(eft, globalScaleFactorIds, nodeScaleFactorIds): ''' Set general followed by node scale factor identifiers. ''' - eft.setNumberOfLocalScaleFactors(len(generalScaleFactorIds) + len(nodeScaleFactorIds)) + eft.setNumberOfLocalScaleFactors(len(globalScaleFactorIds) + len(nodeScaleFactorIds)) s = 1 - for id in generalScaleFactorIds: + for id in globalScaleFactorIds: eft.setScaleFactorType(s, Elementfieldtemplate.SCALE_FACTOR_TYPE_GLOBAL_GENERAL) eft.setScaleFactorIdentifier(s, id) s += 1 diff --git a/tests/test_bladder.py b/tests/test_bladder.py index 47d0ec57..ddea6c1b 100644 --- a/tests/test_bladder.py +++ b/tests/test_bladder.py @@ -87,7 +87,7 @@ def test_bladder1(self): self.assertTrue(region.isValid()) annotationGroups = MeshType_3d_bladder1.generateBaseMesh(region, options) - self.assertEqual(2, len(annotationGroups)) + self.assertEqual(3, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) if annotationGroups is not None: diff --git a/tests/test_cecum.py b/tests/test_cecum.py index 1eeab44c..37a933bd 100644 --- a/tests/test_cecum.py +++ b/tests/test_cecum.py @@ -79,10 +79,10 @@ def test_cecum1(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 65960.02821062482, delta=1.0E-6) + self.assertAlmostEqual(surfaceArea, 65959.34314511667, delta=1.0E-6) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 127893.74708048582, delta=1.0E-6) + self.assertAlmostEqual(volume, 127894.8597203706, delta=1.0E-6) if __name__ == "__main__": unittest.main() diff --git a/tests/test_colon.py b/tests/test_colon.py index b1afded5..6a6671d2 100644 --- a/tests/test_colon.py +++ b/tests/test_colon.py @@ -21,7 +21,7 @@ def test_colon1(self): Test creation of colon scaffold. """ parameterSetNames = MeshType_3d_colon1.getParameterSetNames() - self.assertEqual(parameterSetNames, ["Default", "Human 1", "Human 2", "Mouse 1", "Mouse 2", "Pig 1", "Pig 2"]) + self.assertEqual(parameterSetNames, ["Default", "Cattle 1", "Human 1", "Human 2", "Mouse 1", "Mouse 2", "Pig 1", "Pig 2"]) centralPathDefaultScaffoldPackages = { 'Test line': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { @@ -112,14 +112,14 @@ def test_colon1(self): coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement() self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) - assertAlmostEqualList(self, minimums, [ 108.0250647983898, -36.876103983560014, -25.89741158325083 ], 1.0E-6) - assertAlmostEqualList(self, maximums, [ 185.46457506220003, 48.101157490744, 34.995316052158934 ], 1.0E-6) + assertAlmostEqualList(self, minimums, [ 108.03959866945482, -36.876103983560014, -25.93217903595996 ], 1.0E-6) + assertAlmostEqualList(self, maximums, [ 185.43446761757468, 48.07433517752282, 34.995316052158934 ], 1.0E-6) flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement() self.assertTrue(flatCoordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(flatCoordinates, nodes) assertAlmostEqualList(self, minimums, [ 0.0, 0.0, 0.0 ], 1.0E-6) - assertAlmostEqualList(self, maximums, [ 186.72988844629867, 77.4178187926561, 3.2000000000000006 ], 1.0E-6) + assertAlmostEqualList(self, maximums, [ 186.72988844629867, 75.79769125771575, 3.2000000000000006 ], 1.0E-6) textureCoordinates = fieldmodule.findFieldByName("texture coordinates").castFiniteElement() minimums, maximums = evaluateFieldNodesetRange(textureCoordinates, nodes) @@ -136,10 +136,10 @@ def test_colon1(self): fieldcache = fieldmodule.createFieldcache() result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(surfaceArea, 14612.416788520026, delta=1.0E-6) + self.assertAlmostEqual(surfaceArea, 14623.340352853866, delta=1.0E-6) result, volume = volumeField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 26826.069540921028, delta=1.0E-6) + self.assertAlmostEqual(volume, 26869.74823158621, delta=1.0E-6) def test_mousecolon1(self): """ @@ -175,7 +175,7 @@ def test_mousecolon1(self): fieldcache = fieldmodule.createFieldcache() result, flatSurfaceArea = flatSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(flatSurfaceArea, 652.2383326727861, delta=1.0E-6) + self.assertAlmostEqual(flatSurfaceArea, 629.440514706522, delta=1.0E-6) result, textureSurfaceArea = textureSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) self.assertAlmostEqual(textureSurfaceArea, 1.0, delta=1.0E-6) diff --git a/tests/test_colonsegment.py b/tests/test_colonsegment.py index e7c82570..95cbb77a 100644 --- a/tests/test_colonsegment.py +++ b/tests/test_colonsegment.py @@ -16,7 +16,7 @@ def test_humancolonsegment1(self): Test creation of human colon segment scaffold. """ parameterSetNames = MeshType_3d_colonsegment1.getParameterSetNames() - self.assertEqual(parameterSetNames, [ "Default", "Human 1", "Mouse 1", "Pig 1" ]) + self.assertEqual(parameterSetNames, [ "Default", "Cattle 1", "Human 1", "Mouse 1", "Pig 1" ]) options = MeshType_3d_colonsegment1.getDefaultOptions("Human 1") self.assertEqual(27, len(options)) self.assertEqual(0.0, options.get("Start phase"))