Skip to content

Commit

Permalink
Fix multiple transition node smoothing, add test
Browse files Browse the repository at this point in the history
  • Loading branch information
@rchristie
rchristie committed Oct 10, 2024
1 parent 485839a commit 09c3e60
Show file tree
Hide file tree
Showing 2 changed files with 131 additions and 68 deletions.
136 changes: 68 additions & 68 deletions src/scaffoldmaker/utils/tubenetworkmesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -721,81 +721,81 @@ def _generateCoreCoordinates(self, n2, centre):
cbd1.append(row_d1)
cbd3.append(row_d3)

for i in range(self._elementsCountTransition - 1):
for lst in (ctx, ctd1, ctd3):
lst.append([None] * self._elementsCountAround)
ix = self._rimCoordinates[0][n2][0]
id1 = self._rimCoordinates[1][n2][0]
id3 = self._rimCoordinates[3][n2][0]
start_bn3 = minorBoxSize // 2
topLeft_n1 = minorBoxSize - start_bn3
topRight_n1 = topLeft_n1 + majorBoxSize
bottomRight_n1 = topRight_n1 + minorBoxSize
bottomLeft_n1 = bottomRight_n1 + majorBoxSize
for n1 in range(self._elementsCountAround):
if n1 <= topLeft_n1:
bn1 = 0
bn3 = start_bn3 + n1
if n1 < topLeft_n1:
if self._elementsCountTransition > 1:
for i in range(self._elementsCountTransition - 1):
for lst in (ctx, ctd1, ctd3):
lst.append([None] * self._elementsCountAround)
ix = self._rimCoordinates[0][n2][0]
id1 = self._rimCoordinates[1][n2][0]
id3 = self._rimCoordinates[3][n2][0]
start_bn3 = minorBoxSize // 2
topLeft_n1 = minorBoxSize - start_bn3
topRight_n1 = topLeft_n1 + majorBoxSize
bottomRight_n1 = topRight_n1 + minorBoxSize
bottomLeft_n1 = bottomRight_n1 + majorBoxSize
for n1 in range(self._elementsCountAround):
if n1 <= topLeft_n1:
bn1 = 0
bn3 = start_bn3 + n1
if n1 < topLeft_n1:
start_d1 = cbd3[bn1][bn3]
start_d3 = [-d for d in cbd1[bn1][bn3]]
else:
start_d1 = add(cbd3[bn1][bn3], cbd1[bn1][bn3])
start_d3 = sub(cbd3[bn1][bn3], cbd1[bn1][bn3])
elif n1 <= topRight_n1:
bn1 = n1 - topLeft_n1
bn3 = minorBoxSize
if n1 < topRight_n1:
start_d1 = cbd1[bn1][bn3]
start_d3 = cbd3[bn1][bn3]
else:
start_d1 = sub(cbd1[bn1][bn3], cbd3[bn1][bn3])
start_d3 = add(cbd1[bn1][bn3], cbd3[bn1][bn3])
elif n1 <= bottomRight_n1:
bn1 = majorBoxSize
bn3 = minorBoxSize - (n1 - topRight_n1)
if n1 < bottomRight_n1:
start_d1 = [-d for d in cbd3[bn1][bn3]]
start_d3 = cbd1[bn1][bn3]
else:
start_d1 = [-d for d in add(cbd1[bn1][bn3], cbd3[bn1][bn3])]
start_d3 = sub(cbd1[bn1][bn3], cbd3[bn1][bn3])
elif n1 <= bottomLeft_n1:
bn1 = majorBoxSize - (n1 - bottomRight_n1)
bn3 = 0
if n1 < bottomLeft_n1:
start_d1 = [-d for d in cbd1[bn1][bn3]]
start_d3 = [-d for d in cbd3[bn1][bn3]]
else:
start_d1 = sub(cbd3[bn1][bn3], cbd1[bn1][bn3])
start_d3 = [-d for d in add(cbd1[bn1][bn3], cbd3[bn1][bn3])]
else:
bn1 = 0
bn3 = n1 - bottomLeft_n1
start_d1 = cbd3[bn1][bn3]
start_d3 = [-d for d in cbd1[bn1][bn3]]
else:
start_d1 = add(cbd3[bn1][bn3], cbd1[bn1][bn3])
start_d3 = sub(cbd3[bn1][bn3], cbd1[bn1][bn3])
elif n1 <= topRight_n1:
bn1 = n1 - topLeft_n1
bn3 = minorBoxSize
if n1 < topRight_n1:
start_d1 = cbd1[bn1][bn3]
start_d3 = cbd3[bn1][bn3]
else:
start_d1 = sub(cbd1[bn1][bn3], cbd3[bn1][bn3])
start_d3 = add(cbd1[bn1][bn3], cbd3[bn1][bn3])
elif n1 <= bottomRight_n1:
bn1 = majorBoxSize
bn3 = minorBoxSize - (n1 - topRight_n1)
if n1 < bottomRight_n1:
start_d1 = [-d for d in cbd3[bn1][bn3]]
start_d3 = cbd1[bn1][bn3]
else:
start_d1 = [-d for d in add(cbd1[bn1][bn3], cbd3[bn1][bn3])]
start_d3 = sub(cbd1[bn1][bn3], cbd3[bn1][bn3])
elif n1 <= bottomLeft_n1:
bn1 = majorBoxSize - (n1 - bottomRight_n1)
bn3 = 0
if n1 < bottomLeft_n1:
start_d1 = [-d for d in cbd1[bn1][bn3]]
start_d3 = [-d for d in cbd3[bn1][bn3]]
else:
start_d1 = sub(cbd3[bn1][bn3], cbd1[bn1][bn3])
start_d3 = [-d for d in add(cbd1[bn1][bn3], cbd3[bn1][bn3])]
else:
bn1 = 0
bn3 = n1 - bottomLeft_n1
start_d1 = cbd3[bn1][bn3]
start_d3 = [-d for d in cbd1[bn1][bn3]]
start_x = cbx[bn1][bn3]

