v0.15.0
Loading...
Searching...
No Matches
EntityRefine.cpp
Go to the documentation of this file.
1/* \file EntityRefine.cpp
2 * \brief Tetrahedral refinement algorithm
3
4 It is based on \cite ruprecht1998scheme
5
6 \todo tet refinement should be rewritten, with better error control, and more
7 clear refinement patterns.
8
9*/
10
11
12
13// A scheme for Edge-based Adaptive Tetrahedral Subdivision; Delft Ruprecht
14
15namespace MoFEM {
16
17// TET
18static constexpr int edges_conn[] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
19static constexpr int oposite_edge[] = {5, 3, 4, 1, 2, 0};
20static constexpr int edge_permutations[6][6] = {
21 {0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4},
22 {3, 4, 0, 2, 5, 1}, {4, 5, 1, 0, 3, 2}, {5, 3, 2, 1, 4, 0}};
23static constexpr int edge_mirror_cross[6] = {0, 3, 4, 1, 2, 5};
24static constexpr int edge_mirror_vertical[6] = {0, 4, 3, 2, 1, 5};
25static constexpr int cyclic_node_rotate_face_3[3][4] = {
26 {3, 1, 0, 2}, {0, 1, 2, 3}, {2, 1, 3, 0}}; // 2,0,1
27static constexpr int cyclic_edge_rotate_face_3[3][6] = {
28 {4, 0, 3, 5, 1, 2}, {0, 1, 2, 3, 4, 5}, {1, 4, 5, 2, 0, 3}};
29static constexpr char edge_bits_mark[] = {1, 2, 4, 8, 16, 32};
30
31void tet_type_6(moab::Interface &moab, const EntityHandle *conn,
32 const EntityHandle *edge_new_nodes,
33 EntityHandle *new_tets_conn) {
34 // 0:01 - 4 1:12-5 2:20-6 3:03-7 4:13-8 5:23-9
35 // TET0
36 new_tets_conn[0 * 4 + 0] = conn[0];
37 new_tets_conn[0 * 4 + 1] = edge_new_nodes[0];
38 new_tets_conn[0 * 4 + 2] = edge_new_nodes[2];
39 new_tets_conn[0 * 4 + 3] = edge_new_nodes[3];
40 // TET1
41 new_tets_conn[1 * 4 + 0] = conn[1];
42 new_tets_conn[1 * 4 + 1] = edge_new_nodes[0];
43 new_tets_conn[1 * 4 + 2] = edge_new_nodes[4];
44 new_tets_conn[1 * 4 + 3] = edge_new_nodes[1];
45 // TET2
46 new_tets_conn[2 * 4 + 0] = conn[2];
47 new_tets_conn[2 * 4 + 1] = edge_new_nodes[1];
48 new_tets_conn[2 * 4 + 2] = edge_new_nodes[5];
49 new_tets_conn[2 * 4 + 3] = edge_new_nodes[2];
50 // TET3
51 new_tets_conn[3 * 4 + 0] = conn[3];
52 new_tets_conn[3 * 4 + 1] = edge_new_nodes[3];
53 new_tets_conn[3 * 4 + 2] = edge_new_nodes[5];
54 new_tets_conn[3 * 4 + 3] = edge_new_nodes[4];
55 double coords[6 * 3];
56 moab.get_coords(edge_new_nodes, 6, coords);
57 cblas_daxpy(3, -1, &coords[4 * 3], 1, &coords[2 * 3], 1);
58 cblas_daxpy(3, -1, &coords[3 * 3], 1, &coords[1 * 3], 1);
59 cblas_daxpy(3, -1, &coords[5 * 3], 1, &coords[0 * 3], 1);
60 double L[3] = {cblas_dnrm2(3, &coords[2 * 3], 1),
61 cblas_dnrm2(3, &coords[1 * 3], 1),
62 cblas_dnrm2(3, &coords[0 * 3], 1)};
63 // VARIANT 1 - diag 4-2
64 if (L[0] <= L[1] && L[0] <= L[2]) {
65 // TET4
66 new_tets_conn[4 * 4 + 0] = edge_new_nodes[4];
67 new_tets_conn[4 * 4 + 1] = edge_new_nodes[3];
68 new_tets_conn[4 * 4 + 2] = edge_new_nodes[2];
69 new_tets_conn[4 * 4 + 3] = edge_new_nodes[0];
70 // TET5
71 new_tets_conn[5 * 4 + 0] = edge_new_nodes[4];
72 new_tets_conn[5 * 4 + 1] = edge_new_nodes[2];
73 new_tets_conn[5 * 4 + 2] = edge_new_nodes[3];
74 new_tets_conn[5 * 4 + 3] = edge_new_nodes[5];
75 // TET6
76 new_tets_conn[6 * 4 + 0] = edge_new_nodes[4];
77 new_tets_conn[6 * 4 + 1] = edge_new_nodes[2];
78 new_tets_conn[6 * 4 + 2] = edge_new_nodes[1];
79 new_tets_conn[6 * 4 + 3] = edge_new_nodes[0];
80 // TET7
81 new_tets_conn[7 * 4 + 0] = edge_new_nodes[4];
82 new_tets_conn[7 * 4 + 1] = edge_new_nodes[1];
83 new_tets_conn[7 * 4 + 2] = edge_new_nodes[2];
84 new_tets_conn[7 * 4 + 3] = edge_new_nodes[5];
85 return;
86 }
87 // VARIANT 2 - diag 3-1
88 if (L[1] <= L[0] && L[1] <= L[2]) {
89 // TET4
90 new_tets_conn[4 * 4 + 0] = edge_new_nodes[4];
91 new_tets_conn[4 * 4 + 1] = edge_new_nodes[3];
92 new_tets_conn[4 * 4 + 2] = edge_new_nodes[1];
93 new_tets_conn[4 * 4 + 3] = edge_new_nodes[0];
94 // TET5
95 new_tets_conn[5 * 4 + 0] = edge_new_nodes[4];
96 new_tets_conn[5 * 4 + 1] = edge_new_nodes[1];
97 new_tets_conn[5 * 4 + 2] = edge_new_nodes[3];
98 new_tets_conn[5 * 4 + 3] = edge_new_nodes[5];
99 // TET6
100 new_tets_conn[6 * 4 + 0] = edge_new_nodes[1];
101 new_tets_conn[6 * 4 + 1] = edge_new_nodes[3];
102 new_tets_conn[6 * 4 + 2] = edge_new_nodes[2];
103 new_tets_conn[6 * 4 + 3] = edge_new_nodes[0];
104 // TET7
105 new_tets_conn[7 * 4 + 0] = edge_new_nodes[1];
106 new_tets_conn[7 * 4 + 1] = edge_new_nodes[2];
107 new_tets_conn[7 * 4 + 2] = edge_new_nodes[3];
108 new_tets_conn[7 * 4 + 3] = edge_new_nodes[5];
109 return;
110 }
111 // VARIANT 3 - diag 5-0
112 // TET4
113 new_tets_conn[4 * 4 + 0] = edge_new_nodes[5];
114 new_tets_conn[4 * 4 + 1] = edge_new_nodes[2];
115 new_tets_conn[4 * 4 + 2] = edge_new_nodes[0];
116 new_tets_conn[4 * 4 + 3] = edge_new_nodes[3];
117 // TET5
118 new_tets_conn[5 * 4 + 0] = edge_new_nodes[5];
119 new_tets_conn[5 * 4 + 1] = edge_new_nodes[0];
120 new_tets_conn[5 * 4 + 2] = edge_new_nodes[2];
121 new_tets_conn[5 * 4 + 3] = edge_new_nodes[1];
122 // TET6
123 new_tets_conn[6 * 4 + 0] = edge_new_nodes[5];
124 new_tets_conn[6 * 4 + 1] = edge_new_nodes[0];
125 new_tets_conn[6 * 4 + 2] = edge_new_nodes[4];
126 new_tets_conn[6 * 4 + 3] = edge_new_nodes[3];
127 // TET7
128 new_tets_conn[7 * 4 + 0] = edge_new_nodes[5];
129 new_tets_conn[7 * 4 + 1] = edge_new_nodes[4];
130 new_tets_conn[7 * 4 + 2] = edge_new_nodes[0];
131 new_tets_conn[7 * 4 + 3] = edge_new_nodes[1];
132}
133int tet_type_5(moab::Interface &moab, const EntityHandle *conn,
134 const EntityHandle *edge_new_nodes,
135 EntityHandle *new_tets_conn) {
136 int free_edge = -1;
137 for (int ee = 0; ee < 6; ee++) {
138 if (edge_new_nodes[ee] == no_handle) {
139 free_edge = ee;
140 break;
141 }
142 }
143 int edge0 = oposite_edge[free_edge];
144 EntityHandle conn_[] = {
145 conn[edges_conn[edge0 * 2 + 0]], conn[edges_conn[edge0 * 2 + 1]],
146 conn[edges_conn[free_edge * 2 + 0]], conn[edges_conn[free_edge * 2 + 1]]};
147 const int *edges_ = &edge_permutations[edge0][0];
148 EntityHandle edge_new_nodes_[6];
149 for (int ee = 0; ee < 6; ee++)
150 edge_new_nodes_[ee] = edge_new_nodes[edges_[ee]];
151 bool free_edge_swappped = false;
152 if (conn_[3] < conn_[2]) {
153 free_edge_swappped = true;
154 EntityHandle conn__2_ = conn_[2];
155 conn_[2] = conn_[3];
156 conn_[3] = conn__2_;
157 }
158 assert(conn_[0] != no_handle);
159 assert(conn_[1] != no_handle);
160 assert(conn_[2] != no_handle);
161 assert(conn_[3] != no_handle);
162 assert(edge_new_nodes_[0] != no_handle);
163 assert(edge_new_nodes_[1] != no_handle);
164 assert(edge_new_nodes_[2] != no_handle);
165 assert(edge_new_nodes_[3] != no_handle);
166 assert(edge_new_nodes_[4] != no_handle);
167 // TET0
168 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
169 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[0];
170 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[1];
171 new_tets_conn[0 * 4 + 3] = conn_[1];
172 // TET1
173 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[2];
174 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[0];
175 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[3];
176 new_tets_conn[1 * 4 + 3] = conn_[0];
177 // TET4
178 new_tets_conn[2 * 4 + 0] = conn_[2];
179 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[3];
180 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[4];
181 if (free_edge_swappped) {
182 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
183 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[2];
184 }
185 new_tets_conn[2 * 4 + 3] = conn_[3];
186 double coords[6 * 3];
187 moab.get_coords(edge_new_nodes_, 6, coords);
188 cblas_daxpy(3, -1, &coords[4 * 3], 1, &coords[2 * 3], 1);
189 cblas_daxpy(3, -1, &coords[3 * 3], 1, &coords[1 * 3], 1);
190 double L[2] = {cblas_dnrm2(3, &coords[2 * 3], 1),
191 cblas_dnrm2(3, &coords[1 * 3], 1)};
192 if (L[1] <= L[0]) {
193 // VARIANT 1 diag 4-2
194 // TET2
195 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
196 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[3];
197 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[1];
198 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[0];
199 // TET3
200 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[2];
201 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[1];
202 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[3];
203 new_tets_conn[4 * 4 + 3] = edge_new_nodes_[0];
204 // TET5
205 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[4];
206 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[1];
207 new_tets_conn[5 * 4 + 2] = edge_new_nodes_[3];
208 new_tets_conn[5 * 4 + 3] = conn_[2];
209 // TET6
210 new_tets_conn[6 * 4 + 0] = edge_new_nodes_[1];
211 new_tets_conn[6 * 4 + 1] = edge_new_nodes_[2];
212 new_tets_conn[6 * 4 + 2] = edge_new_nodes_[3];
213 new_tets_conn[6 * 4 + 3] = conn_[2];
214 return 1;
215 }
216 // VARIANT 2 diag 1-3
217 // TET2
218 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
219 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[3];
220 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[2];
221 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[0];
222 // TET3
223 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[2];
224 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[1];
225 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
226 new_tets_conn[4 * 4 + 3] = edge_new_nodes_[0];
227 // TET5
228 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[4];
229 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[1];
230 new_tets_conn[5 * 4 + 2] = edge_new_nodes_[2];
231 new_tets_conn[5 * 4 + 3] = conn_[2];
232 // TET6
233 new_tets_conn[6 * 4 + 0] = edge_new_nodes_[4];
234 new_tets_conn[6 * 4 + 1] = edge_new_nodes_[2];
235 new_tets_conn[6 * 4 + 2] = edge_new_nodes_[3];
236 new_tets_conn[6 * 4 + 3] = conn_[2];
237 return 0;
238}
239int tet_type_4(const EntityHandle *conn, const int *split_edges,
240 const EntityHandle *edge_new_nodes,
241 EntityHandle *new_tets_conn) {
242 char mach_pattern = 0;
243 for (int ee = 0; ee < 4; ee++)
244 mach_pattern |= edge_bits_mark[split_edges[ee]];
245 EntityHandle conn_[4] = {no_handle, no_handle, no_handle, no_handle};
246 int edges_[6];
247 int type = -1;
248 for (int ee = 0; ee < 6; ee++) {
249 char pattern0, pattern1;
250 pattern0 = edge_bits_mark[edge_permutations[ee][0]] |
251 edge_bits_mark[edge_permutations[ee][1]] |
252 edge_bits_mark[edge_permutations[ee][4]] |
253 edge_bits_mark[edge_permutations[ee][3]];
254 pattern1 = edge_bits_mark[edge_permutations[ee][1]] |
255 edge_bits_mark[edge_permutations[ee][4]] |
256 edge_bits_mark[edge_permutations[ee][2]] |
257 edge_bits_mark[edge_permutations[ee][3]];
258 if (pattern0 == mach_pattern || pattern1 == mach_pattern) {
259 int free_edge = oposite_edge[ee];
260 conn_[0] = conn[edges_conn[ee * 2 + 0]];
261 conn_[1] = conn[edges_conn[ee * 2 + 1]];
262 conn_[2] = conn[edges_conn[free_edge * 2 + 0]];
263 conn_[3] = conn[edges_conn[free_edge * 2 + 1]];
264 for (int EE = 0; EE < 6; EE++)
265 edges_[EE] = edge_permutations[ee][EE];
266 if (pattern0 == mach_pattern)
267 type = 0;
268 else if (pattern1 == mach_pattern)
269 type = 1;
270 // printf("no mirror\n");
271 break;
272 }
273 pattern0 = edge_bits_mark[edge_permutations[ee][edge_mirror_cross[0]]] |
274 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[1]]] |
275 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[4]]] |
276 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[3]]];
277 pattern1 = edge_bits_mark[edge_permutations[ee][1]] |
278 edge_bits_mark[edge_permutations[ee][4]] |
279 edge_bits_mark[edge_permutations[ee][2]] |
280 edge_bits_mark[edge_permutations[ee][3]];
281 if (pattern0 == mach_pattern || pattern1 == mach_pattern) {
282 int free_edge = oposite_edge[ee];
283 conn_[0] = conn[edges_conn[ee * 2 + 1]];
284 conn_[1] = conn[edges_conn[ee * 2 + 0]];
285 conn_[2] = conn[edges_conn[free_edge * 2 + 1]];
286 conn_[3] = conn[edges_conn[free_edge * 2 + 0]];
287 for (int EE = 0; EE < 6; EE++)
288 edges_[EE] = edge_permutations[ee][edge_mirror_cross[EE]];
289 if (pattern0 == mach_pattern)
290 type = 0;
291 else if (pattern1 == mach_pattern)
292 type = 1;
293 // printf("mirror\n");
294 break;
295 }
296 }
297 assert(type != -1);
298 EntityHandle edge_new_nodes_[6];
299 for (int ee = 0; ee < 6; ee++)
300 edge_new_nodes_[ee] = edge_new_nodes[edges_[ee]];
301 if (type == 0) {
302 assert(edge_new_nodes_[0] != no_handle);
303 assert(edge_new_nodes_[1] != no_handle);
304 assert(edge_new_nodes_[4] != no_handle);
305 assert(edge_new_nodes_[3] != no_handle);
306 bool free_edge_swappped5 = false;
307 if (conn_[3] < conn_[2]) {
308 free_edge_swappped5 = true;
309 }
310 bool free_edge_swappped2 = false;
311 if (conn_[0] < conn_[2]) {
312 free_edge_swappped2 = true;
313 }
314 // TET0
315 new_tets_conn[0 * 4 + 0] = conn_[1];
316 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
317 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[0];
318 new_tets_conn[0 * 4 + 3] = edge_new_nodes_[4];
319 if (free_edge_swappped5 && (!free_edge_swappped2)) {
320 // TET1
321 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
322 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[3];
323 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
324 new_tets_conn[1 * 4 + 3] = conn_[3];
325 // TET2
326 new_tets_conn[2 * 4 + 0] = conn_[3];
327 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
328 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[3];
329 new_tets_conn[2 * 4 + 3] = conn_[2];
330 // TET3
331 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[0];
332 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[3];
333 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[4];
334 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[1];
335 // TET4
336 new_tets_conn[4 * 4 + 0] = conn_[2];
337 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[0];
338 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[3];
339 new_tets_conn[4 * 4 + 3] = conn_[0];
340 // TET5
341 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[1];
342 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[0];
343 new_tets_conn[5 * 4 + 2] = edge_new_nodes_[3];
344 new_tets_conn[5 * 4 + 3] = conn_[2];
345 return 2;
346 } else if (free_edge_swappped2 && (!free_edge_swappped5)) {
347 // TET1
348 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
349 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
350 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[0];
351 new_tets_conn[1 * 4 + 3] = conn_[0];
352 // TET2
353 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[3];
354 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
355 new_tets_conn[2 * 4 + 2] = conn_[0];
356 new_tets_conn[2 * 4 + 3] = conn_[2];
357 // TET3
358 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[0];
359 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[3];
360 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[4];
361 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[1];
362 // TET4
363 new_tets_conn[4 * 4 + 0] = conn_[2];
364 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
365 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
366 new_tets_conn[4 * 4 + 3] = conn_[3];
367 // TET5
368 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[1];
369 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[3];
370 new_tets_conn[5 * 4 + 2] = edge_new_nodes_[4];
371 new_tets_conn[5 * 4 + 3] = conn_[2];
372 return 3;
373 } else if (free_edge_swappped2 && free_edge_swappped5) {
374 // TET1
375 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
376 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
377 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[0];
378 new_tets_conn[1 * 4 + 3] = conn_[0];
379 // TET2
380 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[3];
381 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
382 new_tets_conn[2 * 4 + 2] = conn_[0];
383 new_tets_conn[2 * 4 + 3] = conn_[2];
384 // TET3
385 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[0];
386 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[3];
387 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[4];
388 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[1];
389 // TET4
390 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[1];
391 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
392 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
393 new_tets_conn[4 * 4 + 3] = conn_[3];
394 // TET5
395 new_tets_conn[5 * 4 + 0] = conn_[3];
396 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[1];
397 new_tets_conn[5 * 4 + 2] = edge_new_nodes_[3];
398 new_tets_conn[5 * 4 + 3] = conn_[2];
399 return 4;
400 } else {
401 // TET1
402 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[0];
403 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[3];
404 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
405 new_tets_conn[1 * 4 + 3] = conn_[2];
406 // TET2
407 new_tets_conn[2 * 4 + 0] = conn_[2];
408 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
409 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[4];
410 new_tets_conn[2 * 4 + 3] = edge_new_nodes_[0];
411 // TET3
412 new_tets_conn[3 * 4 + 0] = conn_[2];
413 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[0];
414 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[3];
415 new_tets_conn[3 * 4 + 3] = conn_[0];
416 // TET4
417 new_tets_conn[4 * 4 + 0] = conn_[2];
418 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
419 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
420 new_tets_conn[4 * 4 + 3] = conn_[3];
421 }
422 } else if (type == 1) {
423 assert(edge_new_nodes_[1] != no_handle);
424 assert(edge_new_nodes_[2] != no_handle);
425 assert(edge_new_nodes_[3] != no_handle);
426 assert(edge_new_nodes_[4] != no_handle);
427 bool free_edge_swappped5 = false;
428 if (conn_[3] < conn_[2]) {
429 free_edge_swappped5 = true;
430 }
431 bool free_edge_swappped0 = false;
432 if (conn_[0] > conn_[1]) {
433 free_edge_swappped0 = true;
434 }
435 if (free_edge_swappped0 && (!free_edge_swappped5)) {
436 // TET0
437 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[3];
438 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[2];
439 new_tets_conn[0 * 4 + 2] = conn_[1];
440 new_tets_conn[0 * 4 + 3] = conn_[0];
441 // TET1
442 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
443 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
444 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
445 new_tets_conn[1 * 4 + 3] = conn_[1];
446 // TET2
447 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[3];
448 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
449 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[1];
450 new_tets_conn[2 * 4 + 3] = conn_[1];
451 // TET3
452 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[1];
453 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
454 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[3];
455 new_tets_conn[3 * 4 + 3] = conn_[2];
456 // TET4
457 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[1];
458 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
459 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
460 new_tets_conn[4 * 4 + 3] = conn_[2];
461 // TET5
462 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[3];
463 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[4];
464 new_tets_conn[5 * 4 + 2] = conn_[2];
465 new_tets_conn[5 * 4 + 3] = conn_[3];
466 return 5;
467 } else if (free_edge_swappped5 && (!free_edge_swappped0)) {
468 // TET0
469 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
470 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
471 new_tets_conn[0 * 4 + 2] = conn_[1];
472 new_tets_conn[0 * 4 + 3] = conn_[0];
473 // TET1
474 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
475 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
476 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
477 new_tets_conn[1 * 4 + 3] = conn_[0];
478 // TET2
479 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[2];
480 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
481 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[3];
482 new_tets_conn[2 * 4 + 3] = conn_[0];
483 // TET3
484 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[1];
485 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
486 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[3];
487 new_tets_conn[3 * 4 + 3] = conn_[3];
488 // TET4
489 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[1];
490 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
491 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
492 new_tets_conn[4 * 4 + 3] = conn_[3];
493 // TET5
494 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[1];
495 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[2];
496 new_tets_conn[5 * 4 + 2] = conn_[3];
497 new_tets_conn[5 * 4 + 3] = conn_[2];
498 return 6;
499 } else if (free_edge_swappped5 && free_edge_swappped0) {
500 // TET0
501 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[3];
502 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[2];
503 new_tets_conn[0 * 4 + 2] = conn_[1];
504 new_tets_conn[0 * 4 + 3] = conn_[0];
505 // TET1
506 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
507 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
508 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
509 new_tets_conn[1 * 4 + 3] = conn_[1];
510 // TET2
511 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[3];
512 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
513 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[1];
514 new_tets_conn[2 * 4 + 3] = conn_[1];
515 // TET3
516 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[1];
517 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
518 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[3];
519 new_tets_conn[3 * 4 + 3] = conn_[3];
520 // TET4
521 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[1];
522 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
523 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
524 new_tets_conn[4 * 4 + 3] = conn_[3];
525 // TET5
526 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[1];
527 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[2];
528 new_tets_conn[5 * 4 + 2] = conn_[3];
529 new_tets_conn[5 * 4 + 3] = conn_[2];
530 return 7;
531 }
532 // TET0
533 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
534 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
535 new_tets_conn[0 * 4 + 2] = conn_[1];
536 new_tets_conn[0 * 4 + 3] = conn_[0];
537 // TET1
538 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[3];
539 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
540 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
541 new_tets_conn[1 * 4 + 3] = conn_[0];
542 // TET2
543 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[2];
544 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
545 new_tets_conn[2 * 4 + 2] = edge_new_nodes_[3];
546 new_tets_conn[2 * 4 + 3] = conn_[0];
547 // TET3
548 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[1];
549 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
550 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[3];
551 new_tets_conn[3 * 4 + 3] = conn_[2];
552 // TET4
553 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[1];
554 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[3];
555 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[4];
556 new_tets_conn[4 * 4 + 3] = conn_[2];
557 // TET5
558 new_tets_conn[5 * 4 + 0] = edge_new_nodes_[3];
559 new_tets_conn[5 * 4 + 1] = edge_new_nodes_[4];
560 new_tets_conn[5 * 4 + 2] = conn_[2];
561 new_tets_conn[5 * 4 + 3] = conn_[3];
562 }
563 return type;
564}
565int tet_type_3(const EntityHandle *conn, const int *split_edges,
566 const EntityHandle *edge_new_nodes,
567 EntityHandle *new_tets_conn) {
568 char mach_pattern = 0;
569 for (int ee = 0; ee < 3; ee++)
570 mach_pattern |= edge_bits_mark[split_edges[ee]];
571 EntityHandle conn_[4];
572 int edges_[6];
573 int type = -1;
574 bool is_rotated = false;
575 for (int ee = 0; ee < 6; ee++) {
576 char pattern0, pattern1, pattern2;
577 pattern0 = edge_bits_mark[edge_permutations[ee][0]] |
578 edge_bits_mark[edge_permutations[ee][1]] |
579 edge_bits_mark[edge_permutations[ee][4]];
580 pattern1 = edge_bits_mark[edge_permutations[ee][1]] |
581 edge_bits_mark[edge_permutations[ee][4]] |
582 edge_bits_mark[edge_permutations[ee][5]];
583 pattern2 = edge_bits_mark[edge_permutations[ee][4]] |
584 edge_bits_mark[edge_permutations[ee][1]] |
585 edge_bits_mark[edge_permutations[ee][2]];
586 if (pattern0 == mach_pattern || pattern1 == mach_pattern ||
587 pattern2 == mach_pattern) {
588 // printf("nothing\n");
589 int free_edge = oposite_edge[ee];
590 conn_[0] = conn[edges_conn[ee * 2 + 0]];
591 conn_[1] = conn[edges_conn[ee * 2 + 1]];
592 conn_[2] = conn[edges_conn[free_edge * 2 + 0]];
593 conn_[3] = conn[edges_conn[free_edge * 2 + 1]];
594 for (int EE = 0; EE < 6; EE++)
595 edges_[EE] = edge_permutations[ee][EE];
596 if (pattern0 == mach_pattern) {
597 // printf("nothing\n");
598 type = 0;
599 } else if (pattern1 == mach_pattern)
600 type = 1;
601 else if (pattern2 == mach_pattern)
602 type = 2;
603 break;
604 }
605 pattern0 = edge_bits_mark[edge_permutations[ee][edge_mirror_cross[0]]] |
606 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[1]]] |
607 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[4]]];
608 pattern1 = edge_bits_mark[edge_permutations[ee][edge_mirror_cross[1]]] |
609 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[4]]] |
610 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[5]]];
611 pattern2 = edge_bits_mark[edge_permutations[ee][edge_mirror_cross[4]]] |
612 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[1]]] |
613 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[2]]];
614 if (pattern0 == mach_pattern || pattern1 == mach_pattern ||
615 pattern2 == mach_pattern) {
616 // printf("edge_mirror_cross\n");
617 int free_edge = oposite_edge[ee];
618 conn_[0] = conn[edges_conn[ee * 2 + 1]];
619 conn_[1] = conn[edges_conn[ee * 2 + 0]];
620 conn_[2] = conn[edges_conn[free_edge * 2 + 1]];
621 conn_[3] = conn[edges_conn[free_edge * 2 + 0]];
622 for (int EE = 0; EE < 6; EE++)
623 edges_[EE] = edge_permutations[ee][edge_mirror_cross[EE]];
624 if (pattern0 == mach_pattern) {
625 // printf("edge_mirror_cross\n");
626 type = 0;
627 } else if (pattern1 == mach_pattern)
628 type = 1;
629 else if (pattern2 == mach_pattern)
630 type = 2;
631 break;
632 }
633 pattern2 = edge_bits_mark[edge_permutations[ee][edge_mirror_vertical[4]]] |
634 edge_bits_mark[edge_permutations[ee][edge_mirror_vertical[1]]] |
635 edge_bits_mark[edge_permutations[ee][edge_mirror_vertical[2]]];
636 if (pattern2 == mach_pattern) {
637 is_rotated = true;
638 // printf("edge_mirror_vertical\n");
639 int free_edge = oposite_edge[ee];
640 conn_[0] = conn[edges_conn[ee * 2 + 0]];
641 conn_[1] = conn[edges_conn[ee * 2 + 1]];
642 conn_[2] = conn[edges_conn[free_edge * 2 + 1]];
643 conn_[3] = conn[edges_conn[free_edge * 2 + 0]];
644 for (int EE = 0; EE < 6; EE++)
645 edges_[EE] = edge_permutations[ee][edge_mirror_vertical[EE]];
646 if (pattern0 == mach_pattern) {
647 // printf("edge_mirror_vertical\n");
648 type = 0;
649 } else if (pattern1 == mach_pattern)
650 type = 1;
651 else if (pattern2 == mach_pattern)
652 type = 2;
653 break;
654 }
655 pattern2 =
656 edge_bits_mark[edge_permutations
657 [ee][edge_mirror_cross[edge_mirror_vertical[4]]]] |
658 edge_bits_mark[edge_permutations
659 [ee][edge_mirror_cross[edge_mirror_vertical[1]]]] |
660 edge_bits_mark
661 [edge_permutations[ee][edge_mirror_cross[edge_mirror_vertical[2]]]];
662 if (pattern2 == mach_pattern) {
663 is_rotated = true;
664 int free_edge = oposite_edge[ee];
665 conn_[0] = conn[edges_conn[ee * 2 + 1]];
666 conn_[1] = conn[edges_conn[ee * 2 + 0]];
667 conn_[2] = conn[edges_conn[free_edge * 2 + 0]];
668 conn_[3] = conn[edges_conn[free_edge * 2 + 1]];
669 for (int EE = 0; EE < 6; EE++)
670 edges_[EE] =
671 edge_permutations[ee][edge_mirror_cross[edge_mirror_vertical[EE]]];
672 if (pattern0 == mach_pattern) {
673 // printf("edge_mirror_cross|edge_mirror_vertical\n");
674 type = 0;
675 } else if (pattern1 == mach_pattern)
676 type = 1;
677 else if (pattern2 == mach_pattern)
678 type = 2;
679 break;
680 }
681 }
682 assert(type != -1);
683 EntityHandle edge_new_nodes_[6];
684 for (int ee = 0; ee < 6; ee++)
685 edge_new_nodes_[ee] = edge_new_nodes[edges_[ee]];
686 if (type == 0) {
687 EntityHandle conn__[4];
688 EntityHandle edge_new_nodes__[6];
689 bcopy(conn_, conn__, 4 * sizeof(EntityHandle));
690 bcopy(edge_new_nodes_, edge_new_nodes__, 6 * sizeof(EntityHandle));
691 for (int rotate_idx = 0; rotate_idx < 3; rotate_idx++) {
692 // fprintf(stderr,"%d\n",rotate_idx);
693 for (int ii = 0; ii < 4; ii++)
694 conn_[ii] = conn__[cyclic_node_rotate_face_3[rotate_idx][ii]];
695 for (int ee = 0; ee < 6; ee++)
696 edge_new_nodes_[ee] =
697 edge_new_nodes__[cyclic_edge_rotate_face_3[rotate_idx][ee]];
698 if ((conn_[0] > conn_[2]) && (conn_[0] > conn_[3]))
699 break;
700 }
701 assert(conn_[0] > conn_[2]);
702 assert(conn_[0] > conn_[3]);
703 assert(edge_new_nodes_[0] != no_handle);
704 assert(edge_new_nodes_[1] != no_handle);
705 assert(edge_new_nodes_[4] != no_handle);
706 // TET0
707 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[0];
708 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
709 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[4];
710 new_tets_conn[0 * 4 + 3] = conn_[1];
711 bool free_edge_swappped5 = false;
712 if (conn_[3] < conn_[2]) {
713 free_edge_swappped5 = true;
714 }
715 if (free_edge_swappped5) {
716 // TET1
717 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
718 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[0];
719 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
720 new_tets_conn[1 * 4 + 3] = conn_[3];
721 // TET2
722 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[0];
723 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[1];
724 new_tets_conn[2 * 4 + 2] = conn_[2];
725 new_tets_conn[2 * 4 + 3] = conn_[3];
726 // TET3
727 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[0];
728 new_tets_conn[3 * 4 + 1] = conn_[0];
729 new_tets_conn[3 * 4 + 2] = conn_[3];
730 new_tets_conn[3 * 4 + 3] = conn_[2];
731 // printf("free_edge_swappped5\n");
732 return 4;
733 }
734 assert(conn_[3] > conn_[2]);
735 // TET1
736 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
737 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[0];
738 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
739 new_tets_conn[1 * 4 + 3] = conn_[2];
740 // TET2
741 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[0];
742 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[4];
743 new_tets_conn[2 * 4 + 2] = conn_[2];
744 new_tets_conn[2 * 4 + 3] = conn_[3];
745 // TET3
746 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[0];
747 new_tets_conn[3 * 4 + 1] = conn_[0];
748 new_tets_conn[3 * 4 + 2] = conn_[3];
749 new_tets_conn[3 * 4 + 3] = conn_[2];
750 // printf("no free_edge_swappped5\n");
751 } else if (type == 1) {
752 assert(edge_new_nodes_[1] != no_handle);
753 assert(edge_new_nodes_[4] != no_handle);
754 assert(edge_new_nodes_[5] != no_handle);
755 // TET0
756 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[1];
757 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[4];
758 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[5];
759 new_tets_conn[0 * 4 + 3] = conn_[0];
760 // TET1
761 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[4];
762 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[5];
763 new_tets_conn[1 * 4 + 2] = conn_[0];
764 new_tets_conn[1 * 4 + 3] = conn_[3];
765 // TET2
766 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[1];
767 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[5];
768 new_tets_conn[2 * 4 + 2] = conn_[2];
769 new_tets_conn[2 * 4 + 3] = conn_[0];
770 // TET3
771 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
772 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[1];
773 new_tets_conn[3 * 4 + 2] = conn_[1];
774 new_tets_conn[3 * 4 + 3] = conn_[0];
775 } else if (type == 2) {
776 assert(edge_new_nodes_[1] != no_handle);
777 assert(edge_new_nodes_[4] != no_handle);
778 assert(edge_new_nodes_[2] != no_handle);
779 bool free_edge_swappped5 = false;
780 if (conn_[3] < conn_[2]) {
781 free_edge_swappped5 = true;
782 }
783 bool free_edge_swappped0 = false;
784 if (conn_[0] > conn_[1]) {
785 free_edge_swappped0 = true;
786 }
787 if (free_edge_swappped5 && (!free_edge_swappped0)) {
788 // TET0
789 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
790 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
791 if (is_rotated) {
792 new_tets_conn[0 * 4 + 2] = conn_[0];
793 new_tets_conn[0 * 4 + 3] = conn_[1];
794 } else {
795 new_tets_conn[0 * 4 + 2] = conn_[1];
796 new_tets_conn[0 * 4 + 3] = conn_[0];
797 }
798 // TET1
799 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
800 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[2];
801 if (is_rotated) {
802 new_tets_conn[1 * 4 + 2] = conn_[3];
803 new_tets_conn[1 * 4 + 3] = conn_[0];
804 } else {
805 new_tets_conn[1 * 4 + 2] = conn_[0];
806 new_tets_conn[1 * 4 + 3] = conn_[3];
807 }
808 // TET2
809 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[1];
810 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
811 if (is_rotated) {
812 new_tets_conn[2 * 4 + 2] = conn_[2];
813 new_tets_conn[2 * 4 + 3] = conn_[3];
814 } else {
815 new_tets_conn[2 * 4 + 2] = conn_[3];
816 new_tets_conn[2 * 4 + 3] = conn_[2];
817 }
818 // TET3
819 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
820 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[1];
821 if (is_rotated) {
822 new_tets_conn[3 * 4 + 2] = conn_[3];
823 new_tets_conn[3 * 4 + 3] = conn_[0];
824 } else {
825 new_tets_conn[3 * 4 + 2] = conn_[0];
826 new_tets_conn[3 * 4 + 3] = conn_[3];
827 }
828 return type;
829 } else if (free_edge_swappped0 && (!free_edge_swappped5)) {
830 // TET0
831 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
832 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
833 if (is_rotated) {
834 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[2];
835 new_tets_conn[0 * 4 + 3] = conn_[1];
836 } else {
837 new_tets_conn[0 * 4 + 2] = conn_[1];
838 new_tets_conn[0 * 4 + 3] = edge_new_nodes_[2];
839 }
840 // TET1
841 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[4];
842 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[1];
843 if (is_rotated) {
844 new_tets_conn[1 * 4 + 2] = conn_[2];
845 new_tets_conn[1 * 4 + 3] = edge_new_nodes_[2];
846 } else {
847 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[2];
848 new_tets_conn[1 * 4 + 3] = conn_[2];
849 }
850 // TET2
851 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[4];
852 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
853 if (is_rotated) {
854 new_tets_conn[2 * 4 + 2] = conn_[0];
855 new_tets_conn[2 * 4 + 3] = conn_[1];
856 } else {
857 new_tets_conn[2 * 4 + 2] = conn_[1];
858 new_tets_conn[2 * 4 + 3] = conn_[0];
859 }
860 // TET3
861 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
862 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
863 if (is_rotated) {
864 new_tets_conn[3 * 4 + 2] = conn_[3];
865 new_tets_conn[3 * 4 + 3] = conn_[0];
866 } else {
867 new_tets_conn[3 * 4 + 2] = conn_[0];
868 new_tets_conn[3 * 4 + 3] = conn_[3];
869 }
870 // TET4
871 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[4];
872 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[2];
873 if (is_rotated) {
874 new_tets_conn[4 * 4 + 2] = conn_[2];
875 new_tets_conn[4 * 4 + 3] = conn_[3];
876 } else {
877 new_tets_conn[4 * 4 + 2] = conn_[3];
878 new_tets_conn[4 * 4 + 3] = conn_[2];
879 }
880 return 5;
881 } else if (free_edge_swappped0 && free_edge_swappped5) {
882 // TET0
883 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
884 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
885 if (is_rotated) {
886 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[2];
887 new_tets_conn[0 * 4 + 3] = conn_[1];
888 } else {
889 new_tets_conn[0 * 4 + 2] = conn_[1];
890 new_tets_conn[0 * 4 + 3] = edge_new_nodes_[2];
891 }
892 // TET1
893 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
894 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[2];
895 if (is_rotated) {
896 new_tets_conn[1 * 4 + 2] = conn_[2];
897 new_tets_conn[1 * 4 + 3] = conn_[3];
898 } else {
899 new_tets_conn[1 * 4 + 2] = conn_[3];
900 new_tets_conn[1 * 4 + 3] = conn_[2];
901 }
902 // TET2
903 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[4];
904 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
905 if (is_rotated) {
906 new_tets_conn[2 * 4 + 2] = conn_[0];
907 new_tets_conn[2 * 4 + 3] = conn_[1];
908 } else {
909 new_tets_conn[2 * 4 + 2] = conn_[1];
910 new_tets_conn[2 * 4 + 3] = conn_[0];
911 }
912 // TET3
913 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
914 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
915 if (is_rotated) {
916 new_tets_conn[3 * 4 + 2] = conn_[3];
917 new_tets_conn[3 * 4 + 3] = conn_[0];
918 } else {
919 new_tets_conn[3 * 4 + 2] = conn_[0];
920 new_tets_conn[3 * 4 + 3] = conn_[3];
921 }
922 // TET4
923 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[4];
924 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[2];
925 if (is_rotated) {
926 new_tets_conn[4 * 4 + 2] = edge_new_nodes_[1];
927 new_tets_conn[4 * 4 + 3] = conn_[3];
928 } else {
929 new_tets_conn[4 * 4 + 2] = conn_[3];
930 new_tets_conn[4 * 4 + 3] = edge_new_nodes_[1];
931 }
932 return 6;
933 }
934 // TET0
935 new_tets_conn[0 * 4 + 0] = edge_new_nodes_[4];
936 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[1];
937 if (is_rotated) {
938 new_tets_conn[0 * 4 + 2] = conn_[0];
939 new_tets_conn[0 * 4 + 3] = conn_[1];
940 } else {
941 new_tets_conn[0 * 4 + 2] = conn_[1];
942 new_tets_conn[0 * 4 + 3] = conn_[0];
943 }
944 // TET1
945 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
946 new_tets_conn[1 * 4 + 1] = edge_new_nodes_[2];
947 if (is_rotated) {
948 new_tets_conn[1 * 4 + 2] = conn_[2];
949 new_tets_conn[1 * 4 + 3] = edge_new_nodes_[4];
950 } else {
951 new_tets_conn[1 * 4 + 2] = edge_new_nodes_[4];
952 new_tets_conn[1 * 4 + 3] = conn_[2];
953 }
954 // TET2
955 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[4];
956 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[2];
957 if (is_rotated) {
958 new_tets_conn[2 * 4 + 2] = conn_[2];
959 new_tets_conn[2 * 4 + 3] = conn_[3];
960 } else {
961 new_tets_conn[2 * 4 + 2] = conn_[3];
962 new_tets_conn[2 * 4 + 3] = conn_[2];
963 }
964 // TET3
965 new_tets_conn[3 * 4 + 0] = edge_new_nodes_[4];
966 new_tets_conn[3 * 4 + 1] = edge_new_nodes_[2];
967 if (is_rotated) {
968 new_tets_conn[3 * 4 + 2] = conn_[0];
969 new_tets_conn[3 * 4 + 3] = edge_new_nodes_[1];
970 } else {
971 new_tets_conn[3 * 4 + 2] = edge_new_nodes_[1];
972 new_tets_conn[3 * 4 + 3] = conn_[0];
973 }
974 // TET4
975 new_tets_conn[4 * 4 + 0] = edge_new_nodes_[4];
976 new_tets_conn[4 * 4 + 1] = edge_new_nodes_[2];
977 if (is_rotated) {
978 new_tets_conn[4 * 4 + 2] = conn_[3];
979 new_tets_conn[4 * 4 + 3] = conn_[0];
980 } else {
981 new_tets_conn[4 * 4 + 2] = conn_[0];
982 new_tets_conn[4 * 4 + 3] = conn_[3];
983 }
984 return 7;
985 }
986 return type;
987}
988int tet_type_2(const EntityHandle *conn, const int *split_edges,
989 const EntityHandle *edge_new_nodes,
990 EntityHandle *new_tets_conn) {
991 if (split_edges[0] == oposite_edge[split_edges[1]]) {
992 // type 2b
993 EntityHandle n0 = conn[edges_conn[2 * split_edges[0] + 0]];
994 EntityHandle n1 = conn[edges_conn[2 * split_edges[0] + 1]];
995 EntityHandle n2 = conn[edges_conn[2 * split_edges[1] + 0]];
996 EntityHandle n3 = conn[edges_conn[2 * split_edges[1] + 1]];
997 // TET0
998 new_tets_conn[0 * 4 + 0] = edge_new_nodes[split_edges[0]];
999 new_tets_conn[0 * 4 + 1] = edge_new_nodes[split_edges[1]];
1000 new_tets_conn[0 * 4 + 2] = n2;
1001 new_tets_conn[0 * 4 + 3] = n0;
1002 // TET1
1003 new_tets_conn[1 * 4 + 0] = edge_new_nodes[split_edges[0]];
1004 new_tets_conn[1 * 4 + 1] = edge_new_nodes[split_edges[1]];
1005 new_tets_conn[1 * 4 + 2] = n1;
1006 new_tets_conn[1 * 4 + 3] = n2;
1007 // TET3
1008 new_tets_conn[2 * 4 + 0] = edge_new_nodes[split_edges[0]];
1009 new_tets_conn[2 * 4 + 1] = edge_new_nodes[split_edges[1]];
1010 new_tets_conn[2 * 4 + 2] = n0;
1011 new_tets_conn[2 * 4 + 3] = n3;
1012 // TET4
1013 new_tets_conn[3 * 4 + 0] = edge_new_nodes[split_edges[0]];
1014 new_tets_conn[3 * 4 + 1] = edge_new_nodes[split_edges[1]];
1015 new_tets_conn[3 * 4 + 2] = n3;
1016 new_tets_conn[3 * 4 + 3] = n1;
1017 return 1;
1018 } else {
1019 int sub_type = 0;
1020 // type 2a
1021 const char mach_pattern =
1022 edge_bits_mark[split_edges[0]] | edge_bits_mark[split_edges[1]];
1023 EntityHandle conn_[4] = {no_handle, no_handle, no_handle, no_handle};
1024 int edges_[6];
1025 for (int ee = 0; ee < 6; ee++) {
1026 char pattern;
1027 pattern = edge_bits_mark[edge_permutations[ee][1]] |
1028 edge_bits_mark[edge_permutations[ee][4]];
1029 if (pattern == mach_pattern) {
1030 int free_edge = oposite_edge[ee];
1031 conn_[0] = conn[edges_conn[ee * 2 + 0]];
1032 conn_[1] = conn[edges_conn[ee * 2 + 1]];
1033 conn_[2] = conn[edges_conn[free_edge * 2 + 0]];
1034 conn_[3] = conn[edges_conn[free_edge * 2 + 1]];
1035 for (int EE = 0; EE < 6; EE++)
1036 edges_[EE] = edge_permutations[ee][EE];
1037 sub_type |= 4;
1038 break;
1039 }
1040 pattern = edge_bits_mark[edge_permutations[ee][edge_mirror_cross[1]]] |
1041 edge_bits_mark[edge_permutations[ee][edge_mirror_cross[4]]];
1042 if (pattern == mach_pattern) {
1043 int free_edge = oposite_edge[ee];
1044 conn_[0] = conn[edges_conn[ee * 2 + 1]];
1045 conn_[1] = conn[edges_conn[ee * 2 + 0]];
1046 conn_[2] = conn[edges_conn[free_edge * 2 + 1]];
1047 conn_[3] = conn[edges_conn[free_edge * 2 + 0]];
1048 for (int EE = 0; EE < 6; EE++)
1049 edges_[EE] = edge_permutations[ee][edge_mirror_cross[EE]];
1050 sub_type |= 8;
1051 break;
1052 }
1053 }
1054 EntityHandle edge_new_nodes_[6];
1055 for (int ee = 0; ee < 6; ee++)
1056 edge_new_nodes_[ee] = edge_new_nodes[edges_[ee]];
1057 assert(edge_new_nodes_[1] != no_handle);
1058 assert(edge_new_nodes_[4] != no_handle);
1059 // TET0
1060 new_tets_conn[0 * 4 + 0] = conn_[1];
1061 new_tets_conn[0 * 4 + 1] = edge_new_nodes_[4];
1062 new_tets_conn[0 * 4 + 2] = edge_new_nodes_[1];
1063 new_tets_conn[0 * 4 + 3] = conn_[0];
1064 bool free_edge_swappped5 = false;
1065 if (conn_[3] < conn_[2]) {
1066 free_edge_swappped5 = true;
1067 sub_type |= 16;
1068 }
1069 if (free_edge_swappped5) {
1070 // TET1
1071 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[1];
1072 new_tets_conn[1 * 4 + 1] = conn_[2];
1073 new_tets_conn[1 * 4 + 2] = conn_[0];
1074 new_tets_conn[1 * 4 + 3] = conn_[3];
1075 // TET2
1076 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[1];
1077 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[4];
1078 new_tets_conn[2 * 4 + 2] = conn_[3];
1079 new_tets_conn[2 * 4 + 3] = conn_[0];
1080 return sub_type;
1081 }
1082 // TET1
1083 new_tets_conn[1 * 4 + 0] = edge_new_nodes_[4];
1084 new_tets_conn[1 * 4 + 1] = conn_[2];
1085 new_tets_conn[1 * 4 + 2] = conn_[0];
1086 new_tets_conn[1 * 4 + 3] = conn_[3];
1087 // TET2
1088 new_tets_conn[2 * 4 + 0] = edge_new_nodes_[1];
1089 new_tets_conn[2 * 4 + 1] = edge_new_nodes_[4];
1090 new_tets_conn[2 * 4 + 2] = conn_[2];
1091 new_tets_conn[2 * 4 + 3] = conn_[0];
1092 return sub_type;
1093 }
1094 return -1;
1095}
1096
1097void hex_to_tet_type_1(const EntityHandle *conn, EntityHandle *new_tets_conn) {
1098/*
1099 v7 ------ v6
1100 |\ |\
1101 | \ | \
1102 | v4 ------ v5
1103 | | | |
1104 v3 -|----- v2 |
1105 \ | \ |
1106 v0 ------ v1
1107*/
1108
1109// TET0
1110new_tets_conn[0 * 4 + 0] = conn[0];
1111new_tets_conn[0 * 4 + 1] = conn[5];
1112new_tets_conn[0 * 4 + 2] = conn[6];
1113new_tets_conn[0 * 4 + 3] = conn[7];
1114
1115// TET1
1116new_tets_conn[1 * 4 + 0] = conn[0];
1117new_tets_conn[1 * 4 + 1] = conn[5];
1118new_tets_conn[1 * 4 + 2] = conn[7];
1119new_tets_conn[1 * 4 + 3] = conn[4];
1120
1121// TET2
1122new_tets_conn[2 * 4 + 0] = conn[0];
1123new_tets_conn[2 * 4 + 1] = conn[3];
1124new_tets_conn[2 * 4 + 2] = conn[7];
1125new_tets_conn[2 * 4 + 3] = conn[6];
1126
1127// TET3
1128new_tets_conn[3 * 4 + 0] = conn[0];
1129new_tets_conn[3 * 4 + 1] = conn[3];
1130new_tets_conn[3 * 4 + 2] = conn[6];
1131new_tets_conn[3 * 4 + 3] = conn[2];
1132
1133// TET4
1134new_tets_conn[4 * 4 + 0] = conn[0];
1135new_tets_conn[4 * 4 + 1] = conn[1];
1136new_tets_conn[4 * 4 + 2] = conn[2];
1137new_tets_conn[4 * 4 + 3] = conn[6];
1138
1139// TET5
1140new_tets_conn[5 * 4 + 0] = conn[0];
1141new_tets_conn[5 * 4 + 1] = conn[1];
1142new_tets_conn[5 * 4 + 2] = conn[6];
1143new_tets_conn[5 * 4 + 3] = conn[5];
1144}
1145
1146void hex_to_tet_type_2(const EntityHandle *conn, EntityHandle *new_tets_conn) {
1147
1148// TET0
1149new_tets_conn[0 * 4 + 0] = conn[1];
1150new_tets_conn[0 * 4 + 1] = conn[6];
1151new_tets_conn[0 * 4 + 2] = conn[7];
1152new_tets_conn[0 * 4 + 3] = conn[4];
1153
1154// TET1
1155new_tets_conn[1 * 4 + 0] = conn[1];
1156new_tets_conn[1 * 4 + 1] = conn[6];
1157new_tets_conn[1 * 4 + 2] = conn[4];
1158new_tets_conn[1 * 4 + 3] = conn[5];
1159
1160// TET2
1161new_tets_conn[2 * 4 + 0] = conn[1];
1162new_tets_conn[2 * 4 + 1] = conn[0];
1163new_tets_conn[2 * 4 + 2] = conn[4];
1164new_tets_conn[2 * 4 + 3] = conn[7];
1165
1166// TET3
1167new_tets_conn[3 * 4 + 0] = conn[1];
1168new_tets_conn[3 * 4 + 1] = conn[0];
1169new_tets_conn[3 * 4 + 2] = conn[7];
1170new_tets_conn[3 * 4 + 3] = conn[3];
1171
1172// TET4
1173new_tets_conn[4 * 4 + 0] = conn[1];
1174new_tets_conn[4 * 4 + 1] = conn[2];
1175new_tets_conn[4 * 4 + 2] = conn[3];
1176new_tets_conn[4 * 4 + 3] = conn[7];
1177
1178// TET5
1179new_tets_conn[5 * 4 + 0] = conn[1];
1180new_tets_conn[5 * 4 + 1] = conn[2];
1181new_tets_conn[5 * 4 + 2] = conn[7];
1182new_tets_conn[5 * 4 + 3] = conn[6];
1183}
1184
1185
1186MoFEMErrorCode hex_to_tet_type_3(const EntityHandle *conn,
1187 const EntityHandle *edge_new_nodes,
1188 EntityHandle *new_tets_conn) {
1189
1190 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1191
1192 const EntityHandle FB = edge_new_nodes[12]; // bottom (v0,v1,v2,v3)
1193 const EntityHandle FT = edge_new_nodes[13]; // top (v4,v5,v6,v7)
1194 const EntityHandle FF = edge_new_nodes[14]; // front (v0,v1,v5,v4)
1195 const EntityHandle FR = edge_new_nodes[15]; // right (v1,v2,v6,v5)
1196 const EntityHandle BK = edge_new_nodes[16]; // back (v2,v3,v7,v6)
1197 const EntityHandle FL = edge_new_nodes[17]; // left (v3,v0,v4,v7)
1198
1199 // Sanity
1200 assert(FB!=no_handle && FT!=no_handle && FF!=no_handle &&
1201 FR!=no_handle && BK!=no_handle && FL!=no_handle);
1202
1203 // Vertex aliases
1204 const EntityHandle v0=conn[0], v1=conn[1], v2=conn[2], v3=conn[3];
1205 const EntityHandle v4=conn[4], v5=conn[5], v6=conn[6], v7=conn[7];
1206
1207 auto put = [&](int t, EntityHandle a, EntityHandle b, EntityHandle c, EntityHandle d){
1208 new_tets_conn[4*t+0]=a; new_tets_conn[4*t+1]=b;
1209 new_tets_conn[4*t+2]=c; new_tets_conn[4*t+3]=d;
1210 };
1211
1212 int t = 0;
1213
1214 // --- 24 face-pyramid tets (apex = face centre), 4 per face ---
1215 // Bottom face (v0,v1,v2,v3) with FB
1216 put(t++, FB, v0, v1, FF); // edge v0-v1 adjacent to front
1217 put(t++, FB, v1, v2, FR); // edge v1-v2 adjacent to right
1218 put(t++, FB, v2, v3, BK); // edge v2-v3 adjacent to back
1219 put(t++, FB, v3, v0, FL); // edge v3-v0 adjacent to left
1220
1221 // Top face (v4,v5,v6,v7) with FT
1222 put(t++, FT, v4, v5, FF); // edge v4-v5 adjacent to front
1223 put(t++, FT, v5, v6, FR);
1224 put(t++, FT, v6, v7, BK);
1225 put(t++, FT, v7, v4, FL);
1226
1227 // Front face (v0,v1,v5,v4) with FF
1228 put(t++, FF, v0, v1, FB); // edge v0-v1 adjacent to bottom
1229 put(t++, FF, v1, v5, FR); // edge v1-v5 adjacent to right
1230 put(t++, FF, v5, v4, FT); // edge v5-v4 adjacent to top
1231 put(t++, FF, v4, v0, FL); // edge v4-v0 adjacent to left
1232
1233 // Right face (v1,v2,v6,v5) with FR
1234 put(t++, FR, v1, v2, FB);
1235 put(t++, FR, v2, v6, BK);
1236 put(t++, FR, v6, v5, FT);
1237 put(t++, FR, v5, v1, FF);
1238
1239 // Back face (v2,v3,v7,v6) with BK
1240 put(t++, BK, v2, v3, FB);
1241 put(t++, BK, v3, v7, FL);
1242 put(t++, BK, v7, v6, FT);
1243 put(t++, BK, v6, v2, FR);
1244
1245 // Left face (v3,v0,v4,v7) with FL
1246 put(t++, FL, v3, v0, FB);
1247 put(t++, FL, v0, v4, FF);
1248 put(t++, FL, v4, v7, FT);
1249 put(t++, FL, v7, v3, BK);
1250
1251 // --- 4 central tets: split the octahedron of {FB,FT,FF,FR,BK,FL} ---
1252 // One consistent 4-tet split (any valid octahedron split into 4 tets works):
1253 put(t++, FB, FF, FR, FT);
1254 put(t++, FB, FR, BK, FT);
1255 put(t++, FB, BK, FL, FT);
1256 put(t++, FB, FL, FF, FT);
1257
1258 assert(t == 28);
1259 MoFEMFunctionReturnHot(0);
1260}
1261
1262void tet_type_1(const EntityHandle *conn, const int split_edge,
1263 const EntityHandle edge_new_node, EntityHandle *new_tets_conn) {
1264 // TET0
1265 new_tets_conn[0 * 4 + 0] = edge_new_node;
1266 new_tets_conn[0 * 4 + 1] = conn[edges_conn[2 * split_edge + 0]];
1267 new_tets_conn[0 * 4 + 2] = conn[edges_conn[2 * oposite_edge[split_edge] + 1]];
1268 new_tets_conn[0 * 4 + 3] = conn[edges_conn[2 * oposite_edge[split_edge] + 0]];
1269 // TET1
1270 new_tets_conn[1 * 4 + 0] = edge_new_node;
1271 new_tets_conn[1 * 4 + 1] = conn[edges_conn[2 * split_edge + 1]];
1272 new_tets_conn[1 * 4 + 2] = conn[edges_conn[2 * oposite_edge[split_edge] + 0]];
1273 new_tets_conn[1 * 4 + 3] = conn[edges_conn[2 * oposite_edge[split_edge] + 1]];
1274}
1275
1276// TRIS
1277MoFEMErrorCode tri_type_3(const EntityHandle *conn,
1278 const EntityHandle *edge_new_nodes,
1279 EntityHandle *new_tris_conn) {
1280 MoFEMFunctionBeginHot;
1281 // TRI0
1282 new_tris_conn[0 * 3 + 0] = conn[0];
1283 new_tris_conn[0 * 3 + 1] = edge_new_nodes[0];
1284 new_tris_conn[0 * 3 + 2] = edge_new_nodes[2];
1285
1286 // TRI1
1287 new_tris_conn[1 * 3 + 0] = edge_new_nodes[0];
1288 new_tris_conn[1 * 3 + 1] = conn[1];
1289 new_tris_conn[1 * 3 + 2] = edge_new_nodes[1];
1290
1291 // TRI2
1292 new_tris_conn[2 * 3 + 0] = edge_new_nodes[2];
1293 new_tris_conn[2 * 3 + 1] = edge_new_nodes[1];
1294 new_tris_conn[2 * 3 + 2] = conn[2];
1295
1296 // TRI3
1297 new_tris_conn[3 * 3 + 0] = edge_new_nodes[0];
1298 new_tris_conn[3 * 3 + 1] = edge_new_nodes[1];
1299 new_tris_conn[3 * 3 + 2] = edge_new_nodes[2];
1300 MoFEMFunctionReturnHot(0);
1301}
1302
1303MoFEMErrorCode tri_type_1(const EntityHandle *conn, const int split_edge,
1304 const EntityHandle edge_new_node,
1305 EntityHandle *new_tris_conn) {
1306
1307 MoFEMFunctionBeginHot;
1308
1309 if (split_edge == 0) {
1310 // TRI0
1311 new_tris_conn[0 * 3 + 0] = conn[0];
1312 new_tris_conn[0 * 3 + 1] = edge_new_node;
1313 new_tris_conn[0 * 3 + 2] = conn[2];
1314 // TRI1
1315 new_tris_conn[1 * 3 + 0] = edge_new_node;
1316 new_tris_conn[1 * 3 + 1] = conn[1];
1317 new_tris_conn[1 * 3 + 2] = conn[2];
1318 } else if (split_edge == 1) {
1319 // TRI0
1320 new_tris_conn[0 * 3 + 0] = conn[0];
1321 new_tris_conn[0 * 3 + 1] = conn[1];
1322 new_tris_conn[0 * 3 + 2] = edge_new_node;
1323 // TRI1
1324 new_tris_conn[1 * 3 + 0] = conn[0];
1325 new_tris_conn[1 * 3 + 1] = edge_new_node;
1326 new_tris_conn[1 * 3 + 2] = conn[2];
1327 } else if (split_edge == 2) {
1328 // TRI0
1329 new_tris_conn[0 * 3 + 0] = conn[0];
1330 new_tris_conn[0 * 3 + 1] = conn[1];
1331 new_tris_conn[0 * 3 + 2] = edge_new_node;
1332 // TRI1
1333 new_tris_conn[1 * 3 + 0] = edge_new_node;;
1334 new_tris_conn[1 * 3 + 1] = conn[1];
1335 new_tris_conn[1 * 3 + 2] = conn[2];
1336 } else {
1337 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Impossible case");
1338 }
1339
1340 MoFEMFunctionReturnHot(0);
1341}
1342
1343MoFEMErrorCode tri_type_2(const EntityHandle *conn, const int *split_edges,
1344 const EntityHandle *edge_new_nodes,
1345 EntityHandle *new_tris_conn) {
1346
1347 MoFEMFunctionBeginHot;
1348
1349 if(split_edges[0] == 0 && split_edges[1] == 1) {
1350 // TRI0
1351 new_tris_conn[0 * 3 + 0] = conn[0];
1352 new_tris_conn[0 * 3 + 1] = edge_new_nodes[0];
1353 new_tris_conn[0 * 3 + 2] = edge_new_nodes[1];
1354 // TRI1
1355 new_tris_conn[1 * 3 + 0] = edge_new_nodes[0];
1356 new_tris_conn[1 * 3 + 1] = conn[1];
1357 new_tris_conn[1 * 3 + 2] = edge_new_nodes[1];
1358 // TRI2
1359 new_tris_conn[2 * 3 + 0] = conn[0];
1360 new_tris_conn[2 * 3 + 1] = edge_new_nodes[1];
1361 new_tris_conn[2 * 3 + 2] = conn[2];
1362 } else if(split_edges[0] == 0 && split_edges[1] == 2) {
1363 // TRI0
1364 new_tris_conn[0 * 3 + 0] = conn[0];
1365 new_tris_conn[0 * 3 + 1] = edge_new_nodes[0];
1366 new_tris_conn[0 * 3 + 2] = edge_new_nodes[2];
1367 // TRI1
1368 new_tris_conn[1 * 3 + 0] = edge_new_nodes[0];
1369 new_tris_conn[1 * 3 + 1] = conn[1];
1370 new_tris_conn[1 * 3 + 2] = edge_new_nodes[2];
1371 // TRI2
1372 new_tris_conn[2 * 3 + 0] = edge_new_nodes[2];
1373 new_tris_conn[2 * 3 + 1] = conn[1];
1374 new_tris_conn[2 * 3 + 2] = conn[2];
1375 } else if(split_edges[0] == 1 && split_edges[1] == 2) {
1376 // TRI0
1377 new_tris_conn[0 * 3 + 0] = conn[0];
1378 new_tris_conn[0 * 3 + 1] = conn[1];
1379 new_tris_conn[0 * 3 + 2] = edge_new_nodes[2];
1380 // TRI1
1381 new_tris_conn[1 * 3 + 0] = edge_new_nodes[2];
1382 new_tris_conn[1 * 3 + 1] = conn[1];
1383 new_tris_conn[1 * 3 + 2] = edge_new_nodes[1];
1384 // TRI2
1385 new_tris_conn[2 * 3 + 0] = edge_new_nodes[2];
1386 new_tris_conn[2 * 3 + 1] = edge_new_nodes[1];
1387 new_tris_conn[2 * 3 + 2] = conn[2];
1388 } else {
1389 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Impossible case");
1390 }
1391
1392 MoFEMFunctionReturnHot(0);
1393}
1394
1395// PRISM
1396MoFEMErrorCode prism_type_1(const EntityHandle *conn,
1397 const BitRefEdges split_edges,
1398 const EntityHandle *edge_new_nodes,
1399 EntityHandle *new_prism_conn) {
1400 MoFEMFunctionBeginHot;
1401 int ee = 0;
1402 for (; ee < 6; ee++) {
1403 if (split_edges.test(ee))
1404 break;
1405 }
1406 if (ee > 2)
1407 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1408 const int cycle_edges[3][6] = {
1409 {0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}};
1410 const int cycle_nodes[3][6] = {
1411 {0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}};
1412 if (edge_new_nodes[cycle_nodes[ee][0]] == no_handle)
1413 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1414 if (edge_new_nodes[cycle_nodes[ee][3]] == no_handle)
1415 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1416 // PRISM0
1417 new_prism_conn[0 * 6 + 0] = edge_new_nodes[cycle_edges[ee][0]];
1418 new_prism_conn[0 * 6 + 1] = conn[cycle_nodes[ee][2]];
1419 new_prism_conn[0 * 6 + 2] = conn[cycle_nodes[ee][0]];
1420 new_prism_conn[0 * 6 + 3] = edge_new_nodes[cycle_edges[ee][3]];
1421 new_prism_conn[0 * 6 + 4] = conn[cycle_nodes[ee][5]];
1422 new_prism_conn[0 * 6 + 5] = conn[cycle_nodes[ee][3]];
1423 // PRISM1
1424 new_prism_conn[1 * 6 + 0] = edge_new_nodes[cycle_edges[ee][0]];
1425 new_prism_conn[1 * 6 + 1] = conn[cycle_nodes[ee][2]];
1426 new_prism_conn[1 * 6 + 2] = conn[cycle_nodes[ee][1]];
1427 new_prism_conn[1 * 6 + 3] = edge_new_nodes[cycle_edges[ee][3]];
1428 new_prism_conn[1 * 6 + 4] = conn[cycle_nodes[ee][5]];
1429 new_prism_conn[1 * 6 + 5] = conn[cycle_nodes[ee][4]];
1430 MoFEMFunctionReturnHot(0);
1431}
1432MoFEMErrorCode prism_type_2(const EntityHandle *conn,
1433 const BitRefEdges split_edges,
1434 const EntityHandle *edge_new_nodes,
1435 EntityHandle *new_prism_conn) {
1436 MoFEMFunctionBeginHot;
1437 const int cycle_edges[3][6] = {
1438 {0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}};
1439 const int cycle_nodes[3][6] = {
1440 {0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}};
1441 int ee = 0;
1442 for (; ee < 3; ee++) {
1443 BitRefEdges mach_pattern(0);
1444 mach_pattern.set(cycle_edges[ee][0]);
1445 mach_pattern.set(cycle_edges[ee][2]);
1446 mach_pattern.set(cycle_edges[ee][3]);
1447 mach_pattern.set(cycle_edges[ee][5]);
1448 if (mach_pattern == split_edges)
1449 break;
1450 }
1451 if (ee > 2)
1452 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1453 EntityHandle _conn_[6], _edge_new_nodes_[6];
1454 int nn = 0;
1455 for (; nn < 6; nn++) {
1456 _conn_[nn] = conn[cycle_nodes[ee][nn]];
1457 _edge_new_nodes_[nn] = edge_new_nodes[cycle_edges[ee][nn]];
1458 }
1459 if (_conn_[1] < _conn_[2]) {
1460 EntityHandle _conn____;
1461 _conn____ = _conn_[2];
1462 _conn_[2] = _conn_[1];
1463 _conn_[1] = _conn____;
1464 _conn____ = _conn_[5];
1465 _conn_[5] = _conn_[4];
1466 _conn_[4] = _conn____;
1467 _conn____ = _edge_new_nodes_[0];
1468 _edge_new_nodes_[0] = _edge_new_nodes_[2];
1469 _edge_new_nodes_[2] = _conn____;
1470 _conn____ = _edge_new_nodes_[6 - 3];
1471 _edge_new_nodes_[6 - 3] = _edge_new_nodes_[8 - 3];
1472 _edge_new_nodes_[8 - 3] = _conn____;
1473 }
1474 if (_edge_new_nodes_[0] == no_handle)
1475 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1476 if (_edge_new_nodes_[2] == no_handle)
1477 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1478 if (_edge_new_nodes_[6 - 3] == no_handle)
1479 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1480 if (_edge_new_nodes_[8 - 3] == no_handle)
1481 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1482 // PRIMS0
1483 new_prism_conn[0 * 6 + 0] = _conn_[0];
1484 new_prism_conn[0 * 6 + 1] = _edge_new_nodes_[0];
1485 new_prism_conn[0 * 6 + 2] = _edge_new_nodes_[2];
1486 new_prism_conn[0 * 6 + 3] = _conn_[3];
1487 new_prism_conn[0 * 6 + 4] = _edge_new_nodes_[6 - 3];
1488 new_prism_conn[0 * 6 + 5] = _edge_new_nodes_[8 - 3];
1489 // PRISM1
1490 new_prism_conn[1 * 6 + 0] = _conn_[1];
1491 new_prism_conn[1 * 6 + 1] = _conn_[2];
1492 new_prism_conn[1 * 6 + 2] = _edge_new_nodes_[0];
1493 new_prism_conn[1 * 6 + 3] = _conn_[4];
1494 new_prism_conn[1 * 6 + 4] = _conn_[5];
1495 new_prism_conn[1 * 6 + 5] = _edge_new_nodes_[6 - 3];
1496 // PRISM2
1497 new_prism_conn[2 * 6 + 0] = _conn_[2];
1498 new_prism_conn[2 * 6 + 1] = _edge_new_nodes_[0];
1499 new_prism_conn[2 * 6 + 2] = _edge_new_nodes_[2];
1500 new_prism_conn[2 * 6 + 3] = _conn_[5];
1501 new_prism_conn[2 * 6 + 4] = _edge_new_nodes_[6 - 3];
1502 new_prism_conn[2 * 6 + 5] = _edge_new_nodes_[8 - 3];
1503 MoFEMFunctionReturnHot(0);
1504}
1505MoFEMErrorCode prism_type_3(const EntityHandle *conn,
1506 const BitRefEdges split_edges,
1507 const EntityHandle *edge_new_nodes,
1508 EntityHandle *new_prism_conn) {
1509 MoFEMFunctionBeginHot;
1510 int ee = 0;
1511 for (; ee < 6; ee++) {
1512 if (!split_edges.test(ee))
1513 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1514 }
1515 // PRISM0
1516 new_prism_conn[0 * 6 + 0] = edge_new_nodes[0];
1517 new_prism_conn[0 * 6 + 1] = edge_new_nodes[2];
1518 new_prism_conn[0 * 6 + 2] = conn[0];
1519 new_prism_conn[0 * 6 + 3] = edge_new_nodes[6 - 3];
1520 new_prism_conn[0 * 6 + 4] = edge_new_nodes[8 - 3];
1521 new_prism_conn[0 * 6 + 5] = conn[3];
1522 // PRISM1
1523 new_prism_conn[1 * 6 + 0] = edge_new_nodes[0];
1524 new_prism_conn[1 * 6 + 1] = edge_new_nodes[1];
1525 new_prism_conn[1 * 6 + 2] = conn[1];
1526 new_prism_conn[1 * 6 + 3] = edge_new_nodes[6 - 3];
1527 new_prism_conn[1 * 6 + 4] = edge_new_nodes[7 - 3];
1528 new_prism_conn[1 * 6 + 5] = conn[4];
1529 // PRISM2
1530 new_prism_conn[2 * 6 + 0] = edge_new_nodes[1];
1531 new_prism_conn[2 * 6 + 1] = edge_new_nodes[2];
1532 new_prism_conn[2 * 6 + 2] = conn[2];
1533 new_prism_conn[2 * 6 + 3] = edge_new_nodes[7 - 3];
1534 new_prism_conn[2 * 6 + 4] = edge_new_nodes[8 - 3];
1535 new_prism_conn[2 * 6 + 5] = conn[5];
1536 // PRISM3
1537 new_prism_conn[3 * 6 + 0] = edge_new_nodes[0];
1538 new_prism_conn[3 * 6 + 1] = edge_new_nodes[1];
1539 new_prism_conn[3 * 6 + 2] = edge_new_nodes[2];
1540 new_prism_conn[3 * 6 + 3] = edge_new_nodes[6 - 3];
1541 new_prism_conn[3 * 6 + 4] = edge_new_nodes[7 - 3];
1542 new_prism_conn[3 * 6 + 5] = edge_new_nodes[8 - 3];
1543 MoFEMFunctionReturnHot(0);
1544}
1545
1546MoFEMErrorCode quad_split_all_edges(const EntityHandle *conn,
1547 const BitRefEdges split_edges,
1548 const EntityHandle *edge_new_nodes,
1549 EntityHandle *new_quad_conn) {
1550
1551 return 0;
1552 }
1553
1554} // namespace MoFEM
a
constexpr double a
Definition approx_sphere.cpp:30
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
FR
FR
Definition free_surface.cpp:624
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:39
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::BitRefEdges
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition Types.hpp:34
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::tet_type_3
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:565
MoFEM::tet_type_6
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:31
MoFEM::hex_to_tet_type_3
MoFEMErrorCode hex_to_tet_type_3(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:1186
MoFEM::tri_type_2
MoFEMErrorCode tri_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1343
MoFEM::quad_split_all_edges
MoFEMErrorCode quad_split_all_edges(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_quad_conn)
MoFEM::hex_to_tet_type_2
void hex_to_tet_type_2(const EntityHandle *conn, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:1146
MoFEM::edge_mirror_cross
static constexpr int edge_mirror_cross[6]
Definition EntityRefine.cpp:23
MoFEM::prism_type_2
MoFEMErrorCode prism_type_2(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1432
MoFEM::edge_bits_mark
static constexpr char edge_bits_mark[]
Definition EntityRefine.cpp:29
MoFEM::tri_type_1
MoFEMErrorCode tri_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1303
MoFEM::prism_type_3
MoFEMErrorCode prism_type_3(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1505
MoFEM::hex_to_tet_type_1
void hex_to_tet_type_1(const EntityHandle *conn, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:1097
MoFEM::oposite_edge
static constexpr int oposite_edge[]
Definition EntityRefine.cpp:19
MoFEM::tri_type_3
MoFEMErrorCode tri_type_3(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1277
MoFEM::prism_type_1
MoFEMErrorCode prism_type_1(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1396
MoFEM::cyclic_edge_rotate_face_3
static constexpr int cyclic_edge_rotate_face_3[3][6]
Definition EntityRefine.cpp:27
MoFEM::edge_mirror_vertical
static constexpr int edge_mirror_vertical[6]
Definition EntityRefine.cpp:24
MoFEM::tet_type_2
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:988
MoFEM::edge_permutations
static constexpr int edge_permutations[6][6]
Definition EntityRefine.cpp:20
MoFEM::tet_type_5
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:133
MoFEM::cyclic_node_rotate_face_3
static constexpr int cyclic_node_rotate_face_3[3][4]
Definition EntityRefine.cpp:25
MoFEM::no_handle
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition Common.hpp:12
MoFEM::tet_type_4
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:239
MoFEM::tet_type_1
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:1262
MoFEM::edges_conn
static constexpr int edges_conn[]
Definition EntityRefine.cpp:18
t
constexpr double t
plate stiffness
Definition plate.cpp:58
Generated by Doxygen 1.9.8 and hosted at University of Glasgow