Skip to content

Add dome ends to tube network mesh #293

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
rchristie wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Add dome ends to tube network mesh #293

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 12 additions & 48 deletions src/scaffoldmaker/meshtypes/meshtype_3d_uterus1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,10 @@
from scaffoldmaker.utils.interpolation import smoothCurveSideCrossDerivatives, smoothCubicHermiteDerivativesLine, \
interpolateCubicHermite, sampleCubicHermiteCurves, computeCubicHermiteDerivativeScaling, \
interpolateLagrangeHermiteDerivative
from scaffoldmaker.utils.networkmesh import NetworkMesh, pathValueLabels
from scaffoldmaker.utils.tubenetworkmesh import TubeNetworkMeshBuilder, TubeNetworkMeshGenerateData, \
PatchTubeNetworkMeshSegment
from scaffoldmaker.utils.zinc_utils import group_add_connected_elements, get_nodeset_path_ordered_field_parameters, \
find_first_node_conditional, get_node_mesh_location
from scaffoldmaker.utils.networkmesh import NetworkMesh
from scaffoldmaker.utils.tubenetworkmesh import BodyTubeNetworkMeshBuilder, TubeNetworkMeshGenerateData
from scaffoldmaker.utils.zinc_utils import (
find_first_node_conditional, get_node_mesh_location, group_add_connected_elements)
import copy
import math

