LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
Loading...
Searching...
No Matches
hybrid2d_refinement_pattern.cc
1
7#include "hybrid2d_refinement_pattern.h"
8
9namespace lf::refinement {
10
11std::ostream &operator<<(std::ostream &o, const RefPat &refpat) {
12 switch (refpat) {
13 case rp_nil:
14 o << "rp_nil";
15 break;
16 case rp_copy:
17 o << "rp_copy";
18 break;
19 case rp_split:
20 o << "rp_split";
21 break;
22 case rp_bisect:
23 o << "rp_bisect";
24 break;
25 case rp_trisect:
26 o << "rp_trisect";
27 break;
28 case rp_trisect_left:
29 o << "rp_trisect_left";
30 break;
31 case rp_quadsect:
32 o << "rp_quadsect";
33 break;
34 case rp_threeedge:
35 o << "rp_threeedge";
36 break;
37 case rp_regular:
38 o << "rp_regular";
39 break;
40 case rp_barycentric:
41 o << "rp_barycentric";
42 break;
43 }
44 return o;
45}
46
47// Implementation for RefinemeentPattern
48// NOLINTNEXTLINE(readability-function-cognitive-complexity)
49lf::base::size_type Hybrid2DRefinementPattern::NumChildren(
50 lf::base::dim_t codim) const {
51 LF_VERIFY_MSG(codim <= ref_el_.Dimension(),
52 "Illegal codim = " << codim << " for " << ref_el_.ToString());
53 // Depending on the type of cell do something different
54 switch (ref_el_) {
55 case lf::base::RefEl::kPoint(): {
56 switch (ref_pat_) {
57 case RefPat::rp_nil: {
58 return 0;
59 }
60 case RefPat::rp_copy: {
61 return 1;
62 }
63 default: {
64 LF_VERIFY_MSG(false, "Illegal refinement pattern for point");
65 }
66 } // end swtich ref_pat
67 break;
68 } // end case of a simple point
69 case lf::base::RefEl::kSegment(): {
70 switch (codim) {
71 case 0: {
72 switch (ref_pat_) {
73 case RefPat::rp_nil: {
74 return 0;
75 }
76 case RefPat::rp_copy: {
77 return 1;
78 }
79 case RefPat::rp_split: {
80 return 2;
81 }
82 default: {
83 LF_VERIFY_MSG(false, "refinement pattern "
84 << (int)ref_pat_
85 << "not implemented for edge!");
86 break;
87 }
88 } // end switch ref_pat
89 break;
90 }
91 case 1: {
92 switch (ref_pat_) {
93 case RefPat::rp_nil:
94 case RefPat::rp_copy: {
95 return 0;
96 }
97 case RefPat::rp_split: {
98 return 1;
99 }
100 default: {
101 LF_VERIFY_MSG(false, "refinement pattern "
102 << (int)ref_pat_
103 << "not implemented for edge!");
104 break;
105 }
106 } // end switch ref_pat
107 break;
108 }
109 default:
110 LF_VERIFY_MSG(false, "invalid codim");
111 } // end switch codim
112 break;
113 } // end case of an edge
114 case lf::base::RefEl::kTria(): {
115 switch (codim) {
116 case 0: {
117 switch (ref_pat_) {
118 case RefPat::rp_nil: {
119 return 0;
120 }
121 case RefPat::rp_copy: {
122 return 1;
123 }
124 case RefPat::rp_bisect: {
125 return 2;
126 }
127 case RefPat::rp_trisect:
128 case RefPat::rp_trisect_left: {
129 return 3;
130 }
131 case RefPat::rp_quadsect:
132 case RefPat::rp_regular: {
133 return 4;
134 }
135 case RefPat::rp_barycentric: {
136 return 6;
137 }
138 default: {
139 LF_VERIFY_MSG(false, "refinement pattern "
140 << (int)ref_pat_
141 << "not implemented for triangle!");
142 break;
143 }
144 } // end switch ref_pat
145 break;
146 default:
147 LF_VERIFY_MSG(false, "invalid codim");
148 } // end case codim = 0
149 case 1: {
150 switch (ref_pat_) {
151 case RefPat::rp_nil:
152 case RefPat::rp_copy: {
153 return 0;
154 }
155 case RefPat::rp_bisect: {
156 return 1;
157 }
158 case RefPat::rp_trisect:
159 case RefPat::rp_trisect_left: {
160 return 2;
161 }
162 case RefPat::rp_quadsect:
163 case RefPat::rp_regular: {
164 return 3;
165 }
166 case RefPat::rp_barycentric: {
167 return 6;
168 }
169 default: {
170 LF_VERIFY_MSG(false, "refinement pattern "
171 << ref_pat_
172 << "not implemented for triangle!");
173 break;
174 }
175 } // end switch ref_pat
176 break;
177 } // end case codim = 1
178 case 2: {
179 if (ref_pat_ == RefPat::rp_barycentric) {
180 return 1;
181 }
182 return 0;
183 }
184 } // end switch codim
185 break;
186 } // end case of a triangle
187 case lf::base::RefEl::kQuad(): {
188 switch (codim) {
189 case 0: {
190 switch (ref_pat_) {
191 case RefPat::rp_nil: {
192 return 0;
193 }
194 case RefPat::rp_copy: {
195 return 1;
196 }
197 case RefPat::rp_trisect: {
198 return 3;
199 }
200 case RefPat::rp_quadsect: {
201 return 4;
202 }
203 case RefPat::rp_bisect:
204 case RefPat::rp_split: {
205 return 2;
206 }
207 case RefPat::rp_threeedge:
208 case RefPat::rp_barycentric:
209 case RefPat::rp_regular: {
210 return 4;
211 }
212 default: {
213 LF_VERIFY_MSG(false, "refinement pattern "
214 << ref_pat_
215 << "not implemented for quadrilateral!");
216 break;
217 }
218 } // end switch ref_pat
219 break;
220 }
221 case 1: {
222 switch (ref_pat_) {
223 case RefPat::rp_nil:
224 case RefPat::rp_copy: {
225 return 0;
226 }
227 case RefPat::rp_trisect: {
228 return 2;
229 }
230 case RefPat::rp_quadsect: {
231 return 3;
232 }
233 case RefPat::rp_bisect:
234 case RefPat::rp_split: {
235 return 1;
236 }
237 case RefPat::rp_threeedge: {
238 return 3;
239 }
240 case RefPat::rp_barycentric:
241 case RefPat::rp_regular: {
242 return 4;
243 }
244 default: {
245 LF_VERIFY_MSG(false, "refinement pattern "
246 << (int)ref_pat_
247 << "not implemented for quadrilateral!");
248 break;
249 }
250 } // end switch ref_pat
251 break;
252 }
253 case 2: {
254 if ((ref_pat_ == RefPat::rp_barycentric) ||
255 (ref_pat_ == RefPat::rp_regular)) {
256 return 1;
257 }
258 return 0;
259 break;
260 }
261 default:
262 LF_VERIFY_MSG(false, "invalid codim");
263 } // end switch codim
264 break;
265 } // end case of a quadrilateral
266 default:
267 LF_VERIFY_MSG(false, "unsupported RefEl");
268 } // end switch cell type
269 return 0;
270}
271
272std::vector<Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>>
273// NOLINTNEXTLINE(readability-function-cognitive-complexity)
274Hybrid2DRefinementPattern::ChildPolygons(lf::base::dim_t codim) const {
275 LF_VERIFY_MSG(codim <= ref_el_.Dimension(),
276 "Illegal codim = " << codim << " for " << ref_el_.ToString());
277 // Local variable for accumulating information about interior entities
278 // created during refinement
279 std::vector<Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>> child_poly{};
280
281 // Lattice point coordinates (lattice_const_ should be a multiple of 6)
282 const int lt_half = base::narrow<int>(lattice_const_ / 2);
283 const int lt_third = base::narrow<int>(lattice_const_ / 3);
284 const int lt_one = base::narrow<int>(lattice_const_);
285 // Depending on the type of cell do something different
286 switch (ref_el_) {
287 case lf::base::RefEl::kPoint(): {
288 if (ref_pat_ != RefPat::rp_nil) {
289 child_poly.emplace_back(0, 1);
290 }
291 break;
292 } // end case of a simple point
293 case lf::base::RefEl::kSegment(): { // case of an edge
294 switch (codim) {
295 case 0: { // child entities are segments as well
296 switch (ref_pat_) {
297 case RefPat::rp_nil: {
298 break;
299 }
300 case RefPat::rp_copy: { // copy of segment, 1 child
301 Eigen::Matrix<int, 1, 2> seg_ref_coord;
302 seg_ref_coord << 0, lt_one;
303 child_poly.emplace_back(seg_ref_coord);
304 break;
305 }
306 case RefPat::rp_split: { // splitting of segment, 2 children
307 Eigen::Matrix<int, 1, 2> part_ref_coord;
308 part_ref_coord << 0, lt_half;
309 child_poly.emplace_back(part_ref_coord);
310 part_ref_coord << lt_half, lt_one;
311 child_poly.emplace_back(part_ref_coord);
312 break;
313 }
314 default: {
315 LF_VERIFY_MSG(false, "refinement pattern "
316 << (int)ref_pat_
317 << "not implemented for edge!");
318 break;
319 }
320 } // end switch ref_pat
321 break;
322 } // case codim = 0
323 case 1: { // child entities are points
324 if (ref_pat_ == RefPat::rp_split) {
325 // child is the midpoint
326 Eigen::Matrix<int, 1, 1> point_ref_coord;
327 point_ref_coord << lt_half;
328 child_poly.emplace_back(point_ref_coord);
329 }
330 break;
331 } // end case codim = 1
332 default:
333 LF_VERIFY_MSG(false, "invalid codim");
334 } // end switch codim
335 break;
336 } // end case of an edge
337 case lf::base::RefEl::kTria(): {
338 // **********************************************************************
339 // Triangle
340 // **********************************************************************
341 // Lattice coordinates of special points in the triangle
342 // Here we use knowledge about the conventions underlying node
343 // and edge numbering in a triangle as defined in RefEl.
344 // This information could also be retrieved from the reference element.
345 Eigen::MatrixXi lt_node_coords(2, 3);
346 lt_node_coords << 0, lt_one, 0, 0, 0, lt_one;
347 Eigen::MatrixXi lt_midpoint_coords(2, 3);
348 lt_midpoint_coords << lt_half, lt_half, 0, 0, lt_half, lt_half;
349
350 // Remap local indices according to anchor values
351 const unsigned int mod_0 = (0 + anchor_) % 3;
352 const unsigned int mod_1 = (1 + anchor_) % 3;
353 const unsigned int mod_2 = (2 + anchor_) % 3;
354
355 switch (codim) {
356 case 0: {
357 // For temporary storage of triangular lattice polygon
358 Eigen::MatrixXi child_coords(2, 3);
359
360 switch (ref_pat_) {
361 case RefPat::rp_nil: {
362 break;
363 }
364 case RefPat::rp_copy: {
365 child_coords = lt_node_coords;
366 child_poly.push_back(child_coords);
367 break;
368 }
369 case RefPat::rp_bisect: {
370 LF_VERIFY_MSG(
371 anchor_set_,
372 "Anchor must be set for bisection refinement of triangle");
373 // Splitting a triangle in two by bisecting anchor edge
374 child_coords.col(0) = lt_node_coords.col(mod_0);
375 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
376 child_coords.col(2) = lt_node_coords.col(mod_2);
377 child_poly.push_back(child_coords);
378
379 child_coords.col(0) = lt_node_coords.col(mod_1);
380 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
381 child_coords.col(2) = lt_node_coords.col(mod_2);
382 child_poly.push_back(child_coords);
383 break;
384 }
385 case RefPat::rp_trisect: {
386 LF_VERIFY_MSG(
387 anchor_set_,
388 "Anchor must be set for trisection refinement of triangle");
389 // Bisect through anchor edge first and then bisect through
390 // edge with the next larger index (mod 3); creates three
391 // child triangles.
392 child_coords.col(0) = lt_node_coords.col(mod_0);
393 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
394 child_coords.col(2) = lt_node_coords.col(mod_2);
395 child_poly.push_back(child_coords);
396
397 child_coords.col(0) = lt_node_coords.col(mod_1);
398 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
399 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
400 child_poly.push_back(child_coords);
401
402 child_coords.col(0) = lt_node_coords.col(mod_2);
403 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
404 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
405 child_poly.push_back(child_coords);
406 break;
407 }
408 case RefPat::rp_trisect_left: {
409 LF_VERIFY_MSG(
410 anchor_set_,
411 "Anchor must be set for trisection refinement of triangle");
412
413 // Bisect through anchor edge first and then bisect through
414 // edge with the next smaller index (mod 3); creates three
415 // child triangles.
416 child_coords.col(0) = lt_node_coords.col(mod_0);
417 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
418 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
419 child_poly.push_back(child_coords);
420
421 child_coords.col(0) = lt_node_coords.col(mod_1);
422 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
423 child_coords.col(2) = lt_node_coords.col(mod_2);
424 child_poly.push_back(child_coords);
425
426 child_coords.col(0) = lt_node_coords.col(mod_2);
427 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
428 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
429 child_poly.push_back(child_coords);
430 break;
431 }
432 case RefPat::rp_quadsect: {
433 LF_VERIFY_MSG(
434 anchor_set_,
435 "Anchor must be set for quadsection refinement of triangle");
436 // Bisect through the anchor edge first and then
437 // through the two remaining edges; creates four child
438 // triangles; every edge is split.
439 child_coords.col(0) = lt_node_coords.col(mod_0);
440 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
441 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
442 child_poly.push_back(child_coords);
443
444 child_coords.col(0) = lt_node_coords.col(mod_1);
445 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
446 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
447 child_poly.push_back(child_coords);
448
449 child_coords.col(0) = lt_node_coords.col(mod_2);
450 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
451 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
452 child_poly.push_back(child_coords);
453
454 child_coords.col(0) = lt_node_coords.col(mod_2);
455 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
456 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
457 child_poly.push_back(child_coords);
458 break;
459 }
460 case RefPat::rp_regular: {
461 // Split triangle into four small congruent triangles
462 child_coords.col(0) = lt_node_coords.col(0);
463 child_coords.col(1) = lt_midpoint_coords.col(0);
464 child_coords.col(2) = lt_midpoint_coords.col(2);
465 child_poly.push_back(child_coords);
466
467 child_coords.col(0) = lt_node_coords.col(1);
468 child_coords.col(1) = lt_midpoint_coords.col(0);
469 child_coords.col(2) = lt_midpoint_coords.col(1);
470 child_poly.push_back(child_coords);
471
472 child_coords.col(0) = lt_node_coords.col(2);
473 child_coords.col(1) = lt_midpoint_coords.col(2);
474 child_coords.col(2) = lt_midpoint_coords.col(1);
475 child_poly.push_back(child_coords);
476
477 child_coords.col(0) = lt_midpoint_coords.col(0);
478 child_coords.col(1) = lt_midpoint_coords.col(1);
479 child_coords.col(2) = lt_midpoint_coords.col(2);
480 child_poly.push_back(child_coords);
481 break;
482 }
483 case RefPat::rp_barycentric: {
484 // Split triangle into 5 smaller triangles by connecting
485 // the center of gravity with the vertices and the midpoints
486 // of the edges.
487
488 // Obtain coordinates of center of gravity
489 const Eigen::Vector2i lt_baryc_coords =
490 Eigen::Vector2i({lt_third, lt_third});
491
492 child_coords.col(0) = lt_node_coords.col(0);
493 child_coords.col(1) = lt_midpoint_coords.col(0);
494 child_coords.col(2) = lt_baryc_coords;
495 child_poly.push_back(child_coords);
496
497 child_coords.col(0) = lt_node_coords.col(1);
498 child_coords.col(1) = lt_midpoint_coords.col(0);
499 child_coords.col(2) = lt_baryc_coords;
500 child_poly.push_back(child_coords);
501
502 child_coords.col(0) = lt_node_coords.col(1);
503 child_coords.col(1) = lt_midpoint_coords.col(1);
504 child_coords.col(2) = lt_baryc_coords;
505 child_poly.push_back(child_coords);
506
507 child_coords.col(0) = lt_node_coords.col(2);
508 child_coords.col(1) = lt_midpoint_coords.col(1);
509 child_coords.col(2) = lt_baryc_coords;
510 child_poly.push_back(child_coords);
511
512 child_coords.col(0) = lt_node_coords.col(2);
513 child_coords.col(1) = lt_midpoint_coords.col(2);
514 child_coords.col(2) = lt_baryc_coords;
515 child_poly.push_back(child_coords);
516
517 child_coords.col(0) = lt_node_coords.col(0);
518 child_coords.col(1) = lt_midpoint_coords.col(2);
519 child_coords.col(2) = lt_baryc_coords;
520 child_poly.push_back(child_coords);
521 break;
522 }
523 default: {
524 LF_VERIFY_MSG(false, "refinement pattern "
525 << (int)ref_pat_
526 << "not implemented for triangle!");
527 break;
528 }
529 } // end switch ref_pat
530 break;
531 } // end case codim = 0
532 case 1: { // codim == 1
533 // Return the location of interior edges of the triangle
534 // This location is defined by two lattice points given by the
535 // rows of the following matrix
536 Eigen::MatrixXi child_coords(2, 2);
537
538 // Different layout of interior edges for different refinement
539 // patterns
540 switch (ref_pat_) {
541 case RefPat::rp_nil:
542 case RefPat::rp_copy: {
543 // No interior edges in these cases
544 break;
545 }
546 case RefPat::rp_bisect: {
547 LF_VERIFY_MSG(
548 anchor_set_,
549 "Anchor must be set for bisection refinement of triangle");
550 // Splitting a triangle in two by bisecting anchor edge
551 // One interior edge is created
552 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
553 child_coords.col(1) = lt_node_coords.col(mod_2);
554 child_poly.push_back(child_coords);
555 break;
556 }
557 case RefPat::rp_trisect: {
558 LF_VERIFY_MSG(
559 anchor_set_,
560 "Anchor must be set for trisection refinement of triangle");
561 // Bisect through anchor edge first and then bisect through
562 // edge with the next larger index (mod 3); creates two
563 // interior edges
564 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
565 child_coords.col(1) = lt_node_coords.col(mod_2);
566 child_poly.push_back(child_coords);
567
568 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
569 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
570 child_poly.push_back(child_coords);
571 break;
572 }
573 case RefPat::rp_trisect_left: {
574 LF_VERIFY_MSG(
575 anchor_set_,
576 "Anchor must be set for trisection refinement of triangle");
577
578 // Bisect through anchor edge first and then bisect through
579 // edge with the next smaller index (mod 3); creates two
580 // interior edges
581 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
582 child_coords.col(1) = lt_node_coords.col(mod_2);
583 child_poly.push_back(child_coords);
584
585 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
586 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
587 child_poly.push_back(child_coords);
588 break;
589 }
590 case RefPat::rp_quadsect: {
591 LF_VERIFY_MSG(
592 anchor_set_,
593 "Anchor must be set for quadsection refinement of triangle");
594 // Bisect through the anchor edge first and then
595 // through the two remaining edges; creates three
596 // interior edges
597 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
598 child_coords.col(1) = lt_node_coords.col(mod_2);
599 child_poly.push_back(child_coords);
600
601 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
602 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
603 child_poly.push_back(child_coords);
604
605 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
606 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
607 child_poly.push_back(child_coords);
608 break;
609 }
610 case RefPat::rp_regular: {
611 // Split triangle into four small congruent triangles
612 // Introduces three interior edges
613 child_coords.col(0) = lt_midpoint_coords.col(0);
614 child_coords.col(1) = lt_midpoint_coords.col(2);
615 child_poly.push_back(child_coords);
616
617 child_coords.col(0) = lt_midpoint_coords.col(0);
618 child_coords.col(1) = lt_midpoint_coords.col(1);
619 child_poly.push_back(child_coords);
620
621 child_coords.col(0) = lt_midpoint_coords.col(2);
622 child_coords.col(1) = lt_midpoint_coords.col(1);
623 child_poly.push_back(child_coords);
624 break;
625 }
626 case RefPat::rp_barycentric: {
627 // Split triangle into 5 smaller triangles by connecting
628 // the center of gravity with the vertices and the midpoints
629 // of the edges. Six interior edges will be created
630
631 // Obtain coordinates of center of gravity
632 const Eigen::Vector2i lt_baryc_coords =
633 Eigen::Vector2i({lt_third, lt_third});
634
635 child_coords.col(0) = lt_node_coords.col(0);
636 child_coords.col(1) = lt_baryc_coords;
637 child_poly.push_back(child_coords);
638
639 child_coords.col(0) = lt_node_coords.col(1);
640 child_coords.col(1) = lt_baryc_coords;
641 child_poly.push_back(child_coords);
642
643 child_coords.col(0) = lt_node_coords.col(2);
644 child_coords.col(1) = lt_baryc_coords;
645 child_poly.push_back(child_coords);
646
647 child_coords.col(0) = lt_midpoint_coords.col(0);
648 child_coords.col(1) = lt_baryc_coords;
649 child_poly.push_back(child_coords);
650
651 child_coords.col(0) = lt_midpoint_coords.col(1);
652 child_coords.col(1) = lt_baryc_coords;
653 child_poly.push_back(child_coords);
654
655 child_coords.col(0) = lt_midpoint_coords.col(2);
656 child_coords.col(1) = lt_baryc_coords;
657 child_poly.push_back(child_coords);
658 break;
659 }
660 default: {
661 LF_VERIFY_MSG(false, "refinement pattern "
662 << (int)ref_pat_
663 << "not implemented for triangle!");
664 break;
665 }
666 } // end switch ref_pat
667 break;
668 } // end codim = 1
669 case 2: {
670 // Return interior points
671 if (ref_pat_ == RefPat::rp_barycentric) {
672 // Only in the case of barycentric refinemnt insert new point
673 // in the center of gravitfy.
674 child_poly.emplace_back(Eigen::Vector2i({lt_third, lt_third}));
675 }
676 break;
677 } // end codim == 2
678 default:
679 LF_VERIFY_MSG(false, "invalid codim");
680 } // end switch codim
681 break;
682 } // end case of a triangle
683 case lf::base::RefEl::kQuad(): {
684 // Lattice coordinates of special points in a quadrilateral
685 // Here we use knowledge about the conventions underlying node
686 // and edge numbering in a quadrilateral as defined in RefEl.
687 // This information could also be retrieved from the reference element.
688 Eigen::MatrixXi lt_node_coords(2, 4);
689 lt_node_coords << 0, lt_one, lt_one, 0, 0, 0, lt_one, lt_one;
690 Eigen::MatrixXi lt_midpoint_coords(2, 4);
691 lt_midpoint_coords << lt_half, lt_one, lt_half, 0, 0, lt_half, lt_one,
692 lt_half;
693
694 // Remap local indices according to anchor values
695 const unsigned int mod_0 = (0 + anchor_) % 4;
696 const unsigned int mod_1 = (1 + anchor_) % 4;
697 const unsigned int mod_2 = (2 + anchor_) % 4;
698 const unsigned int mod_3 = (3 + anchor_) % 4;
699
700 switch (codim) {
701 case 0: {
702 // Return sub-cell location inside quadrilateral
703 // For temporary storage of triangular lattice polygon
704 Eigen::MatrixXi tria_child_coords(2, 3);
705 // For temporary storage of 4-node lattice polygon
706 Eigen::MatrixXi quad_child_coords(2, 4);
707 switch (ref_pat_) {
708 case RefPat::rp_nil: {
709 break;
710 }
711 case RefPat::rp_copy: {
712 child_poly.push_back(lt_node_coords);
713 break;
714 }
715 case RefPat::rp_trisect: {
716 LF_VERIFY_MSG(
717 anchor_set_,
718 "Anchor must be set for trisection refinement of quad");
719
720 // Partition a quad into three triangle, the anchor edge
721 // being split in the process
722 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
723 tria_child_coords.col(1) = lt_node_coords.col(mod_2);
724 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
725 child_poly.push_back(tria_child_coords);
726
727 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
728 tria_child_coords.col(1) = lt_node_coords.col(mod_0);
729 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
730 child_poly.push_back(tria_child_coords);
731
732 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
733 tria_child_coords.col(1) = lt_node_coords.col(mod_1);
734 tria_child_coords.col(2) = lt_node_coords.col(mod_2);
735 child_poly.push_back(tria_child_coords);
736 break;
737 }
738 case RefPat::rp_quadsect: {
739 LF_VERIFY_MSG(
740 anchor_set_,
741 "Anchor must be set for quadsection refinement of triangle");
742 // Partition a quad into four triangle, thus
743 // splitting two edges. The one with the smaller sub index is the
744 // anchor edge
745 tria_child_coords.col(0) = lt_node_coords.col(mod_0);
746 tria_child_coords.col(1) = lt_node_coords.col(mod_3);
747 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_0);
748 child_poly.push_back(tria_child_coords);
749
750 tria_child_coords.col(0) = lt_node_coords.col(mod_1);
751 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
752 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_0);
753 child_poly.push_back(tria_child_coords);
754
755 tria_child_coords.col(0) = lt_node_coords.col(mod_2);
756 tria_child_coords.col(1) = lt_node_coords.col(mod_3);
757 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_1);
758 child_poly.push_back(tria_child_coords);
759
760 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
761 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
762 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
763 child_poly.push_back(tria_child_coords);
764 break;
765 }
766 case RefPat::rp_bisect:
767 case RefPat::rp_split: {
768 LF_VERIFY_MSG(anchor_set_,
769 "Anchor must be set for splitting of quad");
770
771 // Cut a quadrilateral into two
772 quad_child_coords.col(0) = lt_node_coords.col(mod_0);
773 quad_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
774 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_2);
775 quad_child_coords.col(3) = lt_node_coords.col(mod_3);
776 child_poly.push_back(quad_child_coords);
777
778 quad_child_coords.col(0) = lt_node_coords.col(mod_1);
779 quad_child_coords.col(1) = lt_node_coords.col(mod_2);
780 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_2);
781 quad_child_coords.col(3) = lt_midpoint_coords.col(mod_0);
782 child_poly.push_back(quad_child_coords);
783 break;
784 }
785 case RefPat::rp_threeedge: {
786 LF_VERIFY_MSG(
787 anchor_set_,
788 "Anchor must be set for three edge refinement of a quad");
789
790 // A quadrilateral with three split edges is decomposed into one
791 // child quadrilateral and three triangles
792 quad_child_coords.col(0) = lt_node_coords.col(mod_2);
793 quad_child_coords.col(1) = lt_node_coords.col(mod_3);
794 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
795 quad_child_coords.col(3) = lt_midpoint_coords.col(mod_1);
796 child_poly.push_back(quad_child_coords);
797
798 tria_child_coords.col(0) = lt_node_coords.col(mod_0);
799 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
800 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
801 child_poly.push_back(tria_child_coords);
802
803 tria_child_coords.col(0) = lt_node_coords.col(mod_1);
804 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
805 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_1);
806 child_poly.push_back(tria_child_coords);
807
808 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
809 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
810 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
811 child_poly.push_back(tria_child_coords);
812 break;
813 }
814 case RefPat::rp_barycentric:
815 case RefPat::rp_regular: {
816 // Fully symmetric splitting into four quadrilaterals
817 // Obtain coordinates of center of gravity
818 const Eigen::Vector2i lt_baryc_coords =
819 Eigen::Vector2i({lt_half, lt_half});
820
821 quad_child_coords.col(0) = lt_node_coords.col(0);
822 quad_child_coords.col(1) = lt_midpoint_coords.col(0);
823 quad_child_coords.col(2) = lt_baryc_coords;
824 quad_child_coords.col(3) = lt_midpoint_coords.col(3);
825 child_poly.push_back(quad_child_coords);
826
827 quad_child_coords.col(0) = lt_node_coords.col(1);
828 quad_child_coords.col(1) = lt_midpoint_coords.col(1);
829 quad_child_coords.col(2) = lt_baryc_coords;
830 quad_child_coords.col(3) = lt_midpoint_coords.col(0);
831 child_poly.push_back(quad_child_coords);
832
833 quad_child_coords.col(0) = lt_node_coords.col(2);
834 quad_child_coords.col(1) = lt_midpoint_coords.col(1);
835 quad_child_coords.col(2) = lt_baryc_coords;
836 quad_child_coords.col(3) = lt_midpoint_coords.col(2);
837 child_poly.push_back(quad_child_coords);
838
839 quad_child_coords.col(0) = lt_node_coords.col(3);
840 quad_child_coords.col(1) = lt_midpoint_coords.col(2);
841 quad_child_coords.col(2) = lt_baryc_coords;
842 quad_child_coords.col(3) = lt_midpoint_coords.col(3);
843 child_poly.push_back(quad_child_coords);
844 break;
845 }
846 default: {
847 LF_VERIFY_MSG(false, "refinement pattern "
848 << (int)ref_pat_
849 << "not implemented for quadrilateral!");
850 break;
851 }
852 } // end switch ref_pat
853 break;
854 } // end codim == 0
855 case 1: {
856 // Return information about edges created inside a quadrilateral
857 Eigen::MatrixXi child_coords(2, 2);
858 switch (ref_pat_) {
859 case RefPat::rp_nil:
860 case RefPat::rp_copy: {
861 // No internal edges in this case
862 break;
863 }
864 case RefPat::rp_trisect: {
865 LF_VERIFY_MSG(
866 anchor_set_,
867 "Anchor must be set for trisection refinement of quad");
868
869 // Partition a quad into three triangle, the anchor edge
870 // being split in the process. Two interior edges are created
871 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
872 child_coords.col(1) = lt_node_coords.col(mod_2);
873 child_poly.push_back(child_coords);
874
875 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
876 child_coords.col(1) = lt_node_coords.col(mod_3);
877 child_poly.push_back(child_coords);
878 break;
879 }
880 case RefPat::rp_quadsect: {
881 LF_VERIFY_MSG(
882 anchor_set_,
883 "Anchor must be set for quadsection refinement of triangle");
884 // Partition a quad into four triangle, thus
885 // splitting two edges. This spawns three internal edges
886 child_coords.col(0) = lt_node_coords.col(mod_3);
887 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
888 child_poly.push_back(child_coords);
889
890 child_coords.col(0) = lt_midpoint_coords.col(mod_1);
891 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
892 child_poly.push_back(child_coords);
893
894 child_coords.col(0) = lt_node_coords.col(mod_3);
895 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
896 child_poly.push_back(child_coords);
897 break;
898 }
899 case RefPat::rp_bisect:
900 case RefPat::rp_split: {
901 LF_VERIFY_MSG(anchor_set_,
902 "Anchor must be set for splitting of quad");
903
904 // Cut a quadrilateral into two. One interior edge arises
905 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
906 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
907 child_poly.push_back(child_coords);
908 break;
909 }
910 case RefPat::rp_threeedge: {
911 LF_VERIFY_MSG(
912 anchor_set_,
913 "Anchor must be set for three edge refinement of a quad");
914 // Split quadrilateral into one half and three triangles, creating
915 // three interior edges in the process
916 child_coords.col(0) = lt_midpoint_coords.col(mod_3);
917 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
918 child_poly.push_back(child_coords);
919
920 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
921 child_coords.col(1) = lt_midpoint_coords.col(mod_3);
922 child_poly.push_back(child_coords);
923
924 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
925 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
926 child_poly.push_back(child_coords);
927 break;
928 }
929 case RefPat::rp_barycentric:
930 case RefPat::rp_regular: {
931 // Fully symmetric splitting into four quadrilaterals
932 // Four interior edges arise
933 // Obtain coordinates of center of gravity
934 const Eigen::Vector2i lt_baryc_coords =
935 Eigen::Vector2i({lt_half, lt_half});
936
937 child_coords.col(0) = lt_midpoint_coords.col(0);
938 child_coords.col(1) = lt_baryc_coords;
939 child_poly.push_back(child_coords);
940
941 child_coords.col(0) = lt_midpoint_coords.col(1);
942 child_coords.col(1) = lt_baryc_coords;
943 child_poly.push_back(child_coords);
944
945 child_coords.col(0) = lt_midpoint_coords.col(2);
946 child_coords.col(1) = lt_baryc_coords;
947 child_poly.push_back(child_coords);
948
949 child_coords.col(0) = lt_midpoint_coords.col(3);
950 child_coords.col(1) = lt_baryc_coords;
951 child_poly.push_back(child_coords);
952 break;
953 }
954 default: {
955 LF_VERIFY_MSG(false, "refinement pattern "
956 << (int)ref_pat_
957 << "not implemented for quadrilateral!");
958 break;
959 }
960 } // end switch ref_pat
961 break;
962 }
963 case 2: {
964 // Return location of new interior points
965 // Those arise only in the case of regular/barycentric refinement
966 if ((ref_pat_ == RefPat::rp_regular) ||
967 (ref_pat_ == RefPat::rp_barycentric)) {
968 child_poly.emplace_back(Eigen::Vector2i({lt_half, lt_half}));
969 }
970 break;
971 } // end codim == 2
972 default:
973 LF_VERIFY_MSG(false, "invalid codim");
974 } // end switch codim
975 break;
976 } // end case of a quadrilateral
977 default:
978 LF_VERIFY_MSG(false, "Unsupported RefEl");
979 } // end switch cell type
980 return (child_poly);
981}
982
983} // namespace lf::refinement
lf::base::RefEl::kSegment
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition ref_el.h:153
lf::base::RefEl::Dimension
constexpr dim_t Dimension() const
Return the dimension of this reference element.
Definition ref_el.h:204
lf::base::RefEl::kPoint
static constexpr RefEl kPoint()
Returns the (0-dimensional) reference point.
Definition ref_el.h:144
lf::base::RefEl::kTria
static constexpr RefEl kTria()
Returns the reference triangle.
Definition ref_el.h:161
lf::base::RefEl::kQuad
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition ref_el.h:169
lf::base::RefEl::ToString
std::string ToString() const
Return a string representation of this Reference element.
Definition ref_el.h:458
lf::refinement::Hybrid2DRefinementPattern::NumChildren
lf::base::size_type NumChildren(lf::base::dim_t codim) const override
provide number of child entities of a given co-dimension to be created by refinement
Definition hybrid2d_refinement_pattern.cc:49
lf::refinement::Hybrid2DRefinementPattern::ChildPolygons
std::vector< Eigen::Matrix< int, Eigen::Dynamic, Eigen::Dynamic > > ChildPolygons(lf::base::dim_t codim) const override
provide lattice reference coordinates of vertices of child polygons
Definition hybrid2d_refinement_pattern.cc:274
lf::base::size_type
unsigned int size_type
general type for variables related to size of arrays
Definition types.h:20
lf::base::dim_t
unsigned int dim_t
type for dimensions and co-dimensions and numbers derived from them
Definition types.h:32
lf::refinement
tools for regular or local refinement of 2D hybrid meshes
Definition hybrid2d_refinement_pattern.cc:9
lf::refinement::operator<<
std::ostream & operator<<(std::ostream &o, const RefPat &refpat)
Definition hybrid2d_refinement_pattern.cc:11
lf::refinement::RefPat
RefPat
(possibly incomplete) list of refinement patterns for triangles/quadrilaterals
Definition hybrid2d_refinement_pattern.h:35
Generated by doxygen 1.10.0