LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
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
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 } // end switch cell type
267 return 0;
268}
269
270std::vector<Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>>
271Hybrid2DRefinementPattern::ChildPolygons(lf::base::dim_t codim) const {
272 LF_VERIFY_MSG(codim <= ref_el_.Dimension(),
273 "Illegal codim = " << codim << " for " << ref_el_.ToString());
274 // Local variable for accumulating information about interior entities
275 // created during refinement
276 std::vector<Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic>> child_poly{};
277
278 // Lattice point coordinates (lattice_const_ should be a multiple of 6)
279 const unsigned int lt_half = lattice_const_ / 2;
280 const unsigned int lt_third = lattice_const_ / 3;
281 const unsigned int lt_one = lattice_const_;
282 // Depending on the type of cell do something different
283 switch (ref_el_) {
284 case lf::base::RefEl::kPoint(): {
285 if (ref_pat_ != RefPat::rp_nil) {
286 child_poly.emplace_back(0, 1);
287 }
288 break;
289 } // end case of a simple point
290 case lf::base::RefEl::kSegment(): { // case of an edge
291 switch (codim) {
292 case 0: { // child entities are segments as well
293 switch (ref_pat_) {
294 case RefPat::rp_nil: {
295 break;
296 }
297 case RefPat::rp_copy: { // copy of segment, 1 child
298 Eigen::Matrix<int, 1, 2> seg_ref_coord;
299 seg_ref_coord << 0, lt_one;
300 child_poly.emplace_back(seg_ref_coord);
301 break;
302 }
303 case RefPat::rp_split: { // splitting of segment, 2 children
304 Eigen::Matrix<int, 1, 2> part_ref_coord;
305 part_ref_coord << 0, lt_half;
306 child_poly.emplace_back(part_ref_coord);
307 part_ref_coord << lt_half, lt_one;
308 child_poly.emplace_back(part_ref_coord);
309 break;
310 }
311 default: {
312 LF_VERIFY_MSG(false, "refinement pattern "
313 << (int)ref_pat_
314 << "not implemented for edge!");
315 break;
316 }
317 } // end switch ref_pat
318 break;
319 } // case codim = 0
320 case 1: { // child entities are points
321 if (ref_pat_ == RefPat::rp_split) {
322 // child is the midpoint
323 Eigen::Matrix<int, 1, 1> point_ref_coord;
324 point_ref_coord << lt_half;
325 child_poly.emplace_back(point_ref_coord);
326 }
327 break;
328 } // end case codim = 1
329 default:
330 LF_VERIFY_MSG(false, "invalid codim");
331 } // end switch codim
332 break;
333 } // end case of an edge
334 case lf::base::RefEl::kTria(): {
335 // **********************************************************************
336 // Triangle
337 // **********************************************************************
338 // Lattice coordinates of special points in the triangle
339 // Here we use knowledge about the conventions underlying node
340 // and edge numbering in a triangle as defined in RefEl.
341 // This information could also be retrieved from the reference element.
342 Eigen::MatrixXi lt_node_coords(2, 3);
343 lt_node_coords << 0, lt_one, 0, 0, 0, lt_one;
344 Eigen::MatrixXi lt_midpoint_coords(2, 3);
345 lt_midpoint_coords << lt_half, lt_half, 0, 0, lt_half, lt_half;
346
347 // Remap local indices according to anchor values
348 const unsigned int mod_0 = (0 + anchor_) % 3;
349 const unsigned int mod_1 = (1 + anchor_) % 3;
350 const unsigned int mod_2 = (2 + anchor_) % 3;
351
352 switch (codim) {
353 case 0: {
354 // For temporary storage of triangular lattice polygon
355 Eigen::MatrixXi child_coords(2, 3);
356
357 switch (ref_pat_) {
358 case RefPat::rp_nil: {
359 break;
360 }
361 case RefPat::rp_copy: {
362 child_coords = lt_node_coords;
363 child_poly.push_back(child_coords);
364 break;
365 }
366 case RefPat::rp_bisect: {
367 LF_VERIFY_MSG(
368 anchor_set_,
369 "Anchor must be set for bisection refinement of triangle");
370 // Splitting a triangle in two by bisecting anchor edge
371 child_coords.col(0) = lt_node_coords.col(mod_0);
372 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
373 child_coords.col(2) = lt_node_coords.col(mod_2);
374 child_poly.push_back(child_coords);
375
376 child_coords.col(0) = lt_node_coords.col(mod_1);
377 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
378 child_coords.col(2) = lt_node_coords.col(mod_2);
379 child_poly.push_back(child_coords);
380 break;
381 }
382 case RefPat::rp_trisect: {
383 LF_VERIFY_MSG(
384 anchor_set_,
385 "Anchor must be set for trisection refinement of triangle");
386 // Bisect through anchor edge first and then bisect through
387 // edge with the next larger index (mod 3); creates three
388 // child triangles.
389 child_coords.col(0) = lt_node_coords.col(mod_0);
390 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
391 child_coords.col(2) = lt_node_coords.col(mod_2);
392 child_poly.push_back(child_coords);
393
394 child_coords.col(0) = lt_node_coords.col(mod_1);
395 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
396 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
397 child_poly.push_back(child_coords);
398
399 child_coords.col(0) = lt_node_coords.col(mod_2);
400 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
401 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
402 child_poly.push_back(child_coords);
403 break;
404 }
405 case RefPat::rp_trisect_left: {
406 LF_VERIFY_MSG(
407 anchor_set_,
408 "Anchor must be set for trisection refinement of triangle");
409
410 // Bisect through anchor edge first and then bisect through
411 // edge with the next smaller index (mod 3); creates three
412 // child triangles.
413 child_coords.col(0) = lt_node_coords.col(mod_0);
414 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
415 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
416 child_poly.push_back(child_coords);
417
418 child_coords.col(0) = lt_node_coords.col(mod_1);
419 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
420 child_coords.col(2) = lt_node_coords.col(mod_2);
421 child_poly.push_back(child_coords);
422
423 child_coords.col(0) = lt_node_coords.col(mod_2);
424 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
425 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
426 child_poly.push_back(child_coords);
427 break;
428 }
429 case RefPat::rp_quadsect: {
430 LF_VERIFY_MSG(
431 anchor_set_,
432 "Anchor must be set for quadsection refinement of triangle");
433 // Bisect through the anchor edge first and then
434 // through the two remaining edges; creates four child
435 // triangles; every edge is split.
436 child_coords.col(0) = lt_node_coords.col(mod_0);
437 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
438 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
439 child_poly.push_back(child_coords);
440
441 child_coords.col(0) = lt_node_coords.col(mod_1);
442 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
443 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
444 child_poly.push_back(child_coords);
445
446 child_coords.col(0) = lt_node_coords.col(mod_2);
447 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
448 child_coords.col(2) = lt_midpoint_coords.col(mod_1);
449 child_poly.push_back(child_coords);
450
451 child_coords.col(0) = lt_node_coords.col(mod_2);
452 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
453 child_coords.col(2) = lt_midpoint_coords.col(mod_2);
454 child_poly.push_back(child_coords);
455 break;
456 }
457 case RefPat::rp_regular: {
458 // Split triangle into four small congruent triangles
459 child_coords.col(0) = lt_node_coords.col(0);
460 child_coords.col(1) = lt_midpoint_coords.col(0);
461 child_coords.col(2) = lt_midpoint_coords.col(2);
462 child_poly.push_back(child_coords);
463
464 child_coords.col(0) = lt_node_coords.col(1);
465 child_coords.col(1) = lt_midpoint_coords.col(0);
466 child_coords.col(2) = lt_midpoint_coords.col(1);
467 child_poly.push_back(child_coords);
468
469 child_coords.col(0) = lt_node_coords.col(2);
470 child_coords.col(1) = lt_midpoint_coords.col(2);
471 child_coords.col(2) = lt_midpoint_coords.col(1);
472 child_poly.push_back(child_coords);
473
474 child_coords.col(0) = lt_midpoint_coords.col(0);
475 child_coords.col(1) = lt_midpoint_coords.col(1);
476 child_coords.col(2) = lt_midpoint_coords.col(2);
477 child_poly.push_back(child_coords);
478 break;
479 }
480 case RefPat::rp_barycentric: {
481 // Split triangle into 5 smaller triangles by connecting
482 // the center of gravity with the vertices and the midpoints
483 // of the edges.
484
485 // Obtain coordinates of center of gravity
486 Eigen::Vector2i lt_baryc_coords =
487 Eigen::Vector2i({lt_third, lt_third});
488
489 child_coords.col(0) = lt_node_coords.col(0);
490 child_coords.col(1) = lt_midpoint_coords.col(0);
491 child_coords.col(2) = lt_baryc_coords;
492 child_poly.push_back(child_coords);
493
494 child_coords.col(0) = lt_node_coords.col(1);
495 child_coords.col(1) = lt_midpoint_coords.col(0);
496 child_coords.col(2) = lt_baryc_coords;
497 child_poly.push_back(child_coords);
498
499 child_coords.col(0) = lt_node_coords.col(1);
500 child_coords.col(1) = lt_midpoint_coords.col(1);
501 child_coords.col(2) = lt_baryc_coords;
502 child_poly.push_back(child_coords);
503
504 child_coords.col(0) = lt_node_coords.col(2);
505 child_coords.col(1) = lt_midpoint_coords.col(1);
506 child_coords.col(2) = lt_baryc_coords;
507 child_poly.push_back(child_coords);
508
509 child_coords.col(0) = lt_node_coords.col(2);
510 child_coords.col(1) = lt_midpoint_coords.col(2);
511 child_coords.col(2) = lt_baryc_coords;
512 child_poly.push_back(child_coords);
513
514 child_coords.col(0) = lt_node_coords.col(0);
515 child_coords.col(1) = lt_midpoint_coords.col(2);
516 child_coords.col(2) = lt_baryc_coords;
517 child_poly.push_back(child_coords);
518 break;
519 }
520 default: {
521 LF_VERIFY_MSG(false, "refinement pattern "
522 << (int)ref_pat_
523 << "not implemented for triangle!");
524 break;
525 }
526 } // end switch ref_pat
527 break;
528 } // end case codim = 0
529 case 1: { // codim == 1
530 // Return the location of interior edges of the triangle
531 // This location is defined by two lattice points given by the
532 // rows of the following matrix
533 Eigen::MatrixXi child_coords(2, 2);
534
535 // Different layout of interior edges for different refinement
536 // patterns
537 switch (ref_pat_) {
538 case RefPat::rp_nil:
539 case RefPat::rp_copy: {
540 // No interior edges in these cases
541 break;
542 }
543 case RefPat::rp_bisect: {
544 LF_VERIFY_MSG(
545 anchor_set_,
546 "Anchor must be set for bisection refinement of triangle");
547 // Splitting a triangle in two by bisecting anchor edge
548 // One interior edge is created
549 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
550 child_coords.col(1) = lt_node_coords.col(mod_2);
551 child_poly.push_back(child_coords);
552 break;
553 }
554 case RefPat::rp_trisect: {
555 LF_VERIFY_MSG(
556 anchor_set_,
557 "Anchor must be set for trisection refinement of triangle");
558 // Bisect through anchor edge first and then bisect through
559 // edge with the next larger index (mod 3); creates two
560 // interior edges
561 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
562 child_coords.col(1) = lt_node_coords.col(mod_2);
563 child_poly.push_back(child_coords);
564
565 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
566 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
567 child_poly.push_back(child_coords);
568 break;
569 }
570 case RefPat::rp_trisect_left: {
571 LF_VERIFY_MSG(
572 anchor_set_,
573 "Anchor must be set for trisection refinement of triangle");
574
575 // Bisect through anchor edge first and then bisect through
576 // edge with the next smaller index (mod 3); creates two
577 // interior edges
578 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
579 child_coords.col(1) = lt_node_coords.col(mod_2);
580 child_poly.push_back(child_coords);
581
582 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
583 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
584 child_poly.push_back(child_coords);
585 break;
586 }
587 case RefPat::rp_quadsect: {
588 LF_VERIFY_MSG(
589 anchor_set_,
590 "Anchor must be set for quadsection refinement of triangle");
591 // Bisect through the anchor edge first and then
592 // through the two remaining edges; creates three
593 // interior edges
594 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
595 child_coords.col(1) = lt_node_coords.col(mod_2);
596 child_poly.push_back(child_coords);
597
598 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
599 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
600 child_poly.push_back(child_coords);
601
602 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
603 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
604 child_poly.push_back(child_coords);
605 break;
606 }
607 case RefPat::rp_regular: {
608 // Split triangle into four small congruent triangles
609 // Introduces three interior edges
610 child_coords.col(0) = lt_midpoint_coords.col(0);
611 child_coords.col(1) = lt_midpoint_coords.col(2);
612 child_poly.push_back(child_coords);
613
614 child_coords.col(0) = lt_midpoint_coords.col(0);
615 child_coords.col(1) = lt_midpoint_coords.col(1);
616 child_poly.push_back(child_coords);
617
618 child_coords.col(0) = lt_midpoint_coords.col(2);
619 child_coords.col(1) = lt_midpoint_coords.col(1);
620 child_poly.push_back(child_coords);
621 break;
622 }
623 case RefPat::rp_barycentric: {
624 // Split triangle into 5 smaller triangles by connecting
625 // the center of gravity with the vertices and the midpoints
626 // of the edges. Six interior edges will be created
627
628 // Obtain coordinates of center of gravity
629 Eigen::Vector2i lt_baryc_coords =
630 Eigen::Vector2i({lt_third, lt_third});
631
632 child_coords.col(0) = lt_node_coords.col(0);
633 child_coords.col(1) = lt_baryc_coords;
634 child_poly.push_back(child_coords);
635
636 child_coords.col(0) = lt_node_coords.col(1);
637 child_coords.col(1) = lt_baryc_coords;
638 child_poly.push_back(child_coords);
639
640 child_coords.col(0) = lt_node_coords.col(2);
641 child_coords.col(1) = lt_baryc_coords;
642 child_poly.push_back(child_coords);
643
644 child_coords.col(0) = lt_midpoint_coords.col(0);
645 child_coords.col(1) = lt_baryc_coords;
646 child_poly.push_back(child_coords);
647
648 child_coords.col(0) = lt_midpoint_coords.col(1);
649 child_coords.col(1) = lt_baryc_coords;
650 child_poly.push_back(child_coords);
651
652 child_coords.col(0) = lt_midpoint_coords.col(2);
653 child_coords.col(1) = lt_baryc_coords;
654 child_poly.push_back(child_coords);
655 break;
656 }
657 default: {
658 LF_VERIFY_MSG(false, "refinement pattern "
659 << (int)ref_pat_
660 << "not implemented for triangle!");
661 break;
662 }
663 } // end switch ref_pat
664 break;
665 } // end codim = 1
666 case 2: {
667 // Return interior points
668 if (ref_pat_ == RefPat::rp_barycentric) {
669 // Only in the case of barycentric refinemnt insert new point
670 // in the center of gravitfy.
671 child_poly.emplace_back(Eigen::Vector2i({lt_third, lt_third}));
672 }
673 break;
674 } // end codim == 2
675 default:
676 LF_VERIFY_MSG(false, "invalid codim");
677 } // end switch codim
678 break;
679 } // end case of a triangle
680 case lf::base::RefEl::kQuad(): {
681 // Lattice coordinates of special points in a quadrilateral
682 // Here we use knowledge about the conventions underlying node
683 // and edge numbering in a quadrilateral as defined in RefEl.
684 // This information could also be retrieved from the reference element.
685 Eigen::MatrixXi lt_node_coords(2, 4);
686 lt_node_coords << 0, lt_one, lt_one, 0, 0, 0, lt_one, lt_one;
687 Eigen::MatrixXi lt_midpoint_coords(2, 4);
688 lt_midpoint_coords << lt_half, lt_one, lt_half, 0, 0, lt_half, lt_one,
689 lt_half;
690
691 // Remap local indices according to anchor values
692 const unsigned int mod_0 = (0 + anchor_) % 4;
693 const unsigned int mod_1 = (1 + anchor_) % 4;
694 const unsigned int mod_2 = (2 + anchor_) % 4;
695 const unsigned int mod_3 = (3 + anchor_) % 4;
696
697 switch (codim) {
698 case 0: {
699 // Return sub-cell location inside quadrilateral
700 // For temporary storage of triangular lattice polygon
701 Eigen::MatrixXi tria_child_coords(2, 3);
702 // For temporary storage of 4-node lattice polygon
703 Eigen::MatrixXi quad_child_coords(2, 4);
704 switch (ref_pat_) {
705 case RefPat::rp_nil: {
706 break;
707 }
708 case RefPat::rp_copy: {
709 child_poly.push_back(lt_node_coords);
710 break;
711 }
712 case RefPat::rp_trisect: {
713 LF_VERIFY_MSG(
714 anchor_set_,
715 "Anchor must be set for trisection refinement of quad");
716
717 // Partition a quad into three triangle, the anchor edge
718 // being split in the process
719 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
720 tria_child_coords.col(1) = lt_node_coords.col(mod_2);
721 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
722 child_poly.push_back(tria_child_coords);
723
724 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
725 tria_child_coords.col(1) = lt_node_coords.col(mod_0);
726 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
727 child_poly.push_back(tria_child_coords);
728
729 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
730 tria_child_coords.col(1) = lt_node_coords.col(mod_1);
731 tria_child_coords.col(2) = lt_node_coords.col(mod_2);
732 child_poly.push_back(tria_child_coords);
733 break;
734 }
735 case RefPat::rp_quadsect: {
736 LF_VERIFY_MSG(
737 anchor_set_,
738 "Anchor must be set for quadsection refinement of triangle");
739 // Partition a quad into four triangle, thus
740 // splitting two edges. The one with the smaller sub index is the
741 // anchor edge
742 tria_child_coords.col(0) = lt_node_coords.col(mod_0);
743 tria_child_coords.col(1) = lt_node_coords.col(mod_3);
744 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_0);
745 child_poly.push_back(tria_child_coords);
746
747 tria_child_coords.col(0) = lt_node_coords.col(mod_1);
748 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
749 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_0);
750 child_poly.push_back(tria_child_coords);
751
752 tria_child_coords.col(0) = lt_node_coords.col(mod_2);
753 tria_child_coords.col(1) = lt_node_coords.col(mod_3);
754 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_1);
755 child_poly.push_back(tria_child_coords);
756
757 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
758 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
759 tria_child_coords.col(2) = lt_node_coords.col(mod_3);
760 child_poly.push_back(tria_child_coords);
761 break;
762 }
763 case RefPat::rp_bisect:
764 case RefPat::rp_split: {
765 LF_VERIFY_MSG(anchor_set_,
766 "Anchor must be set for splitting of quad");
767
768 // Cut a quadrilateral into two
769 quad_child_coords.col(0) = lt_node_coords.col(mod_0);
770 quad_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
771 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_2);
772 quad_child_coords.col(3) = lt_node_coords.col(mod_3);
773 child_poly.push_back(quad_child_coords);
774
775 quad_child_coords.col(0) = lt_node_coords.col(mod_1);
776 quad_child_coords.col(1) = lt_node_coords.col(mod_2);
777 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_2);
778 quad_child_coords.col(3) = lt_midpoint_coords.col(mod_0);
779 child_poly.push_back(quad_child_coords);
780 break;
781 }
782 case RefPat::rp_threeedge: {
783 LF_VERIFY_MSG(
784 anchor_set_,
785 "Anchor must be set for three edge refinement of a quad");
786
787 // A quadrilateral with three split edges is decomposed into one
788 // child quadrilateral and three triangles
789 quad_child_coords.col(0) = lt_node_coords.col(mod_2);
790 quad_child_coords.col(1) = lt_node_coords.col(mod_3);
791 quad_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
792 quad_child_coords.col(3) = lt_midpoint_coords.col(mod_1);
793 child_poly.push_back(quad_child_coords);
794
795 tria_child_coords.col(0) = lt_node_coords.col(mod_0);
796 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
797 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
798 child_poly.push_back(tria_child_coords);
799
800 tria_child_coords.col(0) = lt_node_coords.col(mod_1);
801 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_0);
802 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_1);
803 child_poly.push_back(tria_child_coords);
804
805 tria_child_coords.col(0) = lt_midpoint_coords.col(mod_0);
806 tria_child_coords.col(1) = lt_midpoint_coords.col(mod_1);
807 tria_child_coords.col(2) = lt_midpoint_coords.col(mod_3);
808 child_poly.push_back(tria_child_coords);
809 break;
810 }
811 case RefPat::rp_barycentric:
812 case RefPat::rp_regular: {
813 // Fully symmetric splitting into four quadrilaterals
814 // Obtain coordinates of center of gravity
815 Eigen::Vector2i lt_baryc_coords =
816 Eigen::Vector2i({lt_half, lt_half});
817
818 quad_child_coords.col(0) = lt_node_coords.col(0);
819 quad_child_coords.col(1) = lt_midpoint_coords.col(0);
820 quad_child_coords.col(2) = lt_baryc_coords;
821 quad_child_coords.col(3) = lt_midpoint_coords.col(3);
822 child_poly.push_back(quad_child_coords);
823
824 quad_child_coords.col(0) = lt_node_coords.col(1);
825 quad_child_coords.col(1) = lt_midpoint_coords.col(1);
826 quad_child_coords.col(2) = lt_baryc_coords;
827 quad_child_coords.col(3) = lt_midpoint_coords.col(0);
828 child_poly.push_back(quad_child_coords);
829
830 quad_child_coords.col(0) = lt_node_coords.col(2);
831 quad_child_coords.col(1) = lt_midpoint_coords.col(1);
832 quad_child_coords.col(2) = lt_baryc_coords;
833 quad_child_coords.col(3) = lt_midpoint_coords.col(2);
834 child_poly.push_back(quad_child_coords);
835
836 quad_child_coords.col(0) = lt_node_coords.col(3);
837 quad_child_coords.col(1) = lt_midpoint_coords.col(2);
838 quad_child_coords.col(2) = lt_baryc_coords;
839 quad_child_coords.col(3) = lt_midpoint_coords.col(3);
840 child_poly.push_back(quad_child_coords);
841 break;
842 }
843 default: {
844 LF_VERIFY_MSG(false, "refinement pattern "
845 << (int)ref_pat_
846 << "not implemented for quadrilateral!");
847 break;
848 }
849 } // end switch ref_pat
850 break;
851 } // end codim == 0
852 case 1: {
853 // Return information about edges created inside a quadrilateral
854 Eigen::MatrixXi child_coords(2, 2);
855 switch (ref_pat_) {
856 case RefPat::rp_nil:
857 case RefPat::rp_copy: {
858 // No internal edges in this case
859 break;
860 }
861 case RefPat::rp_trisect: {
862 LF_VERIFY_MSG(
863 anchor_set_,
864 "Anchor must be set for trisection refinement of quad");
865
866 // Partition a quad into three triangle, the anchor edge
867 // being split in the process. Two interior edges are created
868 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
869 child_coords.col(1) = lt_node_coords.col(mod_2);
870 child_poly.push_back(child_coords);
871
872 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
873 child_coords.col(1) = lt_node_coords.col(mod_3);
874 child_poly.push_back(child_coords);
875 break;
876 }
877 case RefPat::rp_quadsect: {
878 LF_VERIFY_MSG(
879 anchor_set_,
880 "Anchor must be set for quadsection refinement of triangle");
881 // Partition a quad into four triangle, thus
882 // splitting two edges. This spawns three internal edges
883 child_coords.col(0) = lt_node_coords.col(mod_3);
884 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
885 child_poly.push_back(child_coords);
886
887 child_coords.col(0) = lt_midpoint_coords.col(mod_1);
888 child_coords.col(1) = lt_midpoint_coords.col(mod_0);
889 child_poly.push_back(child_coords);
890
891 child_coords.col(0) = lt_node_coords.col(mod_3);
892 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
893 child_poly.push_back(child_coords);
894 break;
895 }
896 case RefPat::rp_bisect:
897 case RefPat::rp_split: {
898 LF_VERIFY_MSG(anchor_set_,
899 "Anchor must be set for splitting of quad");
900
901 // Cut a quadrilateral into two. One interior edge arises
902 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
903 child_coords.col(1) = lt_midpoint_coords.col(mod_2);
904 child_poly.push_back(child_coords);
905 break;
906 }
907 case RefPat::rp_threeedge: {
908 LF_VERIFY_MSG(
909 anchor_set_,
910 "Anchor must be set for three edge refinement of a quad");
911 // Split quadrilateral into one half and three triangles, creating
912 // three interior edges in the process
913 child_coords.col(0) = lt_midpoint_coords.col(mod_3);
914 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
915 child_poly.push_back(child_coords);
916
917 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
918 child_coords.col(1) = lt_midpoint_coords.col(mod_3);
919 child_poly.push_back(child_coords);
920
921 child_coords.col(0) = lt_midpoint_coords.col(mod_0);
922 child_coords.col(1) = lt_midpoint_coords.col(mod_1);
923 child_poly.push_back(child_coords);
924 break;
925 }
926 case RefPat::rp_barycentric:
927 case RefPat::rp_regular: {
928 // Fully symmetric splitting into four quadrilaterals
929 // Four interior edges arise
930 // Obtain coordinates of center of gravity
931 Eigen::Vector2i lt_baryc_coords =
932 Eigen::Vector2i({lt_half, lt_half});
933
934 child_coords.col(0) = lt_midpoint_coords.col(0);
935 child_coords.col(1) = lt_baryc_coords;
936 child_poly.push_back(child_coords);
937
938 child_coords.col(0) = lt_midpoint_coords.col(1);
939 child_coords.col(1) = lt_baryc_coords;
940 child_poly.push_back(child_coords);
941
942 child_coords.col(0) = lt_midpoint_coords.col(2);
943 child_coords.col(1) = lt_baryc_coords;
944 child_poly.push_back(child_coords);
945
946 child_coords.col(0) = lt_midpoint_coords.col(3);
947 child_coords.col(1) = lt_baryc_coords;
948 child_poly.push_back(child_coords);
949 break;
950 }
951 default: {
952 LF_VERIFY_MSG(false, "refinement pattern "
953 << (int)ref_pat_
954 << "not implemented for quadrilateral!");
955 break;
956 }
957 } // end switch ref_pat
958 break;
959 }
960 case 2: {
961 // Return location of new interior points
962 // Those arise only in the case of regular/barycentric refinement
963 if ((ref_pat_ == RefPat::rp_regular) ||
964 (ref_pat_ == RefPat::rp_barycentric)) {
965 child_poly.emplace_back(Eigen::Vector2i({lt_half, lt_half}));
966 }
967 break;
968 } // end codim == 2
969 default:
970 LF_VERIFY_MSG(false, "invalid codim");
971 } // end switch codim
972 break;
973 } // end case of a quadrilateral
974 } // end switch cell type
975 return (child_poly);
976}
977
978} // namespace lf::refinement
lf::base::RefEl::kSegment
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition: ref_el.h:150
lf::base::RefEl::Dimension
constexpr dim_t Dimension() const
Return the dimension of this reference element.
Definition: ref_el.h:201
lf::base::RefEl::kPoint
static constexpr RefEl kPoint()
Returns the (0-dimensional) reference point.
Definition: ref_el.h:141
lf::base::RefEl::kTria
static constexpr RefEl kTria()
Returns the reference triangle.
Definition: ref_el.h:158
lf::base::RefEl::kQuad
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition: ref_el.h:166
lf::base::RefEl::ToString
std::string ToString() const
Return a string representation of this Reference element.
Definition: ref_el.h:455
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:271
lf::base::dim_t
unsigned int dim_t
type for dimensions and co-dimensions and numbers derived from them
Definition: base.h:36
lf::base::size_type
unsigned int size_type
general type for variables related to size of arrays
Definition: base.h:24
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.9.3