Expand All @@ -43,7 +42,7 @@ def __init__(self, region, meshDimension, coordinateFieldName="coordinates",
self._fundusGroup = self.getOrCreateAnnotationGroup(get_uterus_term("fundus of uterus"))
self._bodyGroup = self.getOrCreateAnnotationGroup(get_uterus_term("body of uterus"))
self._bodyNotCervixGroup = self.getOrCreateAnnotationGroup(("body not cervix", ""))
# force these annotation group names in base class
# force these annotation group names in base class with these names
self._leftGroup = self.getOrCreateAnnotationGroup(("left uterus", ""))
self._rightGroup = self.getOrCreateAnnotationGroup(("right uterus", ""))
self._dorsalGroup = self.getOrCreateAnnotationGroup(("dorsal uterus", ""))
Expand All @@ -65,48 +64,23 @@ def getBodyNotCervixMeshGroup(self):
return self._bodyNotCervixGroup.getMeshGroup(self._mesh)


class UterusTubeNetworkMeshBuilder(TubeNetworkMeshBuilder):
class UterusTubeNetworkMeshBuilder(BodyTubeNetworkMeshBuilder):
"""
Adds left, right, dorsal, ventral, fundus, body annotations.
Future: derive from BodyTubeNetworkMeshBuilder to get left/right/dorsal/ventral.
Adds body, cervix, fundus, oviduct annotations.
Derived from BodyTubeNetworkMeshBuilder to get left/right/dorsal/ventral.
"""
def createSegment(self, networkSegment):
if networkSegment.isPatch():
pathParametersList = [get_nodeset_path_ordered_field_parameters(
self._layoutNodes, self._layoutCoordinates, pathValueLabels,
networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions())]
if self._layoutInnerCoordinates:
pathParametersList.append(get_nodeset_path_ordered_field_parameters(
self._layoutNodes, self._layoutInnerCoordinates, pathValueLabels,
networkSegment.getNodeIdentifiers(), networkSegment.getNodeVersions()))
elementsCountAround = self._defaultElementsCountAround
elementsCountCoreBoxMinor = self._defaultElementsCountCoreBoxMinor
coreBoundaryScalingMode = self._defaultCoreBoundaryScalingMode

return PatchTubeNetworkMeshSegment(networkSegment, pathParametersList, elementsCountAround,
self._elementsCountThroughShell, self._isCore, elementsCountCoreBoxMinor,
self._elementsCountTransition, coreBoundaryScalingMode)

return super(UterusTubeNetworkMeshBuilder, self).createSegment(networkSegment)

def generateMesh(self, generateData):
super(UterusTubeNetworkMeshBuilder, self).generateMesh(generateData)
# build temporary left/right dorsal/ventral groups
super(UterusTubeNetworkMeshBuilder, self).generateMesh(generateData, dorsalD3Side=False)
mesh = generateData.getMesh()
leftOviductMeshGroup = generateData.getLeftOviductMeshGroup()
rightOviductMeshGroup = generateData.getRightOviductMeshGroup()
fundusMeshGroup = generateData.getFundusMeshGroup()
bodyMeshGroup = generateData.getBodyMeshGroup()
bodyNotCervixMeshGroup = generateData.getBodyNotCervixMeshGroup()
leftMeshGroup = generateData.getLeftMeshGroup()
rightMeshGroup = generateData.getRightMeshGroup()
dorsalMeshGroup = generateData.getDorsalMeshGroup()
ventralMeshGroup = generateData.getVentralMeshGroup()
for networkSegment in self._networkMesh.getNetworkSegments():
segment = self._segments[networkSegment]
annotationTerms = segment.getAnnotationTerms()
segment.addSideD3ElementsToMeshGroup(True, ventralMeshGroup)
segment.addSideD3ElementsToMeshGroup(False, dorsalMeshGroup)
for annotationTerm in annotationTerms:
if "oviduct" in annotationTerm[0]:
elementsCountRim = segment.getElementsCountRim()
Expand Down Expand Up @@ -140,16 +114,6 @@ def generateMesh(self, generateData):
segment.addAllElementsToMeshGroup(bodyNotCervixMeshGroup)
if "cervix" in annotationTerm[0]:
segment.addAllElementsToMeshGroup(bodyMeshGroup)
if "left" in annotationTerm[0]:
segment.addAllElementsToMeshGroup(leftMeshGroup)
break
if "right" in annotationTerm[0]:
segment.addAllElementsToMeshGroup(rightMeshGroup)
break
else:
# segment on main axis
segment.addSideD2ElementsToMeshGroup(True, leftMeshGroup)
segment.addSideD2ElementsToMeshGroup(False, rightMeshGroup)


class Septum:
Expand Down Expand Up @@ -893,7 +857,7 @@ def getDefaultOptions(cls, parameterSetName="Default"):
options["Structure"] = (
"1-2-3-4-5-6-7-8-9-10-11-12-13-14-15-16-36.1,"
"17-18-19-20-21-22-23-24-25-26-27-28-29-30-31-32-36.2,"
"#-33-34-35-36.3,"
"#33-34-35-36.3,"
"36.4-37-38-39,"
"39-40,"
"40-41-42")
Expand All @@ -918,7 +882,7 @@ def getDefaultOptions(cls, parameterSetName="Default"):
options["Structure"] = (
"1-2-3-4-5-6-7-8-9-31.1,"
"10-11-12-13-14-15-16-17-18-31.2,"
"#-19-20-21-22-23-24-25-26-27-28-29-30-31.3,"
"#19-20-21-22-23-24-25-26-27-28-29-30-31.3,"
"31.4-32-33-34-35-36-37-38-39-40-41-42-43,"
"43-44,"
"44-45-46-47-48")
Expand All @@ -943,7 +907,7 @@ def getDefaultOptions(cls, parameterSetName="Default"):
options["Structure"] = (
"1-2-3-4-5-6-7-8-23.1,"
"9-10-11-12-13-14-15-16-23.2,"
"#-17-18-19-20-21-22-23.3,"
"#17-18-19-20-21-22-23.3,"
"23.4-24-25-26-27-28-29,"
"29-30-31,"
"31-32-33-34-35-36-37-38")
Expand Down
12 changes: 9 additions & 3 deletions src/scaffoldmaker/utils/eft_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1065,8 +1065,12 @@ def determineCubicHermiteSerendipityEft(mesh, nodeParameters, nodeLayouts):
on = n ^ (1 << ed)
if nodeOrder.index(on) > nodeOrder.index(n):
derivativeMagnitude = magnitude(elementDerivative)
newMagnitude = 2.0 / ((1.0 / oldDeltaMagnitude) + (1.0 / derivativeMagnitude))
elementDerivative = mult(elementDerivative, newMagnitude / derivativeMagnitude)
if (oldDeltaMagnitude > 0.0) and (derivativeMagnitude > 0.0):
newMagnitude = 2.0 / ((1.0 / oldDeltaMagnitude) + (1.0 / derivativeMagnitude))
else:
newMagnitude = (oldDeltaMagnitude + derivativeMagnitude) / 2.0
if derivativeMagnitude > 0.0:
elementDerivative = mult(elementDerivative, newMagnitude / derivativeMagnitude)
# update other node delta to work in with this derivative
otherElementDerivative = (
interpolateHermiteLagrangeDerivative(nodeParameters[n][0], elementDerivative, nodeParameters[on][0], 1.0)
Expand Down Expand Up @@ -1144,7 +1148,9 @@ def addTricubicHermiteSerendipityEftParameterScaling(eft, scalefactors, nodePara
xb = nodeParameters[nb][0]
db = nodeParameters[nb][3]
dbScaled = computeCubicHermiteEndDerivative(xa, da, xb, db)
scalefactor = magnitude(dbScaled) / magnitude(db)
mag_db = magnitude(db)
mag_dbScaled = magnitude(dbScaled)
scalefactor = mag_dbScaled / mag_db if ((mag_db > 0.0) and (mag_dbScaled > 0.0)) else 0.0
addScalefactors.append(scalefactor)
if version == 1:
newScalefactors = scalefactors + addScalefactors if scalefactors else addScalefactors
Expand Down
23 changes: 17 additions & 6 deletions src/scaffoldmaker/utils/interpolation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,11 +187,14 @@ def computeCubicHermiteEndDerivative(v1, d1, v2, d2_in):
:return: Scaled d2
"""
d1_mag = magnitude(d1)
d2 = set_magnitude(d2_in, 0.5 * d1_mag)
d2_in_mag = magnitude(d2_in)
if d2_in_mag == 0.0:
return copy.copy(d2_in)
d2 = mult(d2_in, 0.5 * d1_mag / d2_in_mag)
for iters in range(100):
arcLength = getCubicHermiteArcLength(v1, d1, v2, d2)
d2_mag = 2.0 * arcLength - d1_mag
d2 = set_magnitude(d2_in, d2_mag)
d2 = mult(d2_in, d2_mag / d2_in_mag)
if math.fabs(2.0 * arcLength - d1_mag - d2_mag) < (1.0E-6 * arcLength):
break
else:
Expand Down Expand Up @@ -267,6 +270,10 @@ def getCubicHermiteTrimmedCurvesLengths(cx, cd1, startLocation=None, endLocation
:return: Length before start, length between start and end, length after end, list of lengths to nodes.
"""
elementsCount = len(cx) - 1
if startLocation:
assert 0 <= startLocation[0] < elementsCount
if endLocation:
assert 0 <= endLocation[0] < elementsCount
lengthToNode = [0.0]
length = 0.0
for e in range(elementsCount):
Expand Down Expand Up @@ -849,7 +856,8 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
if not fixStartDerivative:
if fixAllDirections or fixStartDirection:
mag = 2.0*arcLengths[0] - magnitude(lastmd1[1])
md1[0] = set_magnitude(nd1[0], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ]
old_mag = magnitude(nd1[0])
md1[0] = mult(nd1[0], mag / old_mag) if ((old_mag > 0.0) and (mag > 0.0)) else [0.0] * componentsCount
else:
md1[0] = interpolateLagrangeHermiteDerivative(nx[0], nx[1], lastmd1[1], 0.0)
# middle
Expand All @@ -876,12 +884,14 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
mag = 0.5 * (dnm + dn)
else: # harmonicMeanMagnitude
mag = 2.0 / (1.0 / dnm + 1.0 / dn)
md1[n] = set_magnitude(md1[n], mag)
old_mag = magnitude(md1[n])
md1[n] = mult(md1[n], mag / old_mag) if (old_mag > 0.0) else [0.0] * componentsCount
# end
if not fixEndDerivative:
if fixAllDirections or fixEndDirection:
mag = 2.0*arcLengths[-1] - magnitude(lastmd1[-2])
md1[-1] = set_magnitude(nd1[-1], mag) if (mag > 0.0) else [ 0.0, 0.0, 0.0 ]
old_mag = magnitude(nd1[-1])
md1[-1] = mult(nd1[-1], mag / old_mag) if ((old_mag > 0.0) and (mag > 0.0)) else [0.0] * componentsCount
else:
md1[-1] = interpolateHermiteLagrangeDerivative(nx[-2], lastmd1[-2], nx[-1], 1.0)
if instrument:
Expand Down Expand Up @@ -1105,7 +1115,8 @@ def smoothCubicHermiteDerivativesLoop(nx, nd1,
mag = 0.5*(arcLengths[nm] + arcLengths[n])
else: # harmonicMeanMagnitude
mag = 2.0/(1.0/arcLengths[nm] + 1.0/arcLengths[n])
md1[n] = set_magnitude(md1[n], mag)
old_mag = magnitude(md1[n])
md1[n] = mult(md1[n], mag / old_mag) if (old_mag > 0.0) else [0.0] * componentsCount
if instrument:
print('iter', iter + 1, md1)
dtol = tol*sum(arcLengths)/len(arcLengths)
Expand Down
Loading
Loading