From 81190129465930b7f56dda3ce1778353a6f02139 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 17 Sep 2020 09:52:59 +1200 Subject: [PATCH 01/10] Check if two vectors are parallel. --- src/scaffoldmaker/utils/vector.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py index 996ee833..7401b9fa 100644 --- a/src/scaffoldmaker/utils/vector.py +++ b/src/scaffoldmaker/utils/vector.py @@ -42,3 +42,15 @@ def addVectors(v1,v2,s1=1.0,s2=1.0): :return: Vector s1*v1+s2*v2 ''' return [(s1 * v1[c] + s2 * v2[c]) for c in range(len(v1))] + + +def isVectorsParallel(v1, v2): + """ + :return: True if the vectors are parallel. + """ + assert (len(v2) == len(v1)), 'Vectors lengths are not the same.' + s = magnitude(v2)/magnitude(v1) + for c in range(len(v2)): + if v2[c] != s * v1[c]: + return False + return True From 45d7105e29b95174d00d68bc7d4bdfa1e7130c26 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 17 Sep 2020 10:00:40 +1200 Subject: [PATCH 02/10] Add a boolean flag to turn on d2 and d3 --- .../meshtypes/meshtype_1d_path1.py | 35 +++++++++++++++---- 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 4198ae4b..8bad298a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -24,6 +24,8 @@ def getName(): def getDefaultOptions(parameterSetName='Default'): return { 'Coordinate dimensions' : 3, + 'D2 derivative': False, + 'D3 derivative': False, 'Length' : 1.0, 'Number of elements' : 1 } @@ -32,18 +34,25 @@ def getDefaultOptions(parameterSetName='Default'): def getOrderedOptionNames(): return [ 'Coordinate dimensions', + 'D2 derivative', + 'D3 derivative', 'Length', 'Number of elements' ] @staticmethod def checkOptions(options): + dependentChanges = False if (options['Coordinate dimensions'] < 1) : options['Coordinate dimensions'] = 1 elif (options['Coordinate dimensions'] > 3) : options['Coordinate dimensions'] = 3 if (options['Number of elements'] < 1) : options['Number of elements'] = 1 + if options['D3 derivative']: + options['D2 derivative'] = True + dependentChanges = True + return dependentChanges @classmethod def generateBaseMesh(cls, region, options): @@ -53,6 +62,8 @@ def generateBaseMesh(cls, region, options): :return: [] empty list of AnnotationGroup """ coordinateDimensions = options['Coordinate dimensions'] + d2derivative = options['D2 derivative'] + d3derivative = options['D3 derivative'] length = options['Length'] elementsCount = options['Number of elements'] @@ -69,22 +80,34 @@ def generateBaseMesh(cls, region, options): nodetemplate.defineField(coordinates) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) - nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + if d2derivative: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) + if d3derivative: + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) nodeIdentifier = 1 x = [ 0.0, 0.0, 0.0 ] dx_ds1 = [ length/elementsCount, 0.0, 0.0 ] - dx_ds2 = [ 0.0, 1.0, 0.0 ] - d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ] + if d2derivative: + dx_ds2 = [ 0.0, 1.0, 0.0 ] + d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ] + if d3derivative: + dx_ds3 = [ 0.0, 0.0, 1.0 ] + d2x_ds1ds3 = [ 0.0, 0.0, 0.0 ] for n in range(elementsCount + 1): x[0] = length*n/elementsCount node = nodes.createNode(nodeIdentifier, nodetemplate) cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) - coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2) + if d2derivative: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2) + if d3derivative: + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d2x_ds1ds3) nodeIdentifier = nodeIdentifier + 1 ################# From 0b43143af26f09769e297643b98951fd2f96860a Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 17 Sep 2020 18:02:17 +1200 Subject: [PATCH 03/10] Add functions to find vector projection and rejection --- src/scaffoldmaker/utils/vector.py | 35 +++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py index 7401b9fa..ec46b2ce 100644 --- a/src/scaffoldmaker/utils/vector.py +++ b/src/scaffoldmaker/utils/vector.py @@ -44,6 +44,41 @@ def addVectors(v1,v2,s1=1.0,s2=1.0): return [(s1 * v1[c] + s2 * v2[c]) for c in range(len(v1))] +def scalarProjectionOfV1OnV2(v1, v2): + """ + :return: Scalar projection of v1 onto v2. + """ + return dotproduct(v1, normalise(v2)) + + +def vectorProjectionOfV1OnV2(v1, v2): + """ + Calculate vector projection of v1 on v2 + :return: A projection vector. + """ + s1 = scalarProjectionOfV1OnV2(v1, v2) + return scalarProduct(s1, normalise(v2)) + + +def vectorRejectionOfV1OnV2(v1, v2): + """ + Calculate vector rejection of v1 on v2 + :return: A rejection vector. + """ + v1p = vectorProjectionOfV1OnV2(v1, v2) + return addVectors(v1, v1p, 1.0, -1.0) + + +def scalarProduct(s, v): + """ + Calculate s * v + :param s: Scalar. + :param v: Vector. + :return: + """ + return [s * v[c] for c in range(len(v))] + + def isVectorsParallel(v1, v2): """ :return: True if the vectors are parallel. From 28c0671f58befe9a3762f7f51769b5f4f58ec29f Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 17 Sep 2020 18:10:10 +1200 Subject: [PATCH 04/10] Make d2 and d3 normal to d1 --- .../meshtypes/meshtype_1d_path1.py | 53 +++++++++++++++---- 1 file changed, 44 insertions(+), 9 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 8bad298a..4f645179 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -11,6 +11,9 @@ from opencmiss.zinc.node import Node from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.utils.interpolation import DerivativeScalingMode, smoothCubicHermiteDerivativesLine +from scaffoldmaker.utils import vector +from opencmiss.zinc.result import RESULT_OK + class MeshType_1d_path1(Scaffold_base): ''' @@ -25,7 +28,7 @@ def getDefaultOptions(parameterSetName='Default'): return { 'Coordinate dimensions' : 3, 'D2 derivative': False, - 'D3 derivative': False, + 'D3 and D2 derivatives': False, 'Length' : 1.0, 'Number of elements' : 1 } @@ -35,7 +38,7 @@ def getOrderedOptionNames(): return [ 'Coordinate dimensions', 'D2 derivative', - 'D3 derivative', + 'D3 and D2 derivatives', 'Length', 'Number of elements' ] @@ -49,7 +52,7 @@ def checkOptions(options): options['Coordinate dimensions'] = 3 if (options['Number of elements'] < 1) : options['Number of elements'] = 1 - if options['D3 derivative']: + if options['D3 and D2 derivatives']: options['D2 derivative'] = True dependentChanges = True return dependentChanges @@ -63,7 +66,7 @@ def generateBaseMesh(cls, region, options): """ coordinateDimensions = options['Coordinate dimensions'] d2derivative = options['D2 derivative'] - d3derivative = options['D3 derivative'] + d3d2derivatives = options['D3 and D2 derivatives'] length = options['Length'] elementsCount = options['Number of elements'] @@ -83,7 +86,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - if d3derivative: + if d3d2derivatives: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) @@ -93,7 +96,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: dx_ds2 = [ 0.0, 1.0, 0.0 ] d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ] - if d3derivative: + if d3d2derivatives: dx_ds3 = [ 0.0, 0.0, 1.0 ] d2x_ds1ds3 = [ 0.0, 0.0, 0.0 ] for n in range(elementsCount + 1): @@ -105,7 +108,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2) - if d3derivative: + if d3d2derivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d2x_ds1ds3) nodeIdentifier = nodeIdentifier + 1 @@ -135,6 +138,26 @@ def smoothPath(cls, region, options, mode : DerivativeScalingMode): d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode) setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ]) + @classmethod + def makeD2NoramlToD1(cls, region, options): + if not options['D2 derivative']: + return + x, d1, d2 = extractPathParametersFromRegion(region)[0:3] + for c in range(len(d1)): + td2 = vector.vectorRejectionOfV1OnV2(d2[c], d1[c]) + d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c])) + setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2]) + + @classmethod + def makeD3ND2NorlamToD1(cls, region, options): + if not options['D3 and D2 derivatives']: + return + cls.makeD2NoramlToD1(region, options) + x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) + for c in range(len(d1)): + d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) + setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) + @classmethod def getInteractiveFunctions(cls): """ @@ -142,7 +165,10 @@ def getInteractiveFunctions(cls): """ return [ ("Smooth d1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)), - ("Smooth d1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)) ] + ("Smooth d1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)), + ("Make D2 normal to D1", lambda region, options: cls.makeD2NoramlToD1(region, options)), + ("Make D3 and D2 normal to D1", lambda region, options: cls.makeD3ND2NorlamToD1(region, options)) + ] def extractPathParametersFromRegion(region, groupName=None): @@ -163,6 +189,8 @@ def extractPathParametersFromRegion(region, groupName=None): cd1 = [] cd2 = [] cd12 = [] + cd3 = [] + cd13 = [] nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) if groupName: group = fieldmodule.findFieldByName(groupName).castGroup() @@ -179,17 +207,24 @@ def extractPathParametersFromRegion(region, groupName=None): result, d1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, componentsCount) result, d2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, componentsCount) result, d12 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, componentsCount) + result, d3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, componentsCount) + result, d13 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) for c in range(componentsCount, 3): x.append(0.0) d1.append(0.0) d2.append(0.0) d12.append(0.0) + d3.append(0.0) + d13.append(0.0) cx.append(x) cd1.append(d1) cd2.append(d2) cd12.append(d12) + cd3.append(d3) + cd13.append(d13) node = nodeIter.next() - return cx, cd1, cd2, cd12 + + return cx, cd1, cd2, cd12, cd3, cd13 def setPathParameters(region, nodeValueLabels, nodeValues): From 8ac30734f9c5f718ede2a9da8a24c1d465557279 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 17 Sep 2020 18:21:39 +1200 Subject: [PATCH 05/10] Update colon, intestine and cecum central path --- src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py | 3 ++- src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py | 8 +++++++- .../meshtypes/meshtype_3d_smallintestine1.py | 3 ++- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index 0be3df24..f8f4b2c3 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -41,6 +41,7 @@ class MeshType_3d_cecum1(Scaffold_base): 'Pig 1': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, + 'D2 derivative': True, 'Length': 120.0, 'Number of elements': 3 }, @@ -307,7 +308,7 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion) + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') del tmpRegion diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 96babcc2..065463ef 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -30,6 +30,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Human 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, + 'D2 derivative': True, 'Length' : 1.0, 'Number of elements' : 8 }, @@ -48,6 +49,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Human 2' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, + 'D2 derivative': True, 'Length' : 1.0, 'Number of elements' : 8 }, @@ -66,6 +68,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, + 'D2 derivative': True, 'Length' : 1.0, 'Number of elements' : 7 }, @@ -83,6 +86,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Mouse 2' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, + 'D2 derivative': True, 'Length' : 1.0, 'Number of elements' : 4 }, @@ -97,6 +101,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Pig 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, + 'D2 derivative': True, 'Length' : 1.0, 'Number of elements' : 36 }, @@ -143,6 +148,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Pig 2': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, + 'D2 derivative': True, 'Length': 90.0, 'Number of elements': 3 }, @@ -391,7 +397,7 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion) + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') del tmpRegion diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index f6a7c529..19b26a07 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -29,6 +29,7 @@ class MeshType_3d_smallintestine1(Scaffold_base): centralPathDefaultScaffoldPackages = { 'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { + 'D2 derivative': True, 'Coordinate dimensions' : 3, 'Length' : 1.0, 'Number of elements' : 45 @@ -231,7 +232,7 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion) + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i],',', cd12[i], '],') del tmpRegion From 43d03958d9453b0a4f67e93dbab736183c31a91d Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Mon, 21 Sep 2020 14:57:37 +1200 Subject: [PATCH 06/10] D3 flag instead of D3D2 flag --- .../meshtypes/meshtype_1d_path1.py | 43 +++++++++++-------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 4f645179..c4f4899a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -28,7 +28,7 @@ def getDefaultOptions(parameterSetName='Default'): return { 'Coordinate dimensions' : 3, 'D2 derivative': False, - 'D3 and D2 derivatives': False, + 'D3 derivative': False, 'Length' : 1.0, 'Number of elements' : 1 } @@ -38,7 +38,7 @@ def getOrderedOptionNames(): return [ 'Coordinate dimensions', 'D2 derivative', - 'D3 and D2 derivatives', + 'D3 derivative', 'Length', 'Number of elements' ] @@ -52,9 +52,7 @@ def checkOptions(options): options['Coordinate dimensions'] = 3 if (options['Number of elements'] < 1) : options['Number of elements'] = 1 - if options['D3 and D2 derivatives']: - options['D2 derivative'] = True - dependentChanges = True + return dependentChanges @classmethod @@ -66,7 +64,7 @@ def generateBaseMesh(cls, region, options): """ coordinateDimensions = options['Coordinate dimensions'] d2derivative = options['D2 derivative'] - d3d2derivatives = options['D3 and D2 derivatives'] + d3derivative = options['D3 derivative'] length = options['Length'] elementsCount = options['Number of elements'] @@ -86,7 +84,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - if d3d2derivatives: + if d3derivative: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) @@ -96,7 +94,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: dx_ds2 = [ 0.0, 1.0, 0.0 ] d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ] - if d3d2derivatives: + if d3derivative: dx_ds3 = [ 0.0, 0.0, 1.0 ] d2x_ds1ds3 = [ 0.0, 0.0, 0.0 ] for n in range(elementsCount + 1): @@ -108,7 +106,7 @@ def generateBaseMesh(cls, region, options): if d2derivative: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2) - if d3d2derivatives: + if d3derivative: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d2x_ds1ds3) nodeIdentifier = nodeIdentifier + 1 @@ -139,7 +137,7 @@ def smoothPath(cls, region, options, mode : DerivativeScalingMode): setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ]) @classmethod - def makeD2NoramlToD1(cls, region, options): + def makeD2Normal(cls, region, options): if not options['D2 derivative']: return x, d1, d2 = extractPathParametersFromRegion(region)[0:3] @@ -149,14 +147,21 @@ def makeD2NoramlToD1(cls, region, options): setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2]) @classmethod - def makeD3ND2NorlamToD1(cls, region, options): - if not options['D3 and D2 derivatives']: + def makeD3Normal(cls, region, options): + if not options['D3 derivative']: return - cls.makeD2NoramlToD1(region, options) - x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) - for c in range(len(d1)): - d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) - setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) + if options['D2 derivative']: + x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) + for c in range(len(d1)): + d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) + setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) + else: + x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) + for c in range(len(d1)): + td3 = vector.vectorRejectionOfV1OnV2(d3[c], d1[c]) + d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c])) + setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) + @classmethod def getInteractiveFunctions(cls): @@ -166,8 +171,8 @@ def getInteractiveFunctions(cls): return [ ("Smooth d1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)), ("Smooth d1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)), - ("Make D2 normal to D1", lambda region, options: cls.makeD2NoramlToD1(region, options)), - ("Make D3 and D2 normal to D1", lambda region, options: cls.makeD3ND2NorlamToD1(region, options)) + ("Make D2 normal", lambda region, options: cls.makeD2Normal(region, options)), + ("Make D3 normal", lambda region, options: cls.makeD3Normal(region, options)) ] From a4d0ec663f296be8b4318ac71126c486fe229259 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Wed, 23 Sep 2020 13:30:44 +1200 Subject: [PATCH 07/10] Minor changes. Pass valueLabels as inpur argument. --- .../meshtypes/meshtype_1d_path1.py | 68 ++++++++++++------- .../meshtypes/meshtype_3d_cecum1.py | 6 +- .../meshtypes/meshtype_3d_colon1.py | 16 +++-- .../meshtypes/meshtype_3d_smallintestine1.py | 6 +- src/scaffoldmaker/utils/vector.py | 21 +++--- 5 files changed, 70 insertions(+), 47 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index c4f4899a..ea25b422 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -27,8 +27,8 @@ def getName(): def getDefaultOptions(parameterSetName='Default'): return { 'Coordinate dimensions' : 3, - 'D2 derivative': False, - 'D3 derivative': False, + 'D2 derivatives': False, + 'D3 derivatives': False, 'Length' : 1.0, 'Number of elements' : 1 } @@ -37,8 +37,8 @@ def getDefaultOptions(parameterSetName='Default'): def getOrderedOptionNames(): return [ 'Coordinate dimensions', - 'D2 derivative', - 'D3 derivative', + 'D2 derivatives', + 'D3 derivatives', 'Length', 'Number of elements' ] @@ -63,8 +63,8 @@ def generateBaseMesh(cls, region, options): :return: [] empty list of AnnotationGroup """ coordinateDimensions = options['Coordinate dimensions'] - d2derivative = options['D2 derivative'] - d3derivative = options['D3 derivative'] + d2derivatives = options['D2 derivatives'] + d3derivatives = options['D3 derivatives'] length = options['Length'] elementsCount = options['Number of elements'] @@ -81,20 +81,20 @@ def generateBaseMesh(cls, region, options): nodetemplate.defineField(coordinates) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) - if d2derivative: + if d2derivatives: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1) - if d3derivative: + if d3derivatives: nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1) nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1) nodeIdentifier = 1 x = [ 0.0, 0.0, 0.0 ] dx_ds1 = [ length/elementsCount, 0.0, 0.0 ] - if d2derivative: + if d2derivatives: dx_ds2 = [ 0.0, 1.0, 0.0 ] d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ] - if d3derivative: + if d3derivatives: dx_ds3 = [ 0.0, 0.0, 1.0 ] d2x_ds1ds3 = [ 0.0, 0.0, 0.0 ] for n in range(elementsCount + 1): @@ -103,10 +103,10 @@ def generateBaseMesh(cls, region, options): cache.setNode(node) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1) - if d2derivative: + if d2derivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2) - if d3derivative: + if d3derivatives: coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3) coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d2x_ds1ds3) nodeIdentifier = nodeIdentifier + 1 @@ -132,33 +132,34 @@ def generateBaseMesh(cls, region, options): @classmethod def smoothPath(cls, region, options, mode : DerivativeScalingMode): - x, d1 = extractPathParametersFromRegion(region)[0:2] + x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1]) d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode) setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ]) @classmethod def makeD2Normal(cls, region, options): - if not options['D2 derivative']: + if not options['D2 derivatives']: return - x, d1, d2 = extractPathParametersFromRegion(region)[0:3] + d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2]) for c in range(len(d1)): - td2 = vector.vectorRejectionOfV1OnV2(d2[c], d1[c]) + td2 = vector.vectorRejection(d2[c], d1[c]) d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c])) setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2]) @classmethod def makeD3Normal(cls, region, options): - if not options['D3 derivative']: + if not options['D3 derivatives']: return - if options['D2 derivative']: - x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) + if options['D2 derivatives']: + d1, d2, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, + Node.VALUE_LABEL_D_DS3]) for c in range(len(d1)): d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) else: - x, d1, d2, d12, d3, d13 = extractPathParametersFromRegion(region) + d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3]) for c in range(len(d1)): - td3 = vector.vectorRejectionOfV1OnV2(d3[c], d1[c]) + td3 = vector.vectorRejection(d3[c], d1[c]) d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c])) setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) @@ -169,19 +170,20 @@ def getInteractiveFunctions(cls): Supply client with functions for smoothing path parameters. """ return [ - ("Smooth d1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)), - ("Smooth d1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)), + ("Smooth D1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)), + ("Smooth D1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)), ("Make D2 normal", lambda region, options: cls.makeD2Normal(region, options)), ("Make D3 normal", lambda region, options: cls.makeD3Normal(region, options)) ] -def extractPathParametersFromRegion(region, groupName=None): +def extractPathParametersFromRegion(region, valueLabels, groupName=None): ''' Returns parameters of all nodes in region in identifier order. Assumes nodes in region have field coordinates (1 to 3 components). Currently limited to nodes with exactly value, d_ds1, d_ds2, d2_ds12, same as path 1 scaffold. + :param valueLabels: List of parameters required as list of node value labels. e.g. [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1]. :param groupName: Optional name of Zinc group to get parameters from. :return: cx, cd1, cd2, cd12 (all padded with zeroes to 3 components) ''' @@ -190,6 +192,7 @@ def extractPathParametersFromRegion(region, groupName=None): componentsCount = coordinates.getNumberOfComponents() assert componentsCount in [ 1, 2, 3 ], 'extractPathParametersFromRegion. Invalid coordinates number of components' cache = fieldmodule.createFieldcache() + cx = [] cd1 = [] cd2 = [] @@ -229,7 +232,22 @@ def extractPathParametersFromRegion(region, groupName=None): cd13.append(d13) node = nodeIter.next() - return cx, cd1, cd2, cd12, cd3, cd13 + result = () + for label in valueLabels: + if label == Node.VALUE_LABEL_VALUE: + result += (cx,) + if label == Node.VALUE_LABEL_D_DS1: + result += (cd1,) + if label == Node.VALUE_LABEL_D_DS2: + result += (cd2,) + if label == Node.VALUE_LABEL_D_DS3: + result += (cd3,) + if label == Node.VALUE_LABEL_D2_DS1DS2: + result += (cd12,) + if label == Node.VALUE_LABEL_D2_DS1DS3: + result += (cd13,) + + return result def setPathParameters(region, nodeValueLabels, nodeValues): diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py index f8f4b2c3..55fe8606 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py @@ -41,7 +41,7 @@ class MeshType_3d_cecum1(Scaffold_base): 'Pig 1': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length': 120.0, 'Number of elements': 3 }, @@ -308,7 +308,9 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion, + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2]) # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') del tmpRegion diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py index 065463ef..cef25e71 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py @@ -30,7 +30,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Human 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length' : 1.0, 'Number of elements' : 8 }, @@ -49,7 +49,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Human 2' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length' : 1.0, 'Number of elements' : 8 }, @@ -68,7 +68,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length' : 1.0, 'Number of elements' : 7 }, @@ -86,7 +86,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Mouse 2' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length' : 1.0, 'Number of elements' : 4 }, @@ -101,7 +101,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Pig 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { 'Coordinate dimensions' : 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length' : 1.0, 'Number of elements' : 36 }, @@ -148,7 +148,7 @@ class MeshType_3d_colon1(Scaffold_base): 'Pig 2': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, - 'D2 derivative': True, + 'D2 derivatives': True, 'Length': 90.0, 'Number of elements': 3 }, @@ -397,7 +397,9 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion, + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2]) # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],') del tmpRegion diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py index 19b26a07..5d0695c4 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py @@ -29,7 +29,7 @@ class MeshType_3d_smallintestine1(Scaffold_base): centralPathDefaultScaffoldPackages = { 'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings' : { - 'D2 derivative': True, + 'D2 derivatives': True, 'Coordinate dimensions' : 3, 'Length' : 1.0, 'Number of elements' : 45 @@ -232,7 +232,9 @@ def generateBaseMesh(cls, region, options): # Central path tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)[0:4] + cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion, + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, + Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2]) # for i in range(len(cx)): # print(i, '[', cx[i], ',', cd1[i], ',', cd2[i],',', cd12[i], '],') del tmpRegion diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py index ec46b2ce..da4f27f9 100644 --- a/src/scaffoldmaker/utils/vector.py +++ b/src/scaffoldmaker/utils/vector.py @@ -44,28 +44,28 @@ def addVectors(v1,v2,s1=1.0,s2=1.0): return [(s1 * v1[c] + s2 * v2[c]) for c in range(len(v1))] -def scalarProjectionOfV1OnV2(v1, v2): +def scalarProjection(v1, v2): """ :return: Scalar projection of v1 onto v2. """ return dotproduct(v1, normalise(v2)) -def vectorProjectionOfV1OnV2(v1, v2): +def vectorProjection(v1, v2): """ Calculate vector projection of v1 on v2 :return: A projection vector. """ - s1 = scalarProjectionOfV1OnV2(v1, v2) + s1 = scalarProjection(v1, v2) return scalarProduct(s1, normalise(v2)) -def vectorRejectionOfV1OnV2(v1, v2): +def vectorRejection(v1, v2): """ Calculate vector rejection of v1 on v2 :return: A rejection vector. """ - v1p = vectorProjectionOfV1OnV2(v1, v2) + v1p = vectorProjection(v1, v2) return addVectors(v1, v1p, 1.0, -1.0) @@ -79,13 +79,12 @@ def scalarProduct(s, v): return [s * v[c] for c in range(len(v))] -def isVectorsParallel(v1, v2): +def parallelVectors(v1, v2): """ :return: True if the vectors are parallel. """ assert (len(v2) == len(v1)), 'Vectors lengths are not the same.' - s = magnitude(v2)/magnitude(v1) - for c in range(len(v2)): - if v2[c] != s * v1[c]: - return False - return True + TOL = 1.0e-6 + if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2))/2.0: + return True + return False From db4e15a77e5874e21012636b0cc7f8df78f2d88d Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 24 Sep 2020 11:30:28 +1200 Subject: [PATCH 08/10] Minor change --- .../meshtypes/meshtype_1d_path1.py | 54 ++++++++++--------- src/scaffoldmaker/utils/vector.py | 4 +- tests/test_colon.py | 3 +- tests/test_smallintestine.py | 3 +- 4 files changed, 36 insertions(+), 28 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index ea25b422..87a09d56 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -132,7 +132,7 @@ def generateBaseMesh(cls, region, options): @classmethod def smoothPath(cls, region, options, mode : DerivativeScalingMode): - x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1]) + x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1])[0:2] d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode) setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ]) @@ -140,7 +140,7 @@ def smoothPath(cls, region, options, mode : DerivativeScalingMode): def makeD2Normal(cls, region, options): if not options['D2 derivatives']: return - d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2]) + d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2])[0:2] for c in range(len(d1)): td2 = vector.vectorRejection(d2[c], d1[c]) d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c])) @@ -152,12 +152,12 @@ def makeD3Normal(cls, region, options): return if options['D2 derivatives']: d1, d2, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, - Node.VALUE_LABEL_D_DS3]) + Node.VALUE_LABEL_D_DS3])[0:3] for c in range(len(d1)): d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) else: - d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3]) + d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3])[0:2] for c in range(len(d1)): td3 = vector.vectorRejection(d3[c], d1[c]) d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c])) @@ -196,9 +196,10 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None): cx = [] cd1 = [] cd2 = [] - cd12 = [] cd3 = [] + cd12 = [] cd13 = [] + cd23 = [] nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) if groupName: group = fieldmodule.findFieldByName(groupName).castGroup() @@ -211,25 +212,28 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None): node = nodeIter.next() while node.isValid(): cache.setNode(node) - result, x = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, componentsCount) - result, d1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, componentsCount) - result, d2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, componentsCount) - result, d12 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, componentsCount) - result, d3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, componentsCount) - result, d13 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) - for c in range(componentsCount, 3): - x.append(0.0) - d1.append(0.0) - d2.append(0.0) - d12.append(0.0) - d3.append(0.0) - d13.append(0.0) - cx.append(x) - cd1.append(d1) - cd2.append(d2) - cd12.append(d12) - cd3.append(d3) - cd13.append(d13) + for label in valueLabels: + if label == Node.VALUE_LABEL_VALUE: + result, x = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, componentsCount) + cx.append(x) + elif label == Node.VALUE_LABEL_D_DS1: + result, d1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, componentsCount) + cd1.append(d1) + elif label == Node.VALUE_LABEL_D_DS2: + result, d2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, componentsCount) + cd2.append(d2) + elif label == Node.VALUE_LABEL_D_DS3: + result, d3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, componentsCount) + cd3.append(d3) + elif label == Node.VALUE_LABEL_D2_DS1DS2: + result, d12 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, componentsCount) + cd12.append(d12) + elif label == Node.VALUE_LABEL_D2_DS1DS3: + result, d13 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) + cd13.append(d13) + elif label == Node.VALUE_LABEL_D2_DS2DS3: + result, d23 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) + cd23.append(d23) node = nodeIter.next() result = () @@ -246,6 +250,8 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None): result += (cd12,) if label == Node.VALUE_LABEL_D2_DS1DS3: result += (cd13,) + if label == Node.VALUE_LABEL_D2_DS2DS3: + result += (cd23,) return result diff --git a/src/scaffoldmaker/utils/vector.py b/src/scaffoldmaker/utils/vector.py index da4f27f9..340753a1 100644 --- a/src/scaffoldmaker/utils/vector.py +++ b/src/scaffoldmaker/utils/vector.py @@ -84,7 +84,7 @@ def parallelVectors(v1, v2): :return: True if the vectors are parallel. """ assert (len(v2) == len(v1)), 'Vectors lengths are not the same.' - TOL = 1.0e-6 - if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2))/2.0: + TOL = 1.0e-6/2.0 + if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2)): return True return False diff --git a/tests/test_colon.py b/tests/test_colon.py index 991b059e..b1afded5 100644 --- a/tests/test_colon.py +++ b/tests/test_colon.py @@ -26,6 +26,7 @@ def test_colon1(self): 'Test line': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, + 'D2 derivatives':True, 'Length': 1.0, 'Number of elements': 1 }, @@ -83,7 +84,7 @@ def test_colon1(self): tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx = extractPathParametersFromRegion(tmpRegion)[0] + cx = extractPathParametersFromRegion(tmpRegion, [Node.VALUE_LABEL_VALUE])[0] self.assertEqual(2, len(cx)) assertAlmostEqualList(self, cx[0], [ 163.7, -25.2, 12.2 ], 1.0E-6) assertAlmostEqualList(self, cx[1], [ 117.2, 32.8, -2.6 ], 1.0E-6) diff --git a/tests/test_smallintestine.py b/tests/test_smallintestine.py index cb4297b4..7cdbc96c 100644 --- a/tests/test_smallintestine.py +++ b/tests/test_smallintestine.py @@ -24,6 +24,7 @@ def test_smallintestine1(self): centralPathDefaultScaffoldPackages = { 'Test line': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { + 'D2 derivatives': True, 'Coordinate dimensions': 3, 'Length': 1.0, 'Number of elements': 3 @@ -62,7 +63,7 @@ def test_smallintestine1(self): tmpRegion = region.createRegion() centralPath.generate(tmpRegion) - cx = extractPathParametersFromRegion(tmpRegion)[0] + cx = extractPathParametersFromRegion(tmpRegion, [Node.VALUE_LABEL_VALUE])[0] self.assertEqual(4, len(cx)) assertAlmostEqualList(self, cx[0], [ -2.3, 18.5, -4.4 ], 1.0E-6) assertAlmostEqualList(self, cx[2], [ -18.3, 12.6, -1.5 ], 1.0E-6) From e92a372acf51ee8e80c14a8d5cc7825f940982e0 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Thu, 24 Sep 2020 15:35:42 +1200 Subject: [PATCH 09/10] Minor changes. --- .../meshtypes/meshtype_1d_path1.py | 55 +++---------------- 1 file changed, 8 insertions(+), 47 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 87a09d56..772e571a 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -193,13 +193,8 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None): assert componentsCount in [ 1, 2, 3 ], 'extractPathParametersFromRegion. Invalid coordinates number of components' cache = fieldmodule.createFieldcache() - cx = [] - cd1 = [] - cd2 = [] - cd3 = [] - cd12 = [] - cd13 = [] - cd23 = [] + valueLabelsCount = len(valueLabels) + returnValues = [[] for i in range(valueLabelsCount)] nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) if groupName: group = fieldmodule.findFieldByName(groupName).castGroup() @@ -212,48 +207,14 @@ def extractPathParametersFromRegion(region, valueLabels, groupName=None): node = nodeIter.next() while node.isValid(): cache.setNode(node) - for label in valueLabels: - if label == Node.VALUE_LABEL_VALUE: - result, x = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, componentsCount) - cx.append(x) - elif label == Node.VALUE_LABEL_D_DS1: - result, d1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, componentsCount) - cd1.append(d1) - elif label == Node.VALUE_LABEL_D_DS2: - result, d2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, componentsCount) - cd2.append(d2) - elif label == Node.VALUE_LABEL_D_DS3: - result, d3 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, componentsCount) - cd3.append(d3) - elif label == Node.VALUE_LABEL_D2_DS1DS2: - result, d12 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, componentsCount) - cd12.append(d12) - elif label == Node.VALUE_LABEL_D2_DS1DS3: - result, d13 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) - cd13.append(d13) - elif label == Node.VALUE_LABEL_D2_DS2DS3: - result, d23 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, componentsCount) - cd23.append(d23) + for i in range(valueLabelsCount): + result, values = coordinates.getNodeParameters(cache, -1, valueLabels[i], 1, componentsCount) + for c in range(componentsCount, 3): + values.append(0.0) + returnValues[i].append(values) node = nodeIter.next() - result = () - for label in valueLabels: - if label == Node.VALUE_LABEL_VALUE: - result += (cx,) - if label == Node.VALUE_LABEL_D_DS1: - result += (cd1,) - if label == Node.VALUE_LABEL_D_DS2: - result += (cd2,) - if label == Node.VALUE_LABEL_D_DS3: - result += (cd3,) - if label == Node.VALUE_LABEL_D2_DS1DS2: - result += (cd12,) - if label == Node.VALUE_LABEL_D2_DS1DS3: - result += (cd13,) - if label == Node.VALUE_LABEL_D2_DS2DS3: - result += (cd23,) - - return result + return returnValues def setPathParameters(region, nodeValueLabels, nodeValues): From 19436447ae122a74ee29339e6d73778ff38d70f5 Mon Sep 17 00:00:00 2001 From: Elias Ghadam Soltani Date: Fri, 25 Sep 2020 10:05:52 +1200 Subject: [PATCH 10/10] Minor changes --- src/scaffoldmaker/meshtypes/meshtype_1d_path1.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py index 772e571a..35618f55 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_path1.py @@ -132,7 +132,7 @@ def generateBaseMesh(cls, region, options): @classmethod def smoothPath(cls, region, options, mode : DerivativeScalingMode): - x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1])[0:2] + x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1]) d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode) setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ]) @@ -140,7 +140,7 @@ def smoothPath(cls, region, options, mode : DerivativeScalingMode): def makeD2Normal(cls, region, options): if not options['D2 derivatives']: return - d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2])[0:2] + d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2]) for c in range(len(d1)): td2 = vector.vectorRejection(d2[c], d1[c]) d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c])) @@ -152,12 +152,12 @@ def makeD3Normal(cls, region, options): return if options['D2 derivatives']: d1, d2, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, - Node.VALUE_LABEL_D_DS3])[0:3] + Node.VALUE_LABEL_D_DS3]) for c in range(len(d1)): d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c])) setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3]) else: - d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3])[0:2] + d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3]) for c in range(len(d1)): td3 = vector.vectorRejection(d3[c], d1[c]) d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c]))