nx = [start_x, ix[n1]]
nd3before = [[self._elementsCountTransition * d for d in start_d3], id3[n1]]
nd3 = [nd3before[0], computeCubicHermiteEndDerivative(nx[0], nd3before[0], nx[1], nd3before[1])]
tx, td3, pe, pxi, psf = sampleCubicHermiteCurvesSmooth(
nx, nd3, self._elementsCountTransition,
derivativeMagnitudeStart=magnitude(nd3[0]) / self._elementsCountTransition,
derivativeMagnitudeEnd=magnitude(nd3[1]) / self._elementsCountTransition)
delta_id1 = sub(id1[n1], start_d1)
td1 = interpolateSampleCubicHermite([start_d1, id1[n1]], [delta_id1, delta_id1], pe, pxi, psf)[0]

for n3 in range(1, self._elementsCountTransition):
ctx[n3 - 1][n1] = tx[n3]
ctd1[n3 - 1][n1] = td1[n3]
ctd3[n3 - 1][n1] = td3[n3]
start_x = cbx[bn1][bn3]

nx = [start_x, ix[n1]]
nd3before = [[self._elementsCountTransition * d for d in start_d3], id3[n1]]
nd3 = [nd3before[0], computeCubicHermiteEndDerivative(nx[0], nd3before[0], nx[1], nd3before[1])]
tx, td3, pe, pxi, psf = sampleCubicHermiteCurvesSmooth(
nx, nd3, self._elementsCountTransition,
derivativeMagnitudeStart=magnitude(nd3[0]) / self._elementsCountTransition,
derivativeMagnitudeEnd=magnitude(nd3[1]) / self._elementsCountTransition)
delta_id1 = sub(id1[n1], start_d1)
td1 = interpolateSampleCubicHermite([start_d1, id1[n1]], [delta_id1, delta_id1], pe, pxi, psf)[0]

