Skip to content

Commit d80af4e

Browse files
authored
Merge pull request ABI-Software#279 from rchristie/ellipsoid
Add ellipsoid scaffold
2 parents a022a33 + e383e9b commit d80af4e

File tree

8 files changed

+3197
-9
lines changed

8 files changed

+3197
-9
lines changed
Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
"""
2+
Generates a solid ellipsoid of hexahedral elements.
3+
"""
4+
import math
5+
from cmlibs.utils.zinc.field import find_or_create_field_coordinates
6+
from scaffoldmaker.annotation.annotationgroup import AnnotationGroup
7+
from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
8+
from scaffoldmaker.utils.ellipsoidmesh import EllipsoidMesh, EllipsoidSurfaceD3Mode
9+
from scaffoldmaker.utils.meshrefinement import MeshRefinement
10+
11+
12+
class MeshType_3d_ellipsoid1(Scaffold_base):
13+
"""
14+
Generates a solid ellipsoid of hexahedral elements.
15+
"""
16+
17+
@classmethod
18+
def getName(cls):
19+
return "3D Ellipsoid 1"
20+
21+
@classmethod
22+
def getDefaultOptions(cls, parameterSetName='Default'):
23+
options = {
24+
"Number of elements across axis 1": 4,
25+
"Number of elements across axis 2": 6,
26+
"Number of elements across axis 3": 8,
27+
"2D surface only": False,
28+
"Number of transition elements": 1,
29+
"Axis length x": 1.0,
30+
"Axis length y": 1.5,
31+
"Axis length z": 2.0,
32+
"Axis 2 x-rotation degrees": 0.0,
33+
"Axis 3 x-rotation degrees": 90.0,
34+
"Advanced n-way derivative factor": 0.6,
35+
"Advanced surface D3 mode": EllipsoidSurfaceD3Mode.SURFACE_NORMAL.value,
36+
"Refine": False,
37+
"Refine number of elements": 4,
38+
}
39+
return options
40+
41+
@classmethod
42+
def getOrderedOptionNames(cls):
43+
return [
44+
"Number of elements across axis 1",
45+
"Number of elements across axis 2",
46+
"Number of elements across axis 3",
47+
"Number of transition elements",
48+
"2D surface only",
49+
"Axis length x",
50+
"Axis length y",
51+
"Axis length z",
52+
"Axis 2 x-rotation degrees",
53+
"Axis 3 x-rotation degrees",
54+
"Advanced n-way derivative factor",
55+
"Advanced surface D3 mode",
56+
"Refine",
57+
"Refine number of elements"
58+
]
59+
60+
@classmethod
61+
def checkOptions(cls, options):
62+
dependent_changes = False
63+
max_transition_count = None
64+
for key in [
65+
"Number of elements across axis 1",
66+
"Number of elements across axis 2",
67+
"Number of elements across axis 3"
68+
]:
69+
if options[key] < 4:
70+
options[key] = 4
71+
elif options[key] % 2:
72+
options[key] += 1
73+
transition_count = (options[key] // 2) - 1
74+
if (max_transition_count is None) or (transition_count < max_transition_count):
75+
max_transition_count = transition_count
76+
77+
if options["Number of transition elements"] < 1:
78+
options["Number of transition elements"] = 1
79+
elif options["Number of transition elements"] > max_transition_count:
80+
options["Number of transition elements"] = max_transition_count
81+
dependent_changes = True
82+
83+
for key in [
84+
"Axis length x",
85+
"Axis length y",
86+
"Axis length z"
87+
]:
88+
if options[key] <= 0.0:
89+
options[key] = 1.0
90+
91+
if options["Advanced n-way derivative factor"] < 0.1:
92+
options["Advanced n-way derivative factor"] = 0.1
93+
elif options["Advanced n-way derivative factor"] > 1.0:
94+
options["Advanced n-way derivative factor"] = 1.0
95+
96+
try:
97+
mode = EllipsoidSurfaceD3Mode(options["Advanced surface D3 mode"])
98+
except ValueError:
99+
options["Advanced surface D3 mode"] = EllipsoidSurfaceD3Mode.SURFACE_NORMAL.value
100+
101+
for key in [
102+
"Refine number of elements"
103+
]:
104+
if options[key] < 1:
105+
options[key] = 1
106+
107+
return dependent_changes
108+
109+
@classmethod
110+
def generateBaseMesh(cls, region, options):
111+
"""
112+
Generate the base tricubic Hermite mesh. See also generateMesh().
113+
:param region: Zinc region to define model in. Must be empty.
114+
:param options: Dict containing options. See getDefaultOptions().
115+
:return: empty list of AnnotationGroup, None
116+
"""
117+
element_counts = [options[key] for key in [
118+
"Number of elements across axis 1", "Number of elements across axis 2", "Number of elements across axis 3"]]
119+
transition_element_count = options["Number of transition elements"]
120+
a = options["Axis length x"]
121+
b = options["Axis length y"]
122+
c = options["Axis length z"]
123+
axis2_x_rotation_radians = math.radians(options["Axis 2 x-rotation degrees"])
124+
axis3_x_rotation_radians = math.radians(options["Axis 3 x-rotation degrees"])
125+
surface_only = options["2D surface only"]
126+
nway_derivative_factor = options["Advanced n-way derivative factor"]
127+
surface_d3_mode = EllipsoidSurfaceD3Mode(options["Advanced surface D3 mode"])
128+
129+
fieldmodule = region.getFieldmodule()
130+
coordinates = find_or_create_field_coordinates(fieldmodule)
131+
132+
ellipsoid = EllipsoidMesh(a, b, c, element_counts, transition_element_count,
133+
axis2_x_rotation_radians, axis3_x_rotation_radians, surface_only)
134+
135+
left_group = AnnotationGroup(region, ("left", ""))
136+
right_group = AnnotationGroup(region, ("right", ""))
137+
back_group = AnnotationGroup(region, ("back", ""))
138+
front_group = AnnotationGroup(region, ("front", ""))
139+
bottom_group = AnnotationGroup(region, ("bottom", ""))
140+
top_group = AnnotationGroup(region, ("top", ""))
141+
annotation_groups = [left_group, right_group, back_group, front_group, bottom_group, top_group]
142+
octant_group_lists = []
143+
for octant in range(8):
144+
octant_group_list = []
145+
octant_group_list.append((right_group if (octant & 1) else left_group).getGroup())
146+
octant_group_list.append((front_group if (octant & 2) else back_group).getGroup())
147+
octant_group_list.append((top_group if (octant & 4) else bottom_group).getGroup())
148+
octant_group_lists.append(octant_group_list)
149+
ellipsoid.set_octant_group_lists(octant_group_lists)
150+
151+
if not surface_only:
152+
box_group = AnnotationGroup(region, ("box", ""))
153+
transition_group = AnnotationGroup(region, ("transition", ""))
154+
annotation_groups += [box_group, transition_group]
155+
ellipsoid.set_box_transition_groups(box_group.getGroup(), transition_group.getGroup())
156+
157+
ellipsoid.set_nway_derivative_factor(nway_derivative_factor)
158+
ellipsoid.set_surface_d3_mode(surface_d3_mode)
159+
160+
ellipsoid.build()
161+
ellipsoid.generate_mesh(fieldmodule, coordinates)
162+
163+
return annotation_groups, None
164+
165+
@classmethod
166+
def refineMesh(cls, meshRefinement, options):
167+
"""
168+
Refine source mesh into separate region, with change of basis.
169+
:param meshRefinement: MeshRefinement, which knows source and target region.
170+
:param options: Dict containing options. See getDefaultOptions().
171+
"""
172+
assert isinstance(meshRefinement, MeshRefinement)
173+
refineElementsCount = options["Refine number of elements"]
174+
meshRefinement.refineAllElementsCubeStandard3d(refineElementsCount, refineElementsCount, refineElementsCount)

src/scaffoldmaker/scaffolds.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
from scaffoldmaker.meshtypes.meshtype_3d_cecum1 import MeshType_3d_cecum1
2323
from scaffoldmaker.meshtypes.meshtype_3d_colon1 import MeshType_3d_colon1
2424
from scaffoldmaker.meshtypes.meshtype_3d_colonsegment1 import MeshType_3d_colonsegment1
25+
from scaffoldmaker.meshtypes.meshtype_3d_ellipsoid1 import MeshType_3d_ellipsoid1
2526
from scaffoldmaker.meshtypes.meshtype_3d_esophagus1 import MeshType_3d_esophagus1
2627
from scaffoldmaker.meshtypes.meshtype_3d_gastrointestinaltract1 import MeshType_3d_gastrointestinaltract1
2728
from scaffoldmaker.meshtypes.meshtype_3d_heart1 import MeshType_3d_heart1
@@ -88,6 +89,7 @@ def __init__(self):
8889
MeshType_3d_cecum1,
8990
MeshType_3d_colon1,
9091
MeshType_3d_colonsegment1,
92+
MeshType_3d_ellipsoid1,
9193
MeshType_3d_esophagus1,
9294
MeshType_3d_gastrointestinaltract1,
9395
MeshType_3d_heart1,

src/scaffoldmaker/utils/eft_utils.py

Lines changed: 70 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -687,6 +687,38 @@ def __init__(self):
687687
[[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 1.0, 0.0],
688688
[-1.0, 0.0, 0.0], [-1.0, -1.0, 0.0], [0.0, -1.0, 0.0], [1.0, -1.0, 0.0],
689689
[0.0, 0.0, -1.0], [0.0, 0.0, 1.0]])
690+
self._nodeLayoutTriplePoint2DQ1 = HermiteNodeLayout(
691+
[[-1.0, 0.0], [0.0, -1.0], [1.0, 1.0]])
692+
self._nodeLayoutTriplePoint2DQ2 = HermiteNodeLayout(
693+
[[0.0, -1.0], [1.0, 0.0], [-1.0, 1.0]])
694+
self._nodeLayoutTriplePoint2DQ3 = HermiteNodeLayout(
695+
[[1.0, 0.0], [0.0, 1.0], [-1.0, -1.0]])
696+
self._nodeLayoutTriplePoint2DQ4 = HermiteNodeLayout(
697+
[[0.0, 1.0], [-1.0, 0.0], [1.0, -1.0]])
698+
self._nodeLayout3WayPoints12 = [
699+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [-1.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]),
700+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]),
701+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]]),
702+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 0.0, 1.0]])]
703+
self._nodeLayout3WayPoints13 = [
704+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
705+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
706+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]),
707+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 0.0], [0.0, 1.0, 0.0]])]
708+
self._nodeLayout3WayPoints23 = [
709+
HermiteNodeLayout([[0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
710+
HermiteNodeLayout([[0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
711+
HermiteNodeLayout([[0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [0.0, -1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]),
712+
HermiteNodeLayout([[0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 1.0], [-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])]
713+
self._nodeLayout4WayPoints = [
714+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]),
715+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, -1.0, 1.0]]),
716+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]]),
717+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]),
718+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, -1.0, 1.0]]),
719+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0], [1.0, -1.0, 1.0]]),
720+
HermiteNodeLayout([[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 1.0, 1.0]]),
721+
HermiteNodeLayout([[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [1.0, 1.0, 1.0]])]
690722
self._nodeLayoutTriplePointTopLeft = HermiteNodeLayout(
691723
[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 1.0]])
692724
self._nodeLayoutTriplePointTopRight = HermiteNodeLayout(
@@ -772,11 +804,48 @@ def getNodeLayoutTriplePoint23Back(self):
772804

773805
return nodeLayout
774806

807+
def getNodeLayoutTriplePoint2D(self):
808+
"""
809+
Get node layout for triple-point corners of 2D quadrants.
810+
:return: List of 4 HermiteNodeLayout.
811+
"""
812+
nodeLayouts = [self._nodeLayoutTriplePoint2DQ1, self._nodeLayoutTriplePoint2DQ2,
813+
self._nodeLayoutTriplePoint2DQ3, self._nodeLayoutTriplePoint2DQ4]
814+
return nodeLayouts
815+
816+
def getNodeLayout3WayPoints12(self):
817+
"""
818+
Get 3-way node layouts for quadrants 12 = NN, NP, PN, PP.
819+
:return: List of 4 HermiteNodeLayout.
820+
"""
821+
return self._nodeLayout3WayPoints12
822+
823+
def getNodeLayout3WayPoints13(self):
824+
"""
825+
Get 3-way node layouts for quadrants 13 = NN, NP, PN, PP.
826+
:return: List of 4 HermiteNodeLayout.
827+
"""
828+
return self._nodeLayout3WayPoints13
829+
830+
def getNodeLayout3WayPoints23(self):
831+
"""
832+
Get 3-way node layouts for quadrants 23 = NN, NP, PN, PP.
833+
:return: List of 4 HermiteNodeLayout.
834+
"""
835+
return self._nodeLayout3WayPoints23
836+
837+
def getNodeLayout4WayPoints(self):
838+
"""
839+
Get node layouts from a regular core for octants 123: NNN, NNP, NPN, NPP, PNN, PNP, PPN, PPP.
840+
:return: List of 8 HermiteNodeLayout.
841+
"""
842+
return self._nodeLayout4WayPoints
843+
775844
def getNodeLayoutTriplePoint(self):
776845
"""
777846
Get node layout for triple-point corners of core box elements. There are four corners (Top Left, Top Right,
778847
Bottom Left, and Bottom Right) each with its specific node layout.
779-
:return: HermiteNodeLayout.
848+
:return: List of 4 HermiteNodeLayout.
780849
"""
781850
nodeLayouts = [self._nodeLayoutTriplePointTopLeft, self._nodeLayoutTriplePointTopRight,
782851
self._nodeLayoutTriplePointBottomLeft, self._nodeLayoutTriplePointBottomRight]

0 commit comments

Comments
 (0)