From 64a60517b0b53507e56370d71f0fad023ad657b7 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Sat, 31 Oct 2020 18:07:52 +1300 Subject: [PATCH 01/18] Add rescale end derivatives to annulusMesh --- .../meshtypes/meshtype_3d_bladderurethra1.py | 28 ++++- src/scaffoldmaker/utils/annulusmesh.py | 102 ++++++++++++++++-- 2 files changed, 119 insertions(+), 11 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index a5d1b61a..2e004752 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -25,7 +25,6 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node - class MeshType_3d_bladderurethra1(Scaffold_base): ''' Generates 3D bladder and urethra meshes with variable numbers @@ -238,6 +237,7 @@ def getOrderedOptionNames(): 'Ureter', 'Ostium position around', 'Ostium position down', + 'Number of elements radially on annulus', 'Use cross derivatives', 'Use linear through wall', 'Refine', @@ -707,11 +707,13 @@ def generateBaseMesh(cls, region, options): nodesOnTrackSurface_d2.append(outerNodes_d2[n2 * elementsCountAround + n1]) trackSurfaceOstium1 = 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 nodesOnTrackSurface2_x = [] nodesOnTrackSurface2_d1 = [] @@ -724,6 +726,7 @@ 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, nodesOnTrackSurface2_d1, nodesOnTrackSurface2_d2) ostium2Position = TrackSurfacePosition(elementsCountAround - ostiumElementPositionAround, @@ -914,11 +917,32 @@ 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(trackSurfaceOstium1.getProportion(o1_Positions[n])) + + endProportions1 = [] + elementsAroundTrackSurface1 = trackSurfaceOstium1.elementsCount1 + elementsAlongTrackSurface1 = trackSurfaceOstium1.elementsCount2 + for n in range(3): + endProportions1.append([(ostiumElementPositionAround - 1 + 1)/elementsAroundTrackSurface1, + (ostiumElementPositionDown - 1 + n + 1)/elementsAlongTrackSurface1]) + for n in range(2): + endProportions1.append([(ostiumElementPositionAround + n + 1) / elementsAroundTrackSurface1, + (ostiumElementPositionDown + 1 + 1) / elementsAlongTrackSurface1]) + for n in range(2): + endProportions1.append([(ostiumElementPositionAround + 1 + 1) / elementsAroundTrackSurface1, + (ostiumElementPositionDown - n + 1) / elementsAlongTrackSurface1]) + endProportions1.append([(ostiumElementPositionAround + 1 )/ elementsAroundTrackSurface1, + (ostiumElementPositionDown - 1 + 1) / elementsAlongTrackSurface1]) + 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=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups, + tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, + rescaleEndDerivatives=True) nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 360f0ae6..895653b3 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -33,7 +33,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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, + rescaleEndDerivatives = False): """ 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 +73,8 @@ 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 rescaleEndDerivatives: If true, rescale end derivatives to arclength between endPoints + and radially adjacent nodes. :return: Final values of nextNodeIdentifier, nextElementIdentifier """ assert (elementsCountRadial >= 1), 'createAnnulusMesh3d: Invalid number of radial elements' @@ -134,6 +137,7 @@ 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 ] @@ -226,10 +230,11 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # scaling end derivatives to arc length gives even curvature along the curve if tracksurface: - mx, md2, md1 = \ + mx, md2, md1, md3, mProportions = \ tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], endProportions[n1][0], endProportions[n1][1], - elementsCountRadial, derivativeStart=ad2, derivativeEnd=bd2)[0:3] + elementsCountRadial, derivativeStart=ad2, derivativeEnd=bd2) + mx, md2, md1 = tracksurface.resampleHermiteCurvePointsSmooth(mx, md2, md1, md3, mProportions)[0:3] else: arcLength = interp.computeCubicHermiteArcLength(ax, ad2, bx, bd2, rescaleDerivatives = False) scaledDerivatives = [ vector.setMagnitude(d2, arcLength) for d2 in [ ad2, bd2 ]] @@ -307,6 +312,24 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n2 in range(elementsCountRadial + 1): pd2[n3][n2][n1] = sd2[n2] + # Calculate arcLength to determine scale factor for remapped d2 + if rescaleEndDerivatives: + scaleFactorMap = [] + for n3 in range(2): + scaleFactorList = [] + for n1 in range(nodesCountAround): + v1 = px[n3][-2][n1] + d1 = pd2[n3][-2][n1] + v2 = px[n3][-1][n1] + d2 = [v2[c] - v1[c] for c in range(3)] + arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) + unscaledEndDerivativeD2 = endDerivativesMap[n3][n1][1] + unscaledArcLength = vector.magnitude([abs(unscaledEndDerivativeD2[0]) * endPointsd1[n3][n1][c] + + abs(unscaledEndDerivativeD2[1]) * endPointsd2[n3][n1][c] for c in range(3)]) + scaleFactorNode = arcLength / unscaledArcLength + scaleFactorList.append(scaleFactorNode) + scaleFactorMap.append(scaleFactorList) + ############## # Create nodes ############## @@ -415,8 +438,25 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, if mapEndLinearDerivativeXi3: eftFactory.setEftLinearDerivative2(eft1, [ 3, 7, 4, 8 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapEndDerivatives: + if rescaleEndDerivatives: + nodeScaleFactorOffset0 = e1 * 1000 + nodeScaleFactorOffset1 = en * 1000 + setEftScaleFactorIds(eft1, [1], + [nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset1 + 1, + nodeScaleFactorOffset1 + 1]) + # nsCheck = [-1, + # nodeScaleFactorOffset0 + 1, + # nodeScaleFactorOffset0 + 1, + # nodeScaleFactorOffset1 + 1, + # nodeScaleFactorOffset1 + 1] + # sfCheck = [ -1.0, + # scaleFactorMap[0][e1], scaleFactorMap[1][e1], + # scaleFactorMap[0][en], scaleFactorMap[1][en]] for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] + sos = 1 if (i == 0) else 3 for n3 in range(2): derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] # handle different d1 on each side of node @@ -425,19 +465,57 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, d3Map = derivativesMap[2] # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] + so = sos + n3 + # if n3 == 0: + # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so) if d1Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) if d2Map is not None: - 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 rescaleEndDerivatives: + # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so, nsCheck[so], d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) + if d2Map == (-1, -1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1, so + 1]), + (Node.VALUE_LABEL_D_DS2, [1, so + 1])]) + elif d2Map == (1, 1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [so + 1]), + (Node.VALUE_LABEL_D_DS2, [so + 1])]) + elif d2Map == (-1, 1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1, so + 1]), + (Node.VALUE_LABEL_D_DS2, [so + 1])]) + elif d2Map == (1, -1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [so + 1]), + (Node.VALUE_LABEL_D_DS2, [1, so + 1])]) + elif d2Map == (-1, 0, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [1, so + 1])]) + elif d2Map == (1, 0, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS1, [so + 1])]) + elif d2Map == (0, -1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS2, [1, so + 1])]) + elif d2Map == (0, 1, 0): + remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, + [(Node.VALUE_LABEL_D_DS2, [so + 1])]) + else: + 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: 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: 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 +524,17 @@ 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 rescaleEndDerivatives: + result3 = element.setScaleFactors(eft1, [ -1.0, + scaleFactorMap[0][e1], scaleFactorMap[1][e1], + scaleFactorMap[0][en], scaleFactorMap[1][en]]) + # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + else: + result3 = element.setScaleFactors(eft1, [ -1.0 ]) + # print('else create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + elementIdentifier += 1 if rowMeshGroups: From 54ad709475140fa56579e77c50d327e5ca0f7fad Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 5 Nov 2020 12:53:35 +1300 Subject: [PATCH 02/18] Allow ureters to move within 2x2 elements --- .../meshtypes/meshtype_3d_bladderurethra1.py | 231 +++++++++++++----- src/scaffoldmaker/utils/annulusmesh.py | 27 +- 2 files changed, 189 insertions(+), 69 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 2e004752..cef9af2b 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -24,6 +24,7 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node +# from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates #km class MeshType_3d_bladderurethra1(Scaffold_base): ''' @@ -192,9 +193,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): '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, + 'Ostium position around': 0.67, # should be on the dorsal part (> 0.5) + 'Ostium position down': 0.83, + 'Number of elements radially on annulus': 2, 'Use cross derivatives': False, 'Use linear through wall': True, 'Refine': False, @@ -208,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['Ostium position around'] = 0.67 # should be on the dorsal part (> 0.5) + options['Ostium position down'] = 0.83 options['Urethra diameter 1'] = 0.75 options['Urethra diameter 2'] = 0.65 options['Urethra wall thickness'] = 0.25 @@ -300,6 +301,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['Ostium position around'] < 0.1: + options['Ostium position around'] = 0.1 + elif options['Ostium position around'] > 0.9: + options['Ostium position around'] = 0.9 + if options['Ostium position down'] < 0.15: + options['Ostium position down'] = 0.15 + elif options['Ostium position down'] > 0.95: + options['Ostium position down'] = 0.95 @classmethod def updateSubScaffoldOptions(cls, options): @@ -693,6 +702,22 @@ def generateBaseMesh(cls, region, options): outerNodes_d1.append(d1List[(2 * n2 + 1) * elementsCountAround + n1]) outerNodes_d2.append(d2List[(2 * n2 + 1) * elementsCountAround + n1]) + #################################################### + # zero = [0.0, 0.0, 0.0] + # fm = region.getFieldmodule() + # fm.beginChange() + # cache = fm.createFieldcache() + # + # # Coordinates field + # coordinates = findOrCreateFieldCoordinates(fm) + # nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + # nodetemplate = nodes.createNodetemplate() + # nodetemplate.defineField(coordinates) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + ########################################################### + if includeUreter: elementsCount1 = elementsCountAround // 2 elementsCount2 = elementsCountAlongBladder @@ -709,8 +734,6 @@ def generateBaseMesh(cls, region, options): nodesOnTrackSurface_d1, nodesOnTrackSurface_d2) ostium1Position = trackSurfaceOstium1.createPositionProportion(ostiumPositionAround, ostiumPositionDown) - ostium1Position.xi1 = 1.0 - ostium1Position.xi2 = 1.0 ostiumElementPositionAround = ostium1Position.e1 ostiumElementPositionDown = ostium1Position.e2 @@ -729,17 +752,32 @@ def generateBaseMesh(cls, region, options): trackSurfaceOstium2 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface2_x, nodesOnTrackSurface2_d1, nodesOnTrackSurface2_d2) - ostium2Position = TrackSurfacePosition(elementsCountAround - ostiumElementPositionAround, - ostiumElementPositionDown - 1, 0.0, 1.0) + ostium2Position = TrackSurfacePosition(elementsCountAround//2 - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0 else 0), + ostiumElementPositionDown, + (1 - ostium1Position.xi1) if ostium1Position.xi1 > 0 else ostium1Position.xi1, + ostium1Position.xi2) annulusMeshGroups = [neckMeshGroup, urinaryBladderMeshGroup] generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaultOptions, - elementsCountAround, elementsCountAroundOstium, trackSurfaceOstium1, + elementsCountAround, elementsCountThroughWall, elementsCountAroundOstium, + trackSurfaceOstium1, ostium1Position, trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, ostiumElementPositionAround, xFinal, d1Final, d2Final, nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadially, annulusMeshGroups) + # nodeIdentifier = 90000 + # for n in range(len(nodesOnTrackSurface2_x)): + # node = nodes.createNode(nodeIdentifier, nodetemplate) + # cache.setNode(node) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodesOnTrackSurface2_x[n]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, zero) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, zero) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, zero) + # print('NodeIdentifier = ', nodeIdentifier) # , xEnd[n]) #, d1[n], d2[n]) + # nodeIdentifier = nodeIdentifier + 1 + # fm.endChange() + fm.endChange() return annotationGroups @@ -811,31 +849,56 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): 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, + elementsCountThroughWall, elementsCountAroundOstium, trackSurfaceOstium1, + ostium1Position, trackSurfaceOstium2, ostium2Position, + ostiumElementPositionDown, ostiumElementPositionAround, + xBladder, d1Bladder, d2Bladder, nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadially, annulusMeshGroups = []): + # #################################################### + # zero = [0.0, 0.0, 0.0] + # fm = region.getFieldmodule() + # fm.beginChange() + # cache = fm.createFieldcache() + # + # # Coordinates field + # coordinates = findOrCreateFieldCoordinates(fm) + # nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + # nodetemplate = nodes.createNodetemplate() + # nodetemplate.defineField(coordinates) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + # ########################################################### + # Create ureters on the surface # Ureter 1 centerUreter1_x, centerUreter1_d1, centerUreter1_d2 = trackSurfaceOstium1.evaluateCoordinates(ostium1Position, 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 * (ostiumElementPositionDown - (1 if ostium1Position.xi2 > 0.5 else 2)) \ + + ostiumElementPositionAround + (1 if ostium1Position.xi1 > 0.5 else 0) + ureter1StartCornerx = xBladder[endPointStartId1 - 1] v1 = [(ureter1StartCornerx[c] - centerUreter1_x[c]) for c in range(3)] + # print(endPointStartId1) ostium1Direction = vector.crossproduct3(td3, v1) nodeIdentifier, elementIdentifier, (o1_x, o1_d1, o1_d2, _, o1_NodeId, o1_Positions) = \ generateOstiumMesh(region, ostiumDefaultOptions, trackSurfaceOstium1, ostium1Position, ostium1Direction, startNodeIdentifier=nextNodeIdentifier, startElementIdentifier=nextElementIdentifier, ostiumMeshGroups=None) + # Ureter 2 centerUreter2_x, centerUreter2_d1, centerUreter2_d2 = trackSurfaceOstium2.evaluateCoordinates(ostium2Position, 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 * \ + (ostiumElementPositionDown - (1 if ostium1Position.xi2 > 0.5 else 2)) \ + + elementsCountAround - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0.5 else 0) + # print(endPointStartId2) + ureter2StartCornerx = xBladder[endPointStartId2 - 1] v2 = [(ureter2StartCornerx[c] - centerUreter2_x[c]) for c in range(3)] ostium2Direction = vector.crossproduct3(ad3, v2) nodeIdentifier, elementIdentifier, (o2_x, o2_d1, o2_d2, _, o2_NodeId, o2_Positions) = \ @@ -853,39 +916,35 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul 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 + 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): 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] @@ -921,20 +980,41 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul for n in range(len(o1_Positions)): startProportions1.append(trackSurfaceOstium1.getProportion(o1_Positions[n])) + startProportions2 = [] + for n in range(len(o2_Positions)): + startProportions2.append(trackSurfaceOstium2.getProportion(o2_Positions[n])) + endProportions1 = [] elementsAroundTrackSurface1 = trackSurfaceOstium1.elementsCount1 elementsAlongTrackSurface1 = trackSurfaceOstium1.elementsCount2 + + endProportions2 = [] + elementsAroundTrackSurface2 = trackSurfaceOstium2.elementsCount1 + elementsAlongTrackSurface2 = trackSurfaceOstium2.elementsCount2 + + firstIdxAround1 = ostiumElementPositionAround + (0 if ostium1Position.xi1 > 0.5 else -1) + firstIdxAlong = ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1) + firstIdxAround2 = elementsCountAround//2 - ostiumElementPositionAround - (2 if ostium1Position.xi1 > 0.5 else 1) + for n in range(3): - endProportions1.append([(ostiumElementPositionAround - 1 + 1)/elementsAroundTrackSurface1, - (ostiumElementPositionDown - 1 + n + 1)/elementsAlongTrackSurface1]) + endProportions1.append([(firstIdxAround1)/elementsAroundTrackSurface1, + (firstIdxAlong + n)/elementsAlongTrackSurface1]) + endProportions2.append([(firstIdxAround2) / elementsAroundTrackSurface2, + (firstIdxAlong + n) / elementsAlongTrackSurface2]) for n in range(2): - endProportions1.append([(ostiumElementPositionAround + n + 1) / elementsAroundTrackSurface1, - (ostiumElementPositionDown + 1 + 1) / elementsAlongTrackSurface1]) + 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([(ostiumElementPositionAround + 1 + 1) / elementsAroundTrackSurface1, - (ostiumElementPositionDown - n + 1) / elementsAlongTrackSurface1]) - endProportions1.append([(ostiumElementPositionAround + 1 )/ elementsAroundTrackSurface1, - (ostiumElementPositionDown - 1 + 1) / elementsAlongTrackSurface1]) + 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, @@ -944,24 +1024,43 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, rescaleEndDerivatives=True) + # print('endNodeId = ', endNode2_Id) 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=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups, + tracksurface=trackSurfaceOstium2, 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 * (ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1)) \ + + ostiumElementPositionAround + (1 if ostium1Position.xi1 > 0.5 else 0) + elementToDeleteStartIdx2 = elementsCountThroughWall * elementsCountAround * (ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1)) \ + + elementsCountAround - ostiumElementPositionAround + (-1 if ostium1Position.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) + + # ####################################################### + # nodeIdentifier = 10000 + # for n in range(len(endPoints2_x[1])): + # node = nodes.createNode(nodeIdentifier, nodetemplate) + # cache.setNode(node) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, endPoints2_x[1][n]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, zero) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, zero) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, zero) + # print('NodeIdentifier = ', nodeIdentifier) #, xEnd[n]) #, d1[n], d2[n]) + # nodeIdentifier = nodeIdentifier + 1 + # fm.endChange() + # ############################################################ + return diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 895653b3..2dd76b2d 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -230,11 +230,20 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # scaling end derivatives to arc length gives even curvature along the curve if tracksurface: + arcLength = interp.computeCubicHermiteArcLength(startPointsx[1][n1], ad2, endPointsx[1][n1], bd2, + rescaleDerivatives=True) + ad2Scaled = vector.setMagnitude(ad2, arcLength / elementsCountRadial * 2 * 0.4) + bd2Scaled = vector.setMagnitude(bd2, arcLength / elementsCountRadial * 2 * 0.6) + mx, md2, md1, md3, mProportions = \ tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], endProportions[n1][0], endProportions[n1][1], - elementsCountRadial, derivativeStart=ad2, derivativeEnd=bd2) - mx, md2, md1 = tracksurface.resampleHermiteCurvePointsSmooth(mx, md2, md1, md3, mProportions)[0:3] + elementsCountRadial, derivativeStart=ad2Scaled, + derivativeEnd=bd2Scaled) + + mx, md2, md1 = tracksurface.resampleHermiteCurvePointsSmooth(mx, md2, md1, md3, mProportions, + derivativeMagnitudeStart= vector.magnitude(ad2Scaled), + derivativeMagnitudeEnd= vector.magnitude(bd2Scaled))[0:3] else: arcLength = interp.computeCubicHermiteArcLength(ax, ad2, bx, bd2, rescaleDerivatives = False) scaledDerivatives = [ vector.setMagnitude(d2, arcLength) for d2 in [ ad2, bd2 ]] @@ -321,7 +330,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, v1 = px[n3][-2][n1] d1 = pd2[n3][-2][n1] v2 = px[n3][-1][n1] - d2 = [v2[c] - v1[c] for c in range(3)] + d2 = [pd1[n3][-1][n1][c] + pd2[n3][-1][n1][c] for c in range(3)] arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) unscaledEndDerivativeD2 = endDerivativesMap[n3][n1][1] unscaledArcLength = vector.magnitude([abs(unscaledEndDerivativeD2[0]) * endPointsd1[n3][n1][c] + @@ -531,6 +540,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, scaleFactorMap[0][e1], scaleFactorMap[1][e1], scaleFactorMap[0][en], scaleFactorMap[1][en]]) # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + else: result3 = element.setScaleFactors(eft1, [ -1.0 ]) # print('else create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) @@ -541,6 +551,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for meshGroup in rowMeshGroups[e2]: meshGroup.addElement(element) + # if rescaleEndDerivatives: + # nodeIdentifier = 90000 + # for n1 in range(len(endPointsx[1])): + # node = nodes.createNode(nodeIdentifier, nodetemplate) + # cache.setNode(node) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, endPointsx[1][n1]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, checkDerivativeStartList[n1]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, checkDerivativeEndList[n1]) + # print(nodeIdentifier) + # nodeIdentifier += 1 + fm.endChange() return nodeIdentifier, elementIdentifier From f1289eb7fec051ff4c4b86c3d4543431eb58dd12 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Fri, 6 Nov 2020 11:27:45 +1300 Subject: [PATCH 03/18] Modify derivativeSignsToExpressionTerms to include scaling by scale factor --- src/scaffoldmaker/utils/annulusmesh.py | 51 ++++++-------------------- 1 file changed, 11 insertions(+), 40 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 2dd76b2d..008a562b 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -14,18 +14,19 @@ 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: List of 1, -1 (scale by scale factor indexed by scaleFactorIdx) or 0 (no term). ''' 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 createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, @@ -482,42 +483,12 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) if d2Map is not None: - if rescaleEndDerivatives: - # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so, nsCheck[so], d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) - if d2Map == (-1, -1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [1, so + 1]), - (Node.VALUE_LABEL_D_DS2, [1, so + 1])]) - elif d2Map == (1, 1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [so + 1]), - (Node.VALUE_LABEL_D_DS2, [so + 1])]) - elif d2Map == (-1, 1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [1, so + 1]), - (Node.VALUE_LABEL_D_DS2, [so + 1])]) - elif d2Map == (1, -1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [so + 1]), - (Node.VALUE_LABEL_D_DS2, [1, so + 1])]) - elif d2Map == (-1, 0, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [1, so + 1])]) - elif d2Map == (1, 0, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS1, [so + 1])]) - elif d2Map == (0, -1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS2, [1, so + 1])]) - elif d2Map == (0, 1, 0): - remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, - [(Node.VALUE_LABEL_D_DS2, [so + 1])]) - else: - 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)) + # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so, nsCheck[so], d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) + 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, so + 1 if rescaleEndDerivatives else None)) if d1Map is not None: 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)) From bce0e46b1377f93ab7c346bf850891b200f08e84 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 11 Nov 2020 09:16:28 +1300 Subject: [PATCH 04/18] Add option to rescaleStartDerivatives --- .../meshtypes/meshtype_3d_bladderurethra1.py | 39 ++++ src/scaffoldmaker/utils/annulusmesh.py | 203 ++++++++++++++---- 2 files changed, 203 insertions(+), 39 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index cef9af2b..38664e6d 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -911,6 +911,8 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul endPoints1_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endNode1_Id = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endDerivativesMap = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] + # startDerivativesMap = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] # KM + # endDerivativesMapTest = [[((1, 0, 0), (0, -1, 0), None)] * elementsCountAroundOstium, [((1, 0, 0), (0, -1, 0), None)] * elementsCountAroundOstium] # KM endPoints2_x = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endPoints2_d1 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endPoints2_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] @@ -950,6 +952,32 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul endPoints2_d1[n3][n1] = d1Bladder[nc2] endPoints2_d2[n3][n1] = d2Bladder[nc2] + # for n1 in range(elementsCountAroundOstium): + # if n1 == 0: + # startDerivativesMap[0][n1] = ((-1, 0, 0), (1, 1, 0), None, (0, 1, 0)) + # startDerivativesMap[1][n1] = ((-1, 0, 0), (1, 1, 0), None, (0, 1, 0)) + # elif n1 == 1: + # startDerivativesMap[0][n1] = ((0, 1, 0), (1, 0, 0), None) + # startDerivativesMap[1][n1] = ((0, 1, 0), (1, 0, 0), None) + # elif n1 == 2: + # startDerivativesMap[0][n1] = ((0, 1, 0), (1, -1, 0), None, (1, 0, 0)) + # startDerivativesMap[1][n1] = ((0, 1, 0), (1, -1, 0), None, (1, 0, 0)) + # elif n1 == 3: + # startDerivativesMap[0][n1] = ((1, 0, 0), (0, -1, 0), None) + # startDerivativesMap[1][n1] = ((1, 0, 0), (0, -1, 0), None) + # elif n1 == 4: + # startDerivativesMap[0][n1] = ((1, 0, 0), (-1, -1, 0), None, (0, -1, 0)) + # startDerivativesMap[1][n1] = ((1, 0, 0), (-1, -1, 0), None, (0, -1, 0)) + # elif n1 == 5: + # startDerivativesMap[0][n1] = ((0, -1, 0), (-1, 0, 0), None) + # startDerivativesMap[1][n1] = ((0, -1, 0), (-1, 0, 0), None) + # elif n1 == 6: + # startDerivativesMap[0][n1] = ((0, -1, 0), (-1, 1, 0), None, (-1, 0, 0)) + # startDerivativesMap[1][n1] = ((0, -1, 0), (-1, 1, 0), None, (-1, 0, 0)) + # else: + # startDerivativesMap[0][n1] = ((-1, 0, 0), (0, 1, 0), None) + # startDerivativesMap[1][n1] = ((-1, 0, 0), (0, 1, 0), None) + for n1 in range(elementsCountAroundOstium): if n1 == 0: endDerivativesMap[0][n1] = ((-1, 0, 0), (-1, -1, 0), None, (0, 1, 0)) @@ -976,6 +1004,7 @@ 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(trackSurfaceOstium1.getProportion(o1_Positions[n])) @@ -1024,6 +1053,16 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, rescaleEndDerivatives=True) + # print('Making test annulus') + # nodeIdentifier, elementIdentifier = createAnnulusMesh3d( + # nodes, mesh, nodeIdentifier, elementIdentifier, + # endPoints1_x, endPoints1_d1, endPoints1_d2, None, endNode1_Id, startDerivativesMap, + # o1_x, o1_d1, o1_d2, None, o1_NodeId, endDerivativesMapTest, + # elementsCountRadial=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups, + # tracksurface=trackSurfaceOstium1, startProportions=endProportions1, endProportions=startProportions1, + # rescaleStartDerivatives=True, rescaleEndDerivatives = True, + # debugTest=True) + # print('endNodeId = ', endNode2_Id) nodeIdentifier, elementIdentifier = createAnnulusMesh3d( nodes, mesh, nodeIdentifier, elementIdentifier, diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 008a562b..ec2d5886 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -12,6 +12,8 @@ from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import vector +from scaffoldmaker.utils import matrix #KM +import math # KM def derivativeSignsToExpressionTerms(valueLabels, signs, scaleFactorIdx = None): @@ -35,7 +37,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, forceStartLinearXi3 = False, forceMidLinearXi3 = False, forceEndLinearXi3 = False, maxStartThickness = None, maxEndThickness = None, useCrossDerivatives = False, elementsCountRadial = 1, meshGroups = None, tracksurface = None, startProportions = None, endProportions = None, - rescaleEndDerivatives = False): + rescaleStartDerivatives = False, rescaleEndDerivatives = False, debugTest = False): """ 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. @@ -74,6 +76,8 @@ 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: If true, rescale start derivatives to arclength between startPoints + and radially adjacent nodes. :param rescaleEndDerivatives: If true, rescale end derivatives to arclength between endPoints and radially adjacent nodes. :return: Final values of nextNodeIdentifier, nextElementIdentifier @@ -127,6 +131,24 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, cache = fm.createFieldcache() coordinates = findOrCreateFieldCoordinates(fm) + ############### MOVE BACK DOWN ################################################## + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) + + if debugTest: + nodeIdentifier = 70000 + ##################################################################################### + # Build arrays of points from start to end px = [ [], [] ] pd1 = [ [], [] ] @@ -233,8 +255,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, if tracksurface: arcLength = interp.computeCubicHermiteArcLength(startPointsx[1][n1], ad2, endPointsx[1][n1], bd2, rescaleDerivatives=True) - ad2Scaled = vector.setMagnitude(ad2, arcLength / elementsCountRadial * 2 * 0.4) - bd2Scaled = vector.setMagnitude(bd2, arcLength / elementsCountRadial * 2 * 0.6) + ad2Scaled = vector.setMagnitude(ad2, arcLength / elementsCountRadial * 2 * (0.4 if vector.magnitude(ad2) < vector.magnitude(bd2) else 0.6)) + bd2Scaled = vector.setMagnitude(bd2, arcLength / elementsCountRadial * 2 * (0.6 if vector.magnitude(ad2) < vector.magnitude(bd2) else 0.4)) mx, md2, md1, md3, mProportions = \ tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], @@ -254,6 +276,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, md1 = interp.interpolateSampleLinear([ ad1, bd1 ], me, mxi) thi = interp.interpolateSampleLinear([ thicknesses[0][n1], thicknesses[-1][n1] ], me, mxi) + # if n1 == 0: + # if debugTest or debug2: + # for i in range(len(mx)): + # node = nodes.createNode(nodeIdentifier, nodetemplate) + # cache.setNode(node) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, mx[i]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, md1[i]) + # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, md2[i]) + # print(nodeIdentifier, vector.normalise(md1[i])) + # nodeIdentifier += 1 + #md2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, fixStartDerivative = True, fixEndDerivative = True) for n2 in range(1, elementsCountRadial): @@ -283,6 +316,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n1 in range(nodesCountAround): normal = vector.normalise(vector.crossproduct3(pd1[1][n2][n1], pd2[1][n2][n1])) + if debugTest: + normal = vector.normalise(vector.crossproduct3(pd2[1][n2][n1], pd1[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) ] @@ -323,8 +358,23 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, pd2[n3][n2][n1] = sd2[n2] # Calculate arcLength to determine scale factor for remapped d2 + if rescaleStartDerivatives: + scaleFactorMapStart = [] + for n3 in range(2): + scaleFactorList = [] + for n1 in range(nodesCountAround): + v1 = px[n3][0][n1] + d1 = pd2[n3][0][n1] + v2 = px[n3][1][n1] + d2 = pd2[n3][1][n1] if elementsCountRadial > 1 else [pd1[n3][-1][n1][c] + pd2[n3][-1][n1][c] for c in range(3)] + arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) + unscaledArcLength = vector.magnitude(pd2[n3][0][1]) + scaleFactorNode = arcLength / unscaledArcLength + scaleFactorList.append(scaleFactorNode) + scaleFactorMapStart.append(scaleFactorList) + if rescaleEndDerivatives: - scaleFactorMap = [] + scaleFactorMapEnd = [] for n3 in range(2): scaleFactorList = [] for n1 in range(nodesCountAround): @@ -338,24 +388,24 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, abs(unscaledEndDerivativeD2[1]) * endPointsd2[n3][n1][c] for c in range(3)]) scaleFactorNode = arcLength / unscaledArcLength scaleFactorList.append(scaleFactorNode) - scaleFactorMap.append(scaleFactorList) + scaleFactorMapEnd.append(scaleFactorList) ############## # Create nodes ############## - nodetemplate = nodes.createNodetemplate() - nodetemplate.defineField(coordinates) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) + # nodetemplate = nodes.createNodetemplate() + # nodetemplate.defineField(coordinates) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + # if useCrossDerivatives: + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + # if useCrossDerivatives: + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) nodetemplateLinearS3 = nodes.createNodetemplate() nodetemplateLinearS3.defineField(coordinates) nodetemplateLinearS3.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) @@ -384,6 +434,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, pd2[n3][n2][n1]) if pd3[n3][n2]: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, pd3[n3][n2][n1]) + # if nodeIdentifier == 402 or nodeIdentifier == 418: + # print('Node', nodeIdentifier, px[n3][n2][n1]) nodeIdentifier = nodeIdentifier + 1 nodeId[n3].append(rowNodeId) @@ -415,15 +467,49 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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] ] - + # if debugTest: + # print('Element', elementIdentifier, nids) if mapDerivatives: eft1 = eftFactory.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) if mapStartLinearDerivativeXi3: eftFactory.setEftLinearDerivative2(eft1, [ 1, 5, 2, 6 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapStartDerivatives: + if rescaleStartDerivatives: + if elementsCountRadial == 1 and rescaleEndDerivatives: + nodeScaleFactorOffset0 = e1 * 1000 + nodeScaleFactorOffset1 = en * 1000 + setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1, + nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) + # sfIDCheck = [1, + # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1, + # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1] + # sfCheck = [-1.0, + # scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], + # scaleFactorMapStart[0][en], scaleFactorMapStart[1][en], + # scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], + # scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]] + # if debugTest: + # print('scaleFactorIds') + # for i in range(len(sfIDCheck)): + # print(sfIDCheck[i]) + # print('scaleFactor') + # for i in range(len(sfCheck)): + # print(sfCheck[i]) + + else: + nodeScaleFactorOffset0 = e1 * 1000 + nodeScaleFactorOffset1 = en * 1000 + setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) + for i in range(2): lns = [ 1, 5 ] if (i == 0) else [ 2, 6 ] + sos = 1 if (i == 0) else 3 for n3 in range(2): derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] # handle different d1 on each side of node @@ -432,41 +518,60 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, d3Map = derivativesMap[2] # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] + so = sos + n3 if d1Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) if d2Map is not None: + # if debugTest: + # # print('mapStartDerivatives - ', ln, 'remap ds2') + # print('rescaleStart', elementIdentifier, ln, sfIDCheck[so], sfCheck[so]) 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)) + derivativeSignsToExpressionTerms( ( Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3 ), d2Map, + so + 1 if rescaleStartDerivatives else None)) if d1Map is not None: + # if debugTest: + # print('mapStartDerivatives - ', ln, 'remap ds1') 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 debugTest: + # print('mapStartDerivatives - ', ln, 'remap ds3') 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: eftFactory.setEftLinearDerivative2(eft1, [ 3, 7, 4, 8 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapEndDerivatives: if rescaleEndDerivatives: - nodeScaleFactorOffset0 = e1 * 1000 - nodeScaleFactorOffset1 = en * 1000 - setEftScaleFactorIds(eft1, [1], - [nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset1 + 1, - nodeScaleFactorOffset1 + 1]) - # nsCheck = [-1, - # nodeScaleFactorOffset0 + 1, - # nodeScaleFactorOffset0 + 1, - # nodeScaleFactorOffset1 + 1, - # nodeScaleFactorOffset1 + 1] - # sfCheck = [ -1.0, - # scaleFactorMap[0][e1], scaleFactorMap[1][e1], - # scaleFactorMap[0][en], scaleFactorMap[1][en]] + if elementsCountRadial == 1 and rescaleStartDerivatives: + pass + else: + nodeScaleFactorOffset0 = e1 * 1000 + nodeScaleFactorOffset1 = en * 1000 + setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) + # sfIDCheck = [1, + # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, + # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1] + # sfCheck = [-1.0, + # scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], + # scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]] + # + # print('rescaleEnd - scaleFactorIds') + # for i in range(len(sfIDCheck)): + # print(sfIDCheck[i]) + # print('scaleFactor') + # for i in range(len(sfCheck)): + # print(sfCheck[i]) for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] - sos = 1 if (i == 0) else 3 + if elementsCountRadial == 1 and rescaleStartDerivatives: + sos = 5 if (i == 0) else 7 + else: + sos = 1 if (i == 0) else 3 for n3 in range(2): derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] # handle different d1 on each side of node @@ -483,16 +588,23 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) if d2Map is not None: - # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so, nsCheck[so], d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) + # if debugTest: + # # print('mapEndDerivatives - ', ln, 'remap ds2') + # #print(elementIdentifier, e1, nodeScaleFactorOffset0, en, nodeScaleFactorOffset1) + # print('rescaleEnd', elementIdentifier, ln, sfIDCheck[so], sfCheck[so]) #, d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) 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, so + 1 if rescaleEndDerivatives else None)) if d1Map is not None: + # if debugTest: + # print('mapEndDerivatives - ', ln, 'remap ds1') 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 debugTest: + # print('mapEndDerivatives - ', ln, 'remap ds3') 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)) @@ -506,10 +618,23 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, result2 = element.setNodesByIdentifier(eft1, nids) if mapDerivatives: - if rescaleEndDerivatives: + if elementsCountRadial == 1 and rescaleStartDerivatives and rescaleEndDerivatives: + result3 = element.setScaleFactors(eft1, [-1.0, + scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], + scaleFactorMapStart[0][en], scaleFactorMapStart[1][en], + scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], + scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) + # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + elif rescaleStartDerivatives: + result3 = element.setScaleFactors(eft1, [-1.0, + scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], + scaleFactorMapStart[0][en], scaleFactorMapStart[1][en]]) + # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + + elif rescaleEndDerivatives: result3 = element.setScaleFactors(eft1, [ -1.0, - scaleFactorMap[0][e1], scaleFactorMap[1][e1], - scaleFactorMap[0][en], scaleFactorMap[1][en]]) + scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], + scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) else: From 9a5b1b68e6aaa1edd12d53050b2ab19e346a9b22 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Wed, 11 Nov 2020 16:53:28 +1300 Subject: [PATCH 05/18] Deal with cases where derivativeMap for d2 is None --- src/scaffoldmaker/utils/annulusmesh.py | 37 +++++++++++++++++--------- 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index ec2d5886..b7bec09d 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -145,8 +145,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) - if debugTest: - nodeIdentifier = 70000 + # if debugTest: + # nodeIdentifier = 70000 ##################################################################################### # Build arrays of points from start to end @@ -363,12 +363,18 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n3 in range(2): scaleFactorList = [] for n1 in range(nodesCountAround): + startDerivativeMapD2 = (0, 1, 0) if startDerivativesMap[n3][n1][1] is None else startDerivativesMap[n3][n1][1] v1 = px[n3][0][n1] - d1 = pd2[n3][0][n1] + d1 = [startDerivativeMapD2[0] * pd1[n3][0][n1][c] + startDerivativeMapD2[1] * pd2[n3][0][n1][c] for c in range(3)] v2 = px[n3][1][n1] - d2 = pd2[n3][1][n1] if elementsCountRadial > 1 else [pd1[n3][-1][n1][c] + pd2[n3][-1][n1][c] for c in range(3)] + d2 = pd2[n3][1][n1] + if elementsCountRadial == 1 and endDerivativesMap is not None: + endDerivativeMapD2 = (0, 1, 0) if endDerivativesMap[n3][n1][1] is None else endDerivativesMap[n3][n1][1] + d2 = [endDerivativeMapD2[0] * pd1[n3][1][n1][c] + endDerivativeMapD2[1] * pd2[n3][1][n1][c] for c in range(3)] arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) - unscaledArcLength = vector.magnitude(pd2[n3][0][1]) + unscaledArcLength = vector.magnitude([startDerivativeMapD2[0] * startPointsd1[n3][n1][c] + + startDerivativeMapD2[1] * startPointsd2[n3][n1][c] for c in + range(3)]) scaleFactorNode = arcLength / unscaledArcLength scaleFactorList.append(scaleFactorNode) scaleFactorMapStart.append(scaleFactorList) @@ -378,14 +384,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n3 in range(2): scaleFactorList = [] for n1 in range(nodesCountAround): + endDerivativesMapD2 = (0, 1, 0) if endDerivativesMap[n3][n1][1] is None else endDerivativesMap[n3][n1][1] v1 = px[n3][-2][n1] d1 = pd2[n3][-2][n1] + if elementsCountRadial == 1 and startDerivativesMap is not None: + startDerivativeMapD2 = (0, 1, 0) if startDerivativesMap[n3][n1][1] is None else startDerivativesMap[n3][n1][1] + d1 = [startDerivativeMapD2[0] * pd1[n3][0][n1][c] + startDerivativeMapD2[1] * pd2[n3][0][n1][c] for c in range(3)] v2 = px[n3][-1][n1] - d2 = [pd1[n3][-1][n1][c] + pd2[n3][-1][n1][c] for c in range(3)] + d2 = [endDerivativesMapD2[0] * pd1[n3][-1][n1][c] + endDerivativesMapD2[1] * pd2[n3][-1][n1][c] for c in range(3)] arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) - unscaledEndDerivativeD2 = endDerivativesMap[n3][n1][1] - unscaledArcLength = vector.magnitude([abs(unscaledEndDerivativeD2[0]) * endPointsd1[n3][n1][c] + - abs(unscaledEndDerivativeD2[1]) * endPointsd2[n3][n1][c] for c in range(3)]) + unscaledArcLength = vector.magnitude([endDerivativesMapD2[0] * endPointsd1[n3][n1][c] + + endDerivativesMapD2[1] * endPointsd2[n3][n1][c] for c in range(3)]) scaleFactorNode = arcLength / unscaledArcLength scaleFactorList.append(scaleFactorNode) scaleFactorMapEnd.append(scaleFactorList) @@ -468,7 +477,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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] ] # if debugTest: - # print('Element', elementIdentifier, nids) + # print('Element', elementIdentifier, nids) if mapDerivatives: eft1 = eftFactory.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) @@ -523,6 +532,8 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) + if d2Map is None: + d2Map = (0, 1, 0) if d2Map is not None: # if debugTest: # # print('mapStartDerivatives - ', ln, 'remap ds2') @@ -559,7 +570,6 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # sfCheck = [-1.0, # scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], # scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]] - # # print('rescaleEnd - scaleFactorIds') # for i in range(len(sfIDCheck)): # print(sfIDCheck[i]) @@ -587,16 +597,19 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) + if d2Map is None: + d2Map = (0, 1, 0) if d2Map is not None: # if debugTest: # # print('mapEndDerivatives - ', ln, 'remap ds2') - # #print(elementIdentifier, e1, nodeScaleFactorOffset0, en, nodeScaleFactorOffset1) + # # print(elementIdentifier, nids[(ln[0] - 1)], sfIDCheck[so], sfCheck[so]) # print('rescaleEnd', elementIdentifier, ln, sfIDCheck[so], sfCheck[so]) #, d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) 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, so + 1 if rescaleEndDerivatives else None)) + if d1Map is not None: # if debugTest: # print('mapEndDerivatives - ', ln, 'remap ds1') From 62a37f0ef4ce16130ca2d160be3d71995648487d Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 12 Nov 2020 15:38:35 +1300 Subject: [PATCH 06/18] Determine scaleFactorID using direction of derivatives --- .../meshtypes/meshtype_3d_bladderurethra1.py | 119 +----------- src/scaffoldmaker/utils/annulusmesh.py | 181 +++++------------- 2 files changed, 57 insertions(+), 243 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 38664e6d..55826edb 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -24,7 +24,6 @@ from opencmiss.zinc.element import Element from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node -# from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates #km class MeshType_3d_bladderurethra1(Scaffold_base): ''' @@ -195,7 +194,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Ureter': copy.deepcopy(ostiumOption), 'Ostium position around': 0.67, # should be on the dorsal part (> 0.5) 'Ostium position down': 0.83, - 'Number of elements radially on annulus': 2, + 'Number of radial elements on annulus': 2, 'Use cross derivatives': False, 'Use linear through wall': True, 'Refine': False, @@ -238,7 +237,7 @@ def getOrderedOptionNames(): 'Ureter', 'Ostium position around', 'Ostium position down', - 'Number of elements radially on annulus', + 'Number of radial elements on annulus', 'Use cross derivatives', 'Use linear through wall', 'Refine', @@ -352,7 +351,7 @@ def generateBaseMesh(cls, region, options): ostiumOptions = options['Ureter'] ostiumDefaultOptions = ostiumOptions.getScaffoldSettings() elementsCountAroundOstium = ostiumDefaultOptions['Number of elements around ostium'] - elementsCountAnnulusRadially = options['Number of elements radially on annulus'] + elementsCountAnnulusRadial = options['Number of radial elements on annulus'] ostiumPositionAround = options['Ostium position around'] ostiumPositionDown = options['Ostium position down'] @@ -702,22 +701,6 @@ def generateBaseMesh(cls, region, options): outerNodes_d1.append(d1List[(2 * n2 + 1) * elementsCountAround + n1]) outerNodes_d2.append(d2List[(2 * n2 + 1) * elementsCountAround + n1]) - #################################################### - # zero = [0.0, 0.0, 0.0] - # fm = region.getFieldmodule() - # fm.beginChange() - # cache = fm.createFieldcache() - # - # # Coordinates field - # coordinates = findOrCreateFieldCoordinates(fm) - # nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - # nodetemplate = nodes.createNodetemplate() - # nodetemplate.defineField(coordinates) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - ########################################################### - if includeUreter: elementsCount1 = elementsCountAround // 2 elementsCount2 = elementsCountAlongBladder @@ -757,27 +740,15 @@ def generateBaseMesh(cls, region, options): (1 - ostium1Position.xi1) if ostium1Position.xi1 > 0 else ostium1Position.xi1, ostium1Position.xi2) annulusMeshGroups = [neckMeshGroup, urinaryBladderMeshGroup] - generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaultOptions, + generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOptions, elementsCountAround, elementsCountThroughWall, elementsCountAroundOstium, trackSurfaceOstium1, ostium1Position, trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, ostiumElementPositionAround, xFinal, d1Final, d2Final, nextNodeIdentifier, - nextElementIdentifier, elementsCountAnnulusRadially, + nextElementIdentifier, elementsCountAnnulusRadial, annulusMeshGroups) - # nodeIdentifier = 90000 - # for n in range(len(nodesOnTrackSurface2_x)): - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, nodesOnTrackSurface2_x[n]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, zero) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, zero) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, zero) - # print('NodeIdentifier = ', nodeIdentifier) # , xEnd[n]) #, d1[n], d2[n]) - # nodeIdentifier = nodeIdentifier + 1 - # fm.endChange() - fm.endChange() return annotationGroups @@ -848,30 +819,14 @@ 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, +def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOptions,elementsCountAround, elementsCountThroughWall, elementsCountAroundOstium, trackSurfaceOstium1, ostium1Position, trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, ostiumElementPositionAround, xBladder, d1Bladder, d2Bladder, - nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadially, + nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadial, annulusMeshGroups = []): - # #################################################### - # zero = [0.0, 0.0, 0.0] - # fm = region.getFieldmodule() - # fm.beginChange() - # cache = fm.createFieldcache() - # - # # Coordinates field - # coordinates = findOrCreateFieldCoordinates(fm) - # nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - # nodetemplate = nodes.createNodetemplate() - # nodetemplate.defineField(coordinates) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - # ########################################################### - # Create ureters on the surface # Ureter 1 centerUreter1_x, centerUreter1_d1, centerUreter1_d2 = trackSurfaceOstium1.evaluateCoordinates(ostium1Position, @@ -882,7 +837,6 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul + ostiumElementPositionAround + (1 if ostium1Position.xi1 > 0.5 else 0) ureter1StartCornerx = xBladder[endPointStartId1 - 1] v1 = [(ureter1StartCornerx[c] - centerUreter1_x[c]) for c in range(3)] - # print(endPointStartId1) ostium1Direction = vector.crossproduct3(td3, v1) nodeIdentifier, elementIdentifier, (o1_x, o1_d1, o1_d2, _, o1_NodeId, o1_Positions) = \ generateOstiumMesh(region, ostiumDefaultOptions, trackSurfaceOstium1, ostium1Position, ostium1Direction, @@ -897,7 +851,6 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul + (elementsCountThroughWall + 1) * elementsCountAround * \ (ostiumElementPositionDown - (1 if ostium1Position.xi2 > 0.5 else 2)) \ + elementsCountAround - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0.5 else 0) - # print(endPointStartId2) ureter2StartCornerx = xBladder[endPointStartId2 - 1] v2 = [(ureter2StartCornerx[c] - centerUreter2_x[c]) for c in range(3)] ostium2Direction = vector.crossproduct3(ad3, v2) @@ -911,8 +864,6 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul endPoints1_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endNode1_Id = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endDerivativesMap = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] - # startDerivativesMap = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] # KM - # endDerivativesMapTest = [[((1, 0, 0), (0, -1, 0), None)] * elementsCountAroundOstium, [((1, 0, 0), (0, -1, 0), None)] * elementsCountAroundOstium] # KM endPoints2_x = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endPoints2_d1 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] endPoints2_d2 = [[None] * elementsCountAroundOstium, [None] * elementsCountAroundOstium] @@ -952,32 +903,6 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul endPoints2_d1[n3][n1] = d1Bladder[nc2] endPoints2_d2[n3][n1] = d2Bladder[nc2] - # for n1 in range(elementsCountAroundOstium): - # if n1 == 0: - # startDerivativesMap[0][n1] = ((-1, 0, 0), (1, 1, 0), None, (0, 1, 0)) - # startDerivativesMap[1][n1] = ((-1, 0, 0), (1, 1, 0), None, (0, 1, 0)) - # elif n1 == 1: - # startDerivativesMap[0][n1] = ((0, 1, 0), (1, 0, 0), None) - # startDerivativesMap[1][n1] = ((0, 1, 0), (1, 0, 0), None) - # elif n1 == 2: - # startDerivativesMap[0][n1] = ((0, 1, 0), (1, -1, 0), None, (1, 0, 0)) - # startDerivativesMap[1][n1] = ((0, 1, 0), (1, -1, 0), None, (1, 0, 0)) - # elif n1 == 3: - # startDerivativesMap[0][n1] = ((1, 0, 0), (0, -1, 0), None) - # startDerivativesMap[1][n1] = ((1, 0, 0), (0, -1, 0), None) - # elif n1 == 4: - # startDerivativesMap[0][n1] = ((1, 0, 0), (-1, -1, 0), None, (0, -1, 0)) - # startDerivativesMap[1][n1] = ((1, 0, 0), (-1, -1, 0), None, (0, -1, 0)) - # elif n1 == 5: - # startDerivativesMap[0][n1] = ((0, -1, 0), (-1, 0, 0), None) - # startDerivativesMap[1][n1] = ((0, -1, 0), (-1, 0, 0), None) - # elif n1 == 6: - # startDerivativesMap[0][n1] = ((0, -1, 0), (-1, 1, 0), None, (-1, 0, 0)) - # startDerivativesMap[1][n1] = ((0, -1, 0), (-1, 1, 0), None, (-1, 0, 0)) - # else: - # startDerivativesMap[0][n1] = ((-1, 0, 0), (0, 1, 0), None) - # startDerivativesMap[1][n1] = ((-1, 0, 0), (0, 1, 0), None) - for n1 in range(elementsCountAroundOstium): if n1 == 0: endDerivativesMap[0][n1] = ((-1, 0, 0), (-1, -1, 0), None, (0, 1, 0)) @@ -1004,7 +929,6 @@ 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(trackSurfaceOstium1.getProportion(o1_Positions[n])) @@ -1049,26 +973,15 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul 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=elementsCountAnnulusRadial, meshGroups=annulusMeshGroups, tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, rescaleEndDerivatives=True) - # print('Making test annulus') - # nodeIdentifier, elementIdentifier = createAnnulusMesh3d( - # nodes, mesh, nodeIdentifier, elementIdentifier, - # endPoints1_x, endPoints1_d1, endPoints1_d2, None, endNode1_Id, startDerivativesMap, - # o1_x, o1_d1, o1_d2, None, o1_NodeId, endDerivativesMapTest, - # elementsCountRadial=elementsCountAnnulusRadially, meshGroups=annulusMeshGroups, - # tracksurface=trackSurfaceOstium1, startProportions=endProportions1, endProportions=startProportions1, - # rescaleStartDerivatives=True, rescaleEndDerivatives = True, - # debugTest=True) - - # print('endNodeId = ', endNode2_Id) 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=elementsCountAnnulusRadial, meshGroups=annulusMeshGroups, tracksurface=trackSurfaceOstium2, startProportions=startProportions2, endProportions=endProportions2, rescaleEndDerivatives=True) @@ -1088,18 +1001,4 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, fm, nodes, mesh, ostiumDefaul mesh_destroy_elements_and_nodes_by_identifiers(mesh, element_identifiers) - # ####################################################### - # nodeIdentifier = 10000 - # for n in range(len(endPoints2_x[1])): - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, endPoints2_x[1][n]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, zero) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, zero) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, zero) - # print('NodeIdentifier = ', nodeIdentifier) #, xEnd[n]) #, d1[n], d2[n]) - # nodeIdentifier = nodeIdentifier + 1 - # fm.endChange() - # ############################################################ - return diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index b7bec09d..53527451 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -12,8 +12,6 @@ from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, setEftScaleFactorIds from scaffoldmaker.utils import interpolation as interp from scaffoldmaker.utils import vector -from scaffoldmaker.utils import matrix #KM -import math # KM def derivativeSignsToExpressionTerms(valueLabels, signs, scaleFactorIdx = None): @@ -37,7 +35,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, forceStartLinearXi3 = False, forceMidLinearXi3 = False, forceEndLinearXi3 = False, maxStartThickness = None, maxEndThickness = None, useCrossDerivatives = False, elementsCountRadial = 1, meshGroups = None, tracksurface = None, startProportions = None, endProportions = None, - rescaleStartDerivatives = False, rescaleEndDerivatives = False, debugTest = False): + rescaleStartDerivatives = False, rescaleEndDerivatives = False): """ 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. @@ -131,24 +129,6 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, cache = fm.createFieldcache() coordinates = findOrCreateFieldCoordinates(fm) - ############### MOVE BACK DOWN ################################################## - nodetemplate = nodes.createNodetemplate() - nodetemplate.defineField(coordinates) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) - if useCrossDerivatives: - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) - - # if debugTest: - # nodeIdentifier = 70000 - ##################################################################################### - # Build arrays of points from start to end px = [ [], [] ] pd1 = [ [], [] ] @@ -276,17 +256,6 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, md1 = interp.interpolateSampleLinear([ ad1, bd1 ], me, mxi) thi = interp.interpolateSampleLinear([ thicknesses[0][n1], thicknesses[-1][n1] ], me, mxi) - # if n1 == 0: - # if debugTest or debug2: - # for i in range(len(mx)): - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, mx[i]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, md1[i]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, md2[i]) - # print(nodeIdentifier, vector.normalise(md1[i])) - # nodeIdentifier += 1 - #md2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, fixStartDerivative = True, fixEndDerivative = True) for n2 in range(1, elementsCountRadial): @@ -316,8 +285,6 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for n1 in range(nodesCountAround): normal = vector.normalise(vector.crossproduct3(pd1[1][n2][n1], pd2[1][n2][n1])) - if debugTest: - normal = vector.normalise(vector.crossproduct3(pd2[1][n2][n1], pd1[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) ] @@ -403,18 +370,18 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # Create nodes ############## - # nodetemplate = nodes.createNodetemplate() - # nodetemplate.defineField(coordinates) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - # if useCrossDerivatives: - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) - # if useCrossDerivatives: - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) - # nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + if useCrossDerivatives: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS2DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D3_DS1DS2DS3, 1) nodetemplateLinearS3 = nodes.createNodetemplate() nodetemplateLinearS3.defineField(coordinates) nodetemplateLinearS3.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) @@ -443,8 +410,6 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, pd2[n3][n2][n1]) if pd3[n3][n2]: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, pd3[n3][n2][n1]) - # if nodeIdentifier == 402 or nodeIdentifier == 418: - # print('Node', nodeIdentifier, px[n3][n2][n1]) nodeIdentifier = nodeIdentifier + 1 nodeId[n3].append(rowNodeId) @@ -476,8 +441,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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] ] - # if debugTest: - # print('Element', elementIdentifier, nids) + if mapDerivatives: eft1 = eftFactory.createEftNoCrossDerivatives() setEftScaleFactorIds(eft1, [1], []) @@ -485,36 +449,13 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, eftFactory.setEftLinearDerivative2(eft1, [ 1, 5, 2, 6 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapStartDerivatives: if rescaleStartDerivatives: - if elementsCountRadial == 1 and rescaleEndDerivatives: - nodeScaleFactorOffset0 = e1 * 1000 - nodeScaleFactorOffset1 = en * 1000 - setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1, - nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) - # sfIDCheck = [1, - # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1, - # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1] - # sfCheck = [-1.0, - # scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], - # scaleFactorMapStart[0][en], scaleFactorMapStart[1][en], - # scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], - # scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]] - # if debugTest: - # print('scaleFactorIds') - # for i in range(len(sfIDCheck)): - # print(sfIDCheck[i]) - # print('scaleFactor') - # for i in range(len(sfCheck)): - # print(sfCheck[i]) - - else: - nodeScaleFactorOffset0 = e1 * 1000 - nodeScaleFactorOffset1 = en * 1000 - setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) + scaleFactorIds = [] + for i in range(2): + for n3 in range(2): + derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] + d2Map = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] + scaleFactorIds.append(getQuadrantID(d2Map)) + setEftScaleFactorIds(eft1, [1], scaleFactorIds) for i in range(2): lns = [ 1, 5 ] if (i == 0) else [ 2, 6 ] @@ -523,7 +464,7 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] # 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 = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] d3Map = derivativesMap[2] # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] @@ -532,25 +473,16 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) - if d2Map is None: - d2Map = (0, 1, 0) if d2Map is not None: - # if debugTest: - # # print('mapStartDerivatives - ', ln, 'remap ds2') - # print('rescaleStart', elementIdentifier, ln, sfIDCheck[so], sfCheck[so]) 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, so + 1 if rescaleStartDerivatives else None)) if d1Map is not None: - # if debugTest: - # print('mapStartDerivatives - ', ln, 'remap ds1') 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 debugTest: - # print('mapStartDerivatives - ', ln, 'remap ds3') 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: @@ -560,22 +492,14 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, if elementsCountRadial == 1 and rescaleStartDerivatives: pass else: - nodeScaleFactorOffset0 = e1 * 1000 - nodeScaleFactorOffset1 = en * 1000 - setEftScaleFactorIds(eft1, [1], [nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1]) - # sfIDCheck = [1, - # nodeScaleFactorOffset0 + 1, nodeScaleFactorOffset0 + 1, - # nodeScaleFactorOffset1 + 1, nodeScaleFactorOffset1 + 1] - # sfCheck = [-1.0, - # scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], - # scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]] - # print('rescaleEnd - scaleFactorIds') - # for i in range(len(sfIDCheck)): - # print(sfIDCheck[i]) - # print('scaleFactor') - # for i in range(len(sfCheck)): - # print(sfCheck[i]) + scaleFactorIds = [] + for i in range(2): + for n3 in range(2): + derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] + d2Map = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] + scaleFactorIds.append(getQuadrantID(d2Map)) + setEftScaleFactorIds(eft1, [1], scaleFactorIds) + for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] if elementsCountRadial == 1 and rescaleStartDerivatives: @@ -586,24 +510,17 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] # 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 = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] d3Map = derivativesMap[2] + # use temporary to safely swap DS1 and DS2: ln = [ lns[n3] ] so = sos + n3 - # if n3 == 0: - # print('i =', i, ', n3 =', n3, elementIdentifier, ln, so) if d1Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS1, [ ( Node.VALUE_LABEL_D2_DS1DS2, [] ) ]) if d3Map is not None: remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS3, [ ( Node.VALUE_LABEL_D2_DS2DS3, [] ) ]) - if d2Map is None: - d2Map = (0, 1, 0) if d2Map is not None: - # if debugTest: - # # print('mapEndDerivatives - ', ln, 'remap ds2') - # # print(elementIdentifier, nids[(ln[0] - 1)], sfIDCheck[so], sfCheck[so]) - # print('rescaleEnd', elementIdentifier, ln, sfIDCheck[so], sfCheck[so]) #, d2Map) # sfCheck[so], scaleFactorMap[n3][e1 if i == 0 else en]) remapEftNodeValueLabel(eft1, ln, Node.VALUE_LABEL_D_DS2, \ derivativeSignsToExpressionTerms((Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, @@ -611,13 +528,9 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, d2Map, so + 1 if rescaleEndDerivatives else None)) if d1Map is not None: - # if debugTest: - # print('mapEndDerivatives - ', ln, 'remap ds1') 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 debugTest: - # print('mapEndDerivatives - ', ln, 'remap ds3') 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)) @@ -637,18 +550,18 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, scaleFactorMapStart[0][en], scaleFactorMapStart[1][en], scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) - # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) elif rescaleStartDerivatives: result3 = element.setScaleFactors(eft1, [-1.0, scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], scaleFactorMapStart[0][en], scaleFactorMapStart[1][en]]) - # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) elif rescaleEndDerivatives: result3 = element.setScaleFactors(eft1, [ -1.0, scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) - # print('rescale - create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) + # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) else: result3 = element.setScaleFactors(eft1, [ -1.0 ]) @@ -660,17 +573,19 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, for meshGroup in rowMeshGroups[e2]: meshGroup.addElement(element) - # if rescaleEndDerivatives: - # nodeIdentifier = 90000 - # for n1 in range(len(endPointsx[1])): - # node = nodes.createNode(nodeIdentifier, nodetemplate) - # cache.setNode(node) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, endPointsx[1][n1]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, checkDerivativeStartList[n1]) - # coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, checkDerivativeEndList[n1]) - # print(nodeIdentifier) - # nodeIdentifier += 1 - 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)] + for i in range(len(maps)): + if d == maps[i]: + return i + 1 From 63f5626056392e5f7c7d7ba0af950cbe40f62213 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 12 Nov 2020 18:02:13 +1300 Subject: [PATCH 07/18] Modify central path for pig distal colon --- .../meshtypes/meshtype_3d_bladderurethra1.py | 1 + .../meshtypes/meshtype_3d_colon1.py | 79 ++++++++++--------- 2 files changed, 42 insertions(+), 38 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 55826edb..a774320d 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -25,6 +25,7 @@ from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node + class MeshType_3d_bladderurethra1(Scaffold_base): ''' Generates 3D bladder and urethra meshes with variable numbers diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index cef25e71..f95fe0e8 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -103,47 +103,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': { From 42e8c4d35ed4996046ef41db31ddb68f4b6198db Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 12 Nov 2020 18:49:16 +1300 Subject: [PATCH 08/18] Add parameter set for cow colon and colonsegment --- .../meshtypes/meshtype_3d_colon1.py | 85 ++++++++++++++++++- .../meshtypes/meshtype_3d_colonsegment1.py | 15 +++- 2 files changed, 97 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index f95fe0e8..b9a7906d 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 = { + 'Cow 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, @@ -172,6 +235,7 @@ def getName(): def getParameterSetNames(): return [ 'Default', + 'Cow 1', 'Human 1', 'Human 2', 'Mouse 1', @@ -183,6 +247,8 @@ def getParameterSetNames(): def getDefaultOptions(cls, parameterSetName='Default'): if 'Human 2' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 2'] + elif 'Cow 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Cow 1'] elif 'Mouse 1' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 1'] elif 'Mouse 2' in parameterSetName: @@ -193,7 +259,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 2'] else: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 1'] - if 'Mouse' in parameterSetName: + if 'Cow' in parameterSetName: + segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Cow 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') @@ -220,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 'Cow 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..b688d633 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', + 'Cow 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 'Cow' 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 From a0bab2751b516aededded89c58f58f298c17cdb3 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 12 Nov 2020 18:59:52 +1300 Subject: [PATCH 09/18] Update unit tests --- tests/test_bladder.py | 2 +- tests/test_cecum.py | 4 ++-- tests/test_colon.py | 14 +++++++------- tests/test_colonsegment.py | 2 +- 4 files changed, 11 insertions(+), 11 deletions(-) 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..b487cbb4 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, 65960.93476963606, 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, 127902.1860469739, delta=1.0E-6) if __name__ == "__main__": unittest.main() diff --git a/tests/test_colon.py b/tests/test_colon.py index b1afded5..66d403c2 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", "Cow 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..e916f33d 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", "Cow 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")) From be17b3816b7caa3b2299bfb4552fa80d462c58f7 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 17 Nov 2020 10:11:22 +1300 Subject: [PATCH 10/18] Change cow to cattle --- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 14 +++++++------- .../meshtypes/meshtype_3d_colonsegment1.py | 4 ++-- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index b9a7906d..6118a7df 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -27,7 +27,7 @@ class MeshType_3d_colon1(Scaffold_base): ''' centralPathDefaultScaffoldPackages = { - 'Cow 1' : ScaffoldPackage(MeshType_1d_path1, { + 'Cattle 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, 'D2 derivatives': True, @@ -235,7 +235,7 @@ def getName(): def getParameterSetNames(): return [ 'Default', - 'Cow 1', + 'Cattle 1', 'Human 1', 'Human 2', 'Mouse 1', @@ -247,8 +247,8 @@ def getParameterSetNames(): def getDefaultOptions(cls, parameterSetName='Default'): if 'Human 2' in parameterSetName: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 2'] - elif 'Cow 1' in parameterSetName: - centralPathOption = cls.centralPathDefaultScaffoldPackages['Cow 1'] + 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: @@ -259,8 +259,8 @@ def getDefaultOptions(cls, parameterSetName='Default'): centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 2'] else: centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 1'] - if 'Cow' in parameterSetName: - segmentProfileOption = ScaffoldPackage(MeshType_3d_colonsegment1, defaultParameterSetName = 'Cow 1') + 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: @@ -288,7 +288,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Refine number of elements along' : 1, 'Refine number of elements through wall' : 1 } - if 'Cow 1' in parameterSetName: + if 'Cattle 1' in parameterSetName: options['Number of segments'] = 40 options['Proximal length'] = 880.0 options['Transverse length'] = 3500.0 diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py index b688d633..04c05ebe 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colonsegment1.py @@ -43,7 +43,7 @@ def getName(): def getParameterSetNames(): return [ 'Default', - 'Cow 1', + 'Cattle 1', 'Human 1', 'Mouse 1', 'Pig 1'] @@ -79,7 +79,7 @@ def getDefaultOptions(parameterSetName='Default'): 'Refine number of elements along segment' : 1, 'Refine number of elements through wall' : 1 } - if 'Cow' 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 From 6a3b591eff5e4a7ac60567be621926f9a82cb127 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Tue, 17 Nov 2020 10:14:02 +1300 Subject: [PATCH 11/18] Update cow to cattle in unittests --- tests/test_colon.py | 2 +- tests/test_colonsegment.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_colon.py b/tests/test_colon.py index 66d403c2..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", "Cow 1", "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': { diff --git a/tests/test_colonsegment.py b/tests/test_colonsegment.py index e916f33d..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", "Cow 1", "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")) From c5acae3aca2225f7331ab58983d9487d8dc4deaa Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 18 Nov 2020 21:27:41 +1300 Subject: [PATCH 12/18] Improve radial rescaling of annulus mesh --- src/scaffoldmaker/utils/annulusmesh.py | 545 +++++++++++-------------- src/scaffoldmaker/utils/eft_utils.py | 6 +- 2 files changed, 251 insertions(+), 300 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 53527451..639a3971 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -29,13 +29,56 @@ def derivativeSignsToExpressionTerms(valueLabels, signs, scaleFactorIdx = None): 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, - rescaleStartDerivatives = False, rescaleEndDerivatives = False): + 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. @@ -78,6 +121,9 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, and radially adjacent nodes. :param rescaleEndDerivatives: If true, rescale end derivatives to arclength between endPoints and radially adjacent nodes. + :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' @@ -141,230 +187,141 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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] + 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) + + # 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) + + # set scalefactors if rescaling, make same on inside for now + if rescaleStartDerivatives: + scaleFactor = vector.magnitude(md2[0])/vector.magnitude(ad2) + for tn3 in range(2): + scaleFactorMapStart[tn3].append(scaleFactor) + if rescaleEndDerivatives: + scaleFactor = vector.magnitude(md2[-1])/vector.magnitude(bd2) + for tn3 in range(2): + scaleFactorMapEnd[tn3].append(scaleFactor) - 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: - arcLength = interp.computeCubicHermiteArcLength(startPointsx[1][n1], ad2, endPointsx[1][n1], bd2, - rescaleDerivatives=True) - ad2Scaled = vector.setMagnitude(ad2, arcLength / elementsCountRadial * 2 * (0.4 if vector.magnitude(ad2) < vector.magnitude(bd2) else 0.6)) - bd2Scaled = vector.setMagnitude(bd2, arcLength / elementsCountRadial * 2 * (0.6 if vector.magnitude(ad2) < vector.magnitude(bd2) else 0.4)) - - mx, md2, md1, md3, mProportions = \ - tracksurface.createHermiteCurvePoints(startProportions[n1][0], startProportions[n1][1], - endProportions[n1][0], endProportions[n1][1], - elementsCountRadial, derivativeStart=ad2Scaled, - derivativeEnd=bd2Scaled) - - mx, md2, md1 = tracksurface.resampleHermiteCurvePointsSmooth(mx, md2, md1, md3, mProportions, - derivativeMagnitudeStart= vector.magnitude(ad2Scaled), - derivativeMagnitudeEnd= vector.magnitude(bd2Scaled))[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) + for n2 in range(1, elementsCountRadial): + px [n3][n2][n1] = mx [n2] + pd1[n3][n2][n1] = md1[n2] + pd2[n3][n2][n1] = md2[n2] + thicknesses[n2][n1] = thi[n2] - #md2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, fixStartDerivative = True, fixEndDerivative = True) + # 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] + if rescaleStartDerivatives: + md2[0] = [ d*scaleFactorMapStart[n3][n1] for d in md2[0]] + if rescaleEndDerivatives: + md2[-1] = [ d*scaleFactorMapEnd[n3][n1] for d in md2[-1]] + sd2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, + fixAllDirections = True, fixStartDerivative = True, fixEndDerivative = True, + magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) 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] - - # 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] - - # Calculate arcLength to determine scale factor for remapped d2 - if rescaleStartDerivatives: - scaleFactorMapStart = [] - for n3 in range(2): - scaleFactorList = [] - for n1 in range(nodesCountAround): - startDerivativeMapD2 = (0, 1, 0) if startDerivativesMap[n3][n1][1] is None else startDerivativesMap[n3][n1][1] - v1 = px[n3][0][n1] - d1 = [startDerivativeMapD2[0] * pd1[n3][0][n1][c] + startDerivativeMapD2[1] * pd2[n3][0][n1][c] for c in range(3)] - v2 = px[n3][1][n1] - d2 = pd2[n3][1][n1] - if elementsCountRadial == 1 and endDerivativesMap is not None: - endDerivativeMapD2 = (0, 1, 0) if endDerivativesMap[n3][n1][1] is None else endDerivativesMap[n3][n1][1] - d2 = [endDerivativeMapD2[0] * pd1[n3][1][n1][c] + endDerivativeMapD2[1] * pd2[n3][1][n1][c] for c in range(3)] - arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) - unscaledArcLength = vector.magnitude([startDerivativeMapD2[0] * startPointsd1[n3][n1][c] + - startDerivativeMapD2[1] * startPointsd2[n3][n1][c] for c in - range(3)]) - scaleFactorNode = arcLength / unscaledArcLength - scaleFactorList.append(scaleFactorNode) - scaleFactorMapStart.append(scaleFactorList) - - if rescaleEndDerivatives: - scaleFactorMapEnd = [] - for n3 in range(2): - scaleFactorList = [] - for n1 in range(nodesCountAround): - endDerivativesMapD2 = (0, 1, 0) if endDerivativesMap[n3][n1][1] is None else endDerivativesMap[n3][n1][1] - v1 = px[n3][-2][n1] - d1 = pd2[n3][-2][n1] - if elementsCountRadial == 1 and startDerivativesMap is not None: - startDerivativeMapD2 = (0, 1, 0) if startDerivativesMap[n3][n1][1] is None else startDerivativesMap[n3][n1][1] - d1 = [startDerivativeMapD2[0] * pd1[n3][0][n1][c] + startDerivativeMapD2[1] * pd2[n3][0][n1][c] for c in range(3)] - v2 = px[n3][-1][n1] - d2 = [endDerivativesMapD2[0] * pd1[n3][-1][n1][c] + endDerivativesMapD2[1] * pd2[n3][-1][n1][c] for c in range(3)] - arcLength = interp.computeCubicHermiteArcLength(v1, d1, v2, d2, rescaleDerivatives=True) - unscaledArcLength = vector.magnitude([endDerivativesMapD2[0] * endPointsd1[n3][n1][c] + - endDerivativesMapD2[1] * endPointsd2[n3][n1][c] for c in range(3)]) - scaleFactorNode = arcLength / unscaledArcLength - scaleFactorList.append(scaleFactorNode) - scaleFactorMapEnd.append(scaleFactorList) + pd2[n3][n2][n1] = sd2[n2] ############## # Create nodes @@ -432,105 +389,120 @@ 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: - if rescaleStartDerivatives: - scaleFactorIds = [] - for i in range(2): - for n3 in range(2): - derivativesMap = startDerivativesMap[n3][e1] if (i == 0) else startDerivativesMap[n3][en] - d2Map = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] - scaleFactorIds.append(getQuadrantID(d2Map)) - setEftScaleFactorIds(eft1, [1], scaleFactorIds) - for i in range(2): lns = [ 1, 5 ] if (i == 0) else [ 2, 6 ] - sos = 1 if (i == 0) else 3 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 = (0, 1, 0) if derivativesMap[1] is None else 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] ] - so = sos + 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, - so + 1 if rescaleStartDerivatives else None)) - if d1Map is not None: + (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: eftFactory.setEftLinearDerivative2(eft1, [ 3, 7, 4, 8 ], Node.VALUE_LABEL_D_DS3, [ Node.VALUE_LABEL_D2_DS1DS3 ]) if mapEndDerivatives: - if rescaleEndDerivatives: - if elementsCountRadial == 1 and rescaleStartDerivatives: - pass - else: - scaleFactorIds = [] - for i in range(2): - for n3 in range(2): - derivativesMap = endDerivativesMap[n3][e1] if (i == 0) else endDerivativesMap[n3][en] - d2Map = (0, 1, 0) if derivativesMap[1] is None else derivativesMap[1] - scaleFactorIds.append(getQuadrantID(d2Map)) - setEftScaleFactorIds(eft1, [1], scaleFactorIds) - for i in range(2): lns = [ 3, 7 ] if (i == 0) else [ 4, 8 ] - if elementsCountRadial == 1 and rescaleStartDerivatives: - sos = 5 if (i == 0) else 7 - else: - sos = 1 if (i == 0) else 3 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 = (0, 1, 0) if derivativesMap[1] is None else 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] ] - so = sos + 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, so + 1 if rescaleEndDerivatives else None)) - - 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)) @@ -542,31 +514,9 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, element = mesh.createElement(elementIdentifier, elementtemplate1) result2 = element.setNodesByIdentifier(eft1, nids) - - if mapDerivatives: - if elementsCountRadial == 1 and rescaleStartDerivatives and rescaleEndDerivatives: - result3 = element.setScaleFactors(eft1, [-1.0, - scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], - scaleFactorMapStart[0][en], scaleFactorMapStart[1][en], - scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], - scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) - # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) - elif rescaleStartDerivatives: - result3 = element.setScaleFactors(eft1, [-1.0, - scaleFactorMapStart[0][e1], scaleFactorMapStart[1][e1], - scaleFactorMapStart[0][en], scaleFactorMapStart[1][en]]) - # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) - - elif rescaleEndDerivatives: - result3 = element.setScaleFactors(eft1, [ -1.0, - scaleFactorMapEnd[0][e1], scaleFactorMapEnd[1][e1], - scaleFactorMapEnd[0][en], scaleFactorMapEnd[1][en]]) - # print('create element annulus', element.isValid(), elementIdentifier, result2, result3, nids) - - else: - result3 = element.setScaleFactors(eft1, [ -1.0 ]) - # print('else 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: @@ -584,8 +534,9 @@ def getQuadrantID(d): :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)] + 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 From 598f1fb9fba97f1e6c37eb42df77c2fb1b1c1a43 Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Wed, 18 Nov 2020 21:35:27 +1300 Subject: [PATCH 13/18] Fix annulus mesh param documentation --- src/scaffoldmaker/utils/annulusmesh.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index 639a3971..ce5154b2 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -19,7 +19,7 @@ def derivativeSignsToExpressionTerms(valueLabels, signs, scaleFactorIdx = None): 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: List of 1, -1 (scale by scale factor indexed by scaleFactorIdx) 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)): @@ -117,10 +117,11 @@ 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: If true, rescale start derivatives to arclength between startPoints - and radially adjacent nodes. - :param rescaleEndDerivatives: If true, rescale end derivatives to arclength between endPoints - and radially adjacent nodes. + :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. From bf94ef80d9bcbc3b1656a6acd02b0f54ef0062df Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 19 Nov 2020 11:09:27 +1300 Subject: [PATCH 14/18] Update names to ureter position around and ureter position down --- .../meshtypes/meshtype_3d_bladderurethra1.py | 53 ++++++++++--------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index a774320d..ac1ffa2f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -193,9 +193,9 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Urethra wall thickness': 0.5, 'Include ureter': False, 'Ureter': copy.deepcopy(ostiumOption), - 'Ostium position around': 0.67, # should be on the dorsal part (> 0.5) - 'Ostium position down': 0.83, - 'Number of radial elements on annulus': 2, + '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.67 # should be on the dorsal part (> 0.5) - options['Ostium position down'] = 0.83 + 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,9 +236,9 @@ def getOrderedOptionNames(): 'Urethra wall thickness', 'Include ureter', 'Ureter', - 'Ostium position around', - 'Ostium position down', - 'Number of radial elements on annulus', + 'Ureter position around', + 'Ureter position down', + 'Number of elements ureter radial', 'Use cross derivatives', 'Use linear through wall', 'Refine', @@ -293,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']: @@ -301,14 +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['Ostium position around'] < 0.1: - options['Ostium position around'] = 0.1 - elif options['Ostium position around'] > 0.9: - options['Ostium position around'] = 0.9 - if options['Ostium position down'] < 0.15: - options['Ostium position down'] = 0.15 - elif options['Ostium position down'] > 0.95: - options['Ostium position down'] = 0.95 + 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): @@ -352,9 +353,9 @@ def generateBaseMesh(cls, region, options): ostiumOptions = options['Ureter'] ostiumDefaultOptions = ostiumOptions.getScaffoldSettings() elementsCountAroundOstium = ostiumDefaultOptions['Number of elements around ostium'] - elementsCountAnnulusRadial = options['Number of radial elements on annulus'] - ostiumPositionAround = options['Ostium position around'] - ostiumPositionDown = options['Ostium position down'] + elementsCountUreterRadial = options['Number of elements ureter radial'] + ureterPositionAround = options['Ureter position around'] + ureterPositionDown = options['Ureter position down'] firstNodeIdentifier = 1 firstElementIdentifier = 1 @@ -647,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 @@ -717,7 +718,7 @@ def generateBaseMesh(cls, region, options): trackSurfaceOstium1 = TrackSurface(elementsCount1, elementsCount2, nodesOnTrackSurface_x, nodesOnTrackSurface_d1, nodesOnTrackSurface_d2) - ostium1Position = trackSurfaceOstium1.createPositionProportion(ostiumPositionAround, ostiumPositionDown) + ostium1Position = trackSurfaceOstium1.createPositionProportion(ureterPositionAround, ureterPositionDown) ostiumElementPositionAround = ostium1Position.e1 ostiumElementPositionDown = ostium1Position.e2 @@ -747,7 +748,7 @@ def generateBaseMesh(cls, region, options): ostium1Position, trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, ostiumElementPositionAround, xFinal, d1Final, d2Final, nextNodeIdentifier, - nextElementIdentifier, elementsCountAnnulusRadial, + nextElementIdentifier, elementsCountUreterRadial, annulusMeshGroups) fm.endChange() @@ -825,7 +826,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt ostium1Position, trackSurfaceOstium2, ostium2Position, ostiumElementPositionDown, ostiumElementPositionAround, xBladder, d1Bladder, d2Bladder, - nextNodeIdentifier, nextElementIdentifier, elementsCountAnnulusRadial, + nextNodeIdentifier, nextElementIdentifier, elementsCountUreterRadial, annulusMeshGroups = []): # Create ureters on the surface @@ -974,7 +975,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt 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=elementsCountAnnulusRadial, meshGroups=annulusMeshGroups, + elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, rescaleEndDerivatives=True) @@ -982,7 +983,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt 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=elementsCountAnnulusRadial, meshGroups=annulusMeshGroups, + elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, tracksurface=trackSurfaceOstium2, startProportions=startProportions2, endProportions=endProportions2, rescaleEndDerivatives=True) From 3c92f767cdce6fc69b4ba9115c8175a423430dd9 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 19 Nov 2020 11:48:02 +1300 Subject: [PATCH 15/18] Rename ostium to ureter --- .../meshtypes/meshtype_3d_bladderurethra1.py | 136 +++++++++--------- 1 file changed, 68 insertions(+), 68 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index ac1ffa2f..0faa9066 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,7 +192,7 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Urethra diameter 2': 1.0, 'Urethra wall thickness': 0.5, 'Include ureter': False, - 'Ureter': copy.deepcopy(ostiumOption), + '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, @@ -317,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): @@ -350,9 +350,9 @@ 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'] + 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'] @@ -706,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 = [] @@ -715,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(ureterPositionAround, ureterPositionDown) - ostiumElementPositionAround = ostium1Position.e1 - ostiumElementPositionDown = ostium1Position.e2 + ureter1Position = trackSurfaceUreter1.createPositionProportion(ureterPositionAround, ureterPositionDown) + ureterElementPositionAround = ureter1Position.e1 + ureterElementPositionDown = ureter1Position.e2 - # Create trackSurface at the outer layer of the bladder for ostium 2 + # Create trackSurface at the outer layer of the bladder for ureter 2 nodesOnTrackSurface2_x = [] nodesOnTrackSurface2_d1 = [] nodesOnTrackSurface2_d2 = [] @@ -735,18 +735,18 @@ def generateBaseMesh(cls, region, options): 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//2 - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0 else 0), - ostiumElementPositionDown, - (1 - ostium1Position.xi1) if ostium1Position.xi1 > 0 else ostium1Position.xi1, - ostium1Position.xi2) + 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, nodes, mesh, ostiumDefaultOptions, - elementsCountAround, elementsCountThroughWall, elementsCountAroundOstium, - trackSurfaceOstium1, - ostium1Position, trackSurfaceOstium2, ostium2Position, - ostiumElementPositionDown, ostiumElementPositionAround, + generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ureterDefaultOptions, + elementsCountAround, elementsCountThroughWall, elementsCountAroundUreter, + trackSurfaceUreter1, + ureter1Position, trackSurfaceUreter2, ureter2Position, + ureterElementPositionDown, ureterElementPositionAround, xFinal, d1Final, d2Final, nextNodeIdentifier, nextElementIdentifier, elementsCountUreterRadial, annulusMeshGroups) @@ -821,55 +821,55 @@ 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, nodes, mesh, ostiumDefaultOptions,elementsCountAround, - elementsCountThroughWall, elementsCountAroundOstium, trackSurfaceOstium1, - ostium1Position, trackSurfaceOstium2, ostium2Position, - ostiumElementPositionDown, ostiumElementPositionAround, +def generateOstiumsAndAnnulusMeshOnBladder(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]) endPointStartId1 = elementsCountThroughWall + 1 \ - + (elementsCountThroughWall + 1) * elementsCountAround * (ostiumElementPositionDown - (1 if ostium1Position.xi2 > 0.5 else 2)) \ - + ostiumElementPositionAround + (1 if ostium1Position.xi1 > 0.5 else 0) + + (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]) endPointStartId2 = elementsCountThroughWall + 1 \ + (elementsCountThroughWall + 1) * elementsCountAround * \ - (ostiumElementPositionDown - (1 if ostium1Position.xi2 > 0.5 else 2)) \ - + elementsCountAround - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0.5 else 0) + (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] + # 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): @@ -895,7 +895,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt 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] @@ -905,7 +905,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt 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)) @@ -933,23 +933,23 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt startProportions1 = [] for n in range(len(o1_Positions)): - startProportions1.append(trackSurfaceOstium1.getProportion(o1_Positions[n])) + startProportions1.append(trackSurfaceUreter1.getProportion(o1_Positions[n])) startProportions2 = [] for n in range(len(o2_Positions)): - startProportions2.append(trackSurfaceOstium2.getProportion(o2_Positions[n])) + startProportions2.append(trackSurfaceUreter2.getProportion(o2_Positions[n])) endProportions1 = [] - elementsAroundTrackSurface1 = trackSurfaceOstium1.elementsCount1 - elementsAlongTrackSurface1 = trackSurfaceOstium1.elementsCount2 + elementsAroundTrackSurface1 = trackSurfaceUreter1.elementsCount1 + elementsAlongTrackSurface1 = trackSurfaceUreter1.elementsCount2 endProportions2 = [] - elementsAroundTrackSurface2 = trackSurfaceOstium2.elementsCount1 - elementsAlongTrackSurface2 = trackSurfaceOstium2.elementsCount2 + elementsAroundTrackSurface2 = trackSurfaceUreter2.elementsCount1 + elementsAlongTrackSurface2 = trackSurfaceUreter2.elementsCount2 - firstIdxAround1 = ostiumElementPositionAround + (0 if ostium1Position.xi1 > 0.5 else -1) - firstIdxAlong = ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1) - firstIdxAround2 = elementsCountAround//2 - ostiumElementPositionAround - (2 if ostium1Position.xi1 > 0.5 else 1) + 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, @@ -976,7 +976,7 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt o1_x, o1_d1, o1_d2, None, o1_NodeId, None, endPoints1_x, endPoints1_d1, endPoints1_d2, None, endNode1_Id, endDerivativesMap, elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, - tracksurface=trackSurfaceOstium1, startProportions=startProportions1, endProportions=endProportions1, + tracksurface=trackSurfaceUreter1, startProportions=startProportions1, endProportions=endProportions1, rescaleEndDerivatives=True) nodeIdentifier, elementIdentifier = createAnnulusMesh3d( @@ -984,15 +984,15 @@ def generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ostiumDefaultOpt o2_x, o2_d1, o2_d2, None, o2_NodeId, None, endPoints2_x, endPoints2_d1, endPoints2_d2, None, endNode2_Id, endDerivativesMap, elementsCountRadial=elementsCountUreterRadial, meshGroups=annulusMeshGroups, - tracksurface=trackSurfaceOstium2, startProportions=startProportions2, endProportions=endProportions2, + tracksurface=trackSurfaceUreter2, startProportions=startProportions2, endProportions=endProportions2, rescaleEndDerivatives=True) # Delete elements under annulus mesh element_identifiers = [] - elementToDeleteStartIdx1 = elementsCountThroughWall * elementsCountAround * (ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1)) \ - + ostiumElementPositionAround + (1 if ostium1Position.xi1 > 0.5 else 0) - elementToDeleteStartIdx2 = elementsCountThroughWall * elementsCountAround * (ostiumElementPositionDown - (0 if ostium1Position.xi2 > 0.5 else 1)) \ - + elementsCountAround - ostiumElementPositionAround + (-1 if ostium1Position.xi1 > 0.5 else 0) + 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] From 1f8bb2adb5a0ea8357b9288a1e3819f829693d58 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 19 Nov 2020 11:54:00 +1300 Subject: [PATCH 16/18] Rename function to generateUreters --- .../meshtypes/meshtype_3d_bladderurethra1.py | 25 ++++++++----------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py index 0faa9066..fa68b8a0 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py @@ -742,14 +742,11 @@ def generateBaseMesh(cls, region, options): (1 - ureter1Position.xi1) if ureter1Position.xi1 > 0 else ureter1Position.xi1, ureter1Position.xi2) annulusMeshGroups = [neckMeshGroup, urinaryBladderMeshGroup] - generateOstiumsAndAnnulusMeshOnBladder(region, nodes, mesh, ureterDefaultOptions, - elementsCountAround, elementsCountThroughWall, elementsCountAroundUreter, - trackSurfaceUreter1, - ureter1Position, trackSurfaceUreter2, ureter2Position, - ureterElementPositionDown, ureterElementPositionAround, - xFinal, d1Final, d2Final, nextNodeIdentifier, - nextElementIdentifier, elementsCountUreterRadial, - 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 @@ -821,13 +818,11 @@ 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, nodes, mesh, ureterDefaultOptions,elementsCountAround, - elementsCountThroughWall, elementsCountAroundUreter, trackSurfaceUreter1, - ureter1Position, trackSurfaceUreter2, ureter2Position, - ureterElementPositionDown, ureterElementPositionAround, - xBladder, d1Bladder, d2Bladder, - nextNodeIdentifier, nextElementIdentifier, elementsCountUreterRadial, - 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 From 773643d24a85d4157bba566ea4013e83f2653951 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 19 Nov 2020 12:29:13 +1300 Subject: [PATCH 17/18] Modify to calculate scale factor for inside --- src/scaffoldmaker/utils/annulusmesh.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/scaffoldmaker/utils/annulusmesh.py b/src/scaffoldmaker/utils/annulusmesh.py index ce5154b2..f06814c2 100644 --- a/src/scaffoldmaker/utils/annulusmesh.py +++ b/src/scaffoldmaker/utils/annulusmesh.py @@ -256,12 +256,10 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, # set scalefactors if rescaling, make same on inside for now if rescaleStartDerivatives: scaleFactor = vector.magnitude(md2[0])/vector.magnitude(ad2) - for tn3 in range(2): - scaleFactorMapStart[tn3].append(scaleFactor) + scaleFactorMapStart[n3].append(scaleFactor) if rescaleEndDerivatives: scaleFactor = vector.magnitude(md2[-1])/vector.magnitude(bd2) - for tn3 in range(2): - scaleFactorMapEnd[tn3].append(scaleFactor) + scaleFactorMapEnd[n3].append(scaleFactor) for n2 in range(1, elementsCountRadial): px [n3][n2][n1] = mx [n2] @@ -314,13 +312,20 @@ def createAnnulusMesh3d(nodes, mesh, nextNodeIdentifier, nextElementIdentifier, 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: - md2[0] = [ d*scaleFactorMapStart[n3][n1] for d in md2[0]] + scaleFactor = vector.magnitude(sd2[0]) / vector.magnitude(md2[0]) + scaleFactorMapStart[n3].append(scaleFactor) if rescaleEndDerivatives: - md2[-1] = [ d*scaleFactorMapEnd[n3][n1] for d in md2[-1]] - sd2 = interp.smoothCubicHermiteDerivativesLine(mx, md2, - fixAllDirections = True, fixStartDerivative = True, fixEndDerivative = True, - magnitudeScalingMode = interp.DerivativeScalingMode.HARMONIC_MEAN) + 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] From 678a205a4cf2941949fac9e6e03f3e28cee93539 Mon Sep 17 00:00:00 2001 From: Mabelle Lin Date: Thu, 19 Nov 2020 12:35:58 +1300 Subject: [PATCH 18/18] Update unit test --- tests/test_cecum.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_cecum.py b/tests/test_cecum.py index b487cbb4..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.93476963606, 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, 127902.1860469739, delta=1.0E-6) + self.assertAlmostEqual(volume, 127894.8597203706, delta=1.0E-6) if __name__ == "__main__": unittest.main()