for n3 in range(1, self._elementsCountTransition):
ctx[n3 - 1][n1] = tx[n3]
ctd1[n3 - 1][n1] = td1[n3]
ctd3[n3 - 1][n1] = td3[n3]

# smooth td1 around:
for n3 in range(1, self._elementsCountTransition):
ctd1[n3 - 1] = smoothCubicHermiteDerivativesLoop(ctx[n3 - 1], ctd1[n3 - 1], fixAllDirections=False)


return cbx, cbd1, cbd3, ctx, ctd1, ctd3

def _determineCoreD2Derivatives(self, boxx, boxd1, boxd3, transx, transd1, transd3):
Expand Down
63 changes: 63 additions & 0 deletions tests/test_network.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -473,6 +473,69 @@ def test_3d_tube_network_bifurcation_core(self):
self.assertAlmostEqual(volume, 0.09907643906540035, delta=X_TOL)
self.assertAlmostEqual(surfaceArea, 2.0238210110948307, delta=X_TOL)

def test_3d_tube_network_line_core_transition2(self):
"""
Test line 3-D tube network with solid core and 2 transition elements.
"""
scaffoldPackage = ScaffoldPackage(MeshType_3d_tubenetwork1, defaultParameterSetName="Default")
settings = scaffoldPackage.getScaffoldSettings()
networkLayoutScaffoldPackage = settings["Network layout"]
networkLayoutSettings = networkLayoutScaffoldPackage.getScaffoldSettings()
self.assertTrue(networkLayoutSettings["Define inner coordinates"])
self.assertEqual(13, len(settings))
self.assertEqual(8, settings["Number of elements around"])
self.assertEqual(1, settings["Number of elements through shell"])
self.assertEqual([0], settings["Annotation numbers of elements around"])
self.assertEqual(4.0, settings["Target element density along longest segment"])
self.assertEqual([0], settings["Annotation numbers of elements along"])
self.assertFalse(settings["Use linear through shell"])
self.assertFalse(settings["Show trim surfaces"])
self.assertFalse(settings["Core"])
self.assertEqual(2, settings["Number of elements across core box minor"])
self.assertEqual(1, settings["Number of elements across core transition"])
self.assertEqual([0], settings["Annotation numbers of elements across core box minor"])
settings["Core"] = True
settings["Number of elements across core transition"] = 2

context = Context("Test")
region = context.getDefaultRegion()
self.assertTrue(region.isValid())
scaffoldPackage.generate(region)

fieldmodule = region.getFieldmodule()
mesh3d = fieldmodule.findMeshByDimension(3)

self.assertEqual(112, mesh3d.getSize())
nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
self.assertEqual(165, nodes.getSize())
coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
self.assertTrue(coordinates.isValid())

X_TOL = 1.0E-6

minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
assertAlmostEqualList(self, minimums, [0.0, -0.1, -0.1], X_TOL)
assertAlmostEqualList(self, maximums, [1.0, 0.1, 0.1], X_TOL)

with ChangeManager(fieldmodule):
one = fieldmodule.createFieldConstant(1.0)
isExterior = fieldmodule.createFieldIsExterior()
mesh2d = fieldmodule.findMeshByDimension(2)
fieldcache = fieldmodule.createFieldcache()

volumeField = fieldmodule.createFieldMeshIntegral(one, coordinates, mesh3d)
volumeField.setNumbersOfPoints(4)
result, volume = volumeField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)

surfaceAreaField = fieldmodule.createFieldMeshIntegral(isExterior, coordinates, mesh2d)
surfaceAreaField.setNumbersOfPoints(4)
result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)

self.assertAlmostEqual(volume, 0.0313832204833548, delta=X_TOL)
self.assertAlmostEqual(surfaceArea, 0.6907602069977625, delta=X_TOL)

def test_3d_tube_network_sphere_cube(self):
"""
Test sphere cube 3-D tube network is generated correctly.
Expand Down

0 comments on commit 09c3e60

Please sign in to comment.