FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
MeshQuad4.h
Go to the documentation of this file.
1
10#ifndef GOOSEFEM_MESHQUAD4_H
11#define GOOSEFEM_MESHQUAD4_H
12
13#include "Mesh.h"
14#include "config.h"
15
16namespace GooseFEM {
17namespace Mesh {
18
22namespace Quad4 {
23
24// pre-allocation
25namespace Map {
27}
28
32class Regular : public RegularBase2d<Regular> {
33public:
34 Regular() = default;
35
43 Regular(size_t nelx, size_t nely, double h = 1.0)
44 {
45 m_h = h;
46 m_nelx = nelx;
47 m_nely = nely;
48 m_ndim = 2;
49 m_nne = 4;
50
51 GOOSEFEM_ASSERT(m_nelx >= 1);
52 GOOSEFEM_ASSERT(m_nely >= 1);
53
54 m_nnode = (m_nelx + 1) * (m_nely + 1);
55 m_nelem = m_nelx * m_nely;
56 }
57
64 {
65 return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
66 }
67
68private:
69 friend class RegularBase<Regular>;
70 friend class RegularBase2d<Regular>;
71
72 size_t nelx_impl() const
73 {
74 return m_nelx;
75 }
76
77 size_t nely_impl() const
78 {
79 return m_nely;
80 }
81
82 ElementType getElementType_impl() const
83 {
84 return ElementType::Quad4;
85 }
86
87 array_type::tensor<double, 2> coor_impl() const
88 {
89 array_type::tensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
90
91 array_type::tensor<double, 1> x =
92 xt::linspace<double>(0.0, m_h * static_cast<double>(m_nelx), m_nelx + 1);
93 array_type::tensor<double, 1> y =
94 xt::linspace<double>(0.0, m_h * static_cast<double>(m_nely), m_nely + 1);
95
96 size_t inode = 0;
97
98 for (size_t iy = 0; iy < m_nely + 1; ++iy) {
99 for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
100 ret(inode, 0) = x(ix);
101 ret(inode, 1) = y(iy);
102 ++inode;
103 }
104 }
105
106 return ret;
107 }
108
109 array_type::tensor<size_t, 2> conn_impl() const
110 {
111 array_type::tensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
112
113 size_t ielem = 0;
114
115 for (size_t iy = 0; iy < m_nely; ++iy) {
116 for (size_t ix = 0; ix < m_nelx; ++ix) {
117 ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix);
118 ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1);
119 ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix);
120 ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1);
121 ++ielem;
122 }
123 }
124
125 return ret;
126 }
127
128 array_type::tensor<size_t, 1> nodesBottomEdge_impl() const
129 {
130 return xt::arange<size_t>(m_nelx + 1);
131 }
132
133 array_type::tensor<size_t, 1> nodesTopEdge_impl() const
134 {
135 return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
136 }
137
138 array_type::tensor<size_t, 1> nodesLeftEdge_impl() const
139 {
140 return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
141 }
142
143 array_type::tensor<size_t, 1> nodesRightEdge_impl() const
144 {
145 return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
146 }
147
148 array_type::tensor<size_t, 1> nodesBottomOpenEdge_impl() const
149 {
150 return xt::arange<size_t>(1, m_nelx);
151 }
152
153 array_type::tensor<size_t, 1> nodesTopOpenEdge_impl() const
154 {
155 return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
156 }
157
158 array_type::tensor<size_t, 1> nodesLeftOpenEdge_impl() const
159 {
160 return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
161 }
162
163 array_type::tensor<size_t, 1> nodesRightOpenEdge_impl() const
164 {
165 return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
166 }
167
168 size_t nodesBottomLeftCorner_impl() const
169 {
170 return 0;
171 }
172
173 size_t nodesBottomRightCorner_impl() const
174 {
175 return m_nelx;
176 }
177
178 size_t nodesTopLeftCorner_impl() const
179 {
180 return m_nely * (m_nelx + 1);
181 }
182
183 size_t nodesTopRightCorner_impl() const
184 {
185 return m_nely * (m_nelx + 1) + m_nelx;
186 }
187
188 double m_h;
189 size_t m_nelx;
190 size_t m_nely;
191 size_t m_nelem;
192 size_t m_nnode;
193 size_t m_nne;
194 size_t m_ndim;
195};
196
200class FineLayer : public RegularBase2d<FineLayer> {
201public:
202 FineLayer() = default;
203
215 FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1)
216 {
217 this->init(nelx, nely, h, nfine);
218 }
219
228 template <class C, class E, std::enable_if_t<xt::is_xexpression<C>::value, bool> = true>
229 FineLayer(const C& coor, const E& conn)
230 {
231 this->init_by_mapping(coor, conn);
232 }
233
242 {
243 return m_nhx;
244 }
245
254 {
255 return m_nhy;
256 }
257
266 {
267 return m_refine;
268 }
269
278 {
279 return m_layer_nelx;
280 }
281
288 {
289 size_t nely = m_nhy.size();
290 size_t iy = (nely - 1) / 2;
291 return m_startElem(iy) + xt::arange<size_t>(m_layer_nelx(iy));
292 }
293
300 {
301 GOOSEFEM_ASSERT(layer < m_layer_nelx.size());
302 size_t n = m_layer_nelx(layer);
303 if (m_refine(layer) != -1) {
304 n *= 4;
305 }
306 return m_startElem(layer) + xt::arange<size_t>(n);
307 }
308
315 std::vector<size_t> start_stop_rows,
316 std::vector<size_t> start_stop_cols) const
317 {
318 GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2);
319 GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2);
320
321 std::array<size_t, 2> rows;
322 std::array<size_t, 2> cols;
323
324 if (start_stop_rows.size() == 2) {
325 std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin());
326 GOOSEFEM_ASSERT(rows[0] <= this->nely());
327 GOOSEFEM_ASSERT(rows[1] <= this->nely());
328 }
329 else {
330 rows[0] = 0;
331 rows[1] = this->nely();
332 }
333
334 if (start_stop_cols.size() == 2) {
335 std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin());
336 GOOSEFEM_ASSERT(cols[0] <= this->nelx());
337 GOOSEFEM_ASSERT(cols[1] <= this->nelx());
338 }
339 else {
340 cols[0] = 0;
341 cols[1] = this->nelx();
342 }
343
344 if (rows[0] == rows[1] || cols[0] == cols[1]) {
345 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({0});
346 return ret;
347 }
348
349 // Compute dimensions
350
351 auto H = xt::cumsum(m_nhy);
352 size_t yl = 0;
353 if (rows[0] > 0) {
354 yl = xt::argmax(H > rows[0])();
355 }
356 size_t yu = xt::argmax(H >= rows[1])();
357 size_t hx = std::max(m_nhx(yl), m_nhx(yu));
358 size_t xl = (cols[0] - cols[0] % hx) / hx;
359 size_t xu = (cols[1] - cols[1] % hx) / hx;
360
361 // Allocate output
362
363 size_t N = 0;
364
365 for (size_t i = yl; i <= yu; ++i) {
366 // no refinement
367 size_t n = (xu - xl) * hx / m_nhx(i);
368 // refinement
369 if (m_refine(i) != -1) {
370 n *= 4;
371 }
372 N += n;
373 }
374
375 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({N});
376
377 // Write output
378
379 N = 0;
380
381 for (size_t i = yl; i <= yu; ++i) {
382 // no refinement
383 size_t n = (xu - xl) * hx / m_nhx(i);
384 size_t h = hx;
385 // refinement
386 if (m_refine(i) != -1) {
387 n *= 4;
388 h *= 4;
389 }
391 m_startElem(i) + xl * h / m_nhx(i) + xt::arange<size_t>(n);
392 xt::view(ret, xt::range(N, N + n)) = e;
393 N += n;
394 }
395
396 return ret;
397 }
398
409 elementgrid_around_ravel(size_t e, size_t size, bool periodic = true)
410 {
411 GOOSEFEM_WIP_ASSERT(periodic == true);
412
413 size_t iy = xt::argmin(m_startElem <= e)() - 1;
414 size_t nel = m_layer_nelx(iy);
415
416 GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
417
418 auto H = xt::cumsum(m_nhy);
419
420 if (2 * size >= H(H.size() - 1)) {
421 return xt::arange<size_t>(this->nelem());
422 }
423
424 size_t hy = H(iy);
425 size_t l = xt::argmax(H > (hy - size - 1))();
426 size_t u = xt::argmax(H >= (hy + size))();
427 size_t lh = 0;
428 if (l > 0) {
429 lh = H(l - 1);
430 }
431 size_t uh = H(u);
432
433 size_t step = xt::amax(m_nhx)();
434 size_t relx = (e - m_startElem(iy)) % step;
435 size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
436 size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
437 size_t dx = m_nhx(u);
438 size_t xl = 0;
439 size_t xu = nel;
440 if (mid >= size) {
441 xl = mid - size;
442 }
443 if (mid + size < nel) {
444 xu = mid + size + 1;
445 }
446 xl = xl - xl % dx;
447 xu = xu - xu % dx;
448 if (mid - xl < size) {
449 if (xl < dx) {
450 xl = 0;
451 }
452 else {
453 xl -= dx;
454 }
455 }
456 if (xu - mid < size) {
457 if (xu > nel - dx) {
458 xu = nel;
459 }
460 else {
461 xu += dx;
462 }
463 }
464
465 auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu});
466 auto map = this->roll(nroll);
467 return xt::view(map, xt::keep(ret));
468 }
469
480 // -
482 elementgrid_leftright(size_t e, size_t left, size_t right, bool periodic = true)
483 {
484 GOOSEFEM_WIP_ASSERT(periodic == true);
485
486 size_t iy = xt::argmin(m_startElem <= e)() - 1;
487 size_t nel = m_layer_nelx(iy);
488
489 GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
490
491 size_t step = xt::amax(m_nhx)();
492 size_t relx = (e - m_startElem(iy)) % step;
493 size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
494 size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
495 size_t dx = m_nhx(iy);
496 size_t xl = 0;
497 size_t xu = nel;
498 if (mid >= left) {
499 xl = mid - left;
500 }
501 if (mid + right < nel) {
502 xu = mid + right + 1;
503 }
504 xl = xl - xl % dx;
505 xu = xu - xu % dx;
506 if (mid - xl < left) {
507 if (xl < dx) {
508 xl = 0;
509 }
510 else {
511 xl -= dx;
512 }
513 }
514 if (xu - mid < right) {
515 if (xu > nel - dx) {
516 xu = nel;
517 }
518 else {
519 xu += dx;
520 }
521 }
522
523 auto H = xt::cumsum(m_nhy);
524 size_t yl = 0;
525 if (iy > 0) {
526 yl = H(iy - 1);
527 }
528 auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu});
529 auto map = this->roll(nroll);
530 return xt::view(map, xt::keep(ret));
531 }
532
539 {
540 auto conn = this->conn();
541 size_t nely = static_cast<size_t>(m_nhy.size());
542 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({m_nelem});
543
544 // loop over all element layers
545 for (size_t iy = 0; iy < nely; ++iy) {
546
547 // no refinement
548 size_t shift = n * (m_layer_nelx(iy) / m_layer_nelx(0));
549 size_t nel = m_layer_nelx(iy);
550
551 // refinement
552 if (m_refine(iy) != -1) {
553 shift = n * (m_layer_nelx(iy) / m_layer_nelx(0)) * 4;
554 nel = m_layer_nelx(iy) * 4;
555 }
556
557 // element numbers of the layer, and roll them
558 auto e = m_startElem(iy) + xt::arange<size_t>(nel);
559 xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift);
560 }
561
562 return ret;
563 }
564
565private:
566 friend class RegularBase<FineLayer>;
567 friend class RegularBase2d<FineLayer>;
569
570 size_t nelx_impl() const
571 {
572 return xt::amax(m_layer_nelx)();
573 }
574
575 size_t nely_impl() const
576 {
577 return xt::sum(m_nhy)();
578 }
579
580 ElementType getElementType_impl() const
581 {
582 return ElementType::Quad4;
583 }
584
585 array_type::tensor<double, 2> coor_impl() const
586 {
587 // allocate output
588 array_type::tensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
589
590 // current node, number of element layers
591 size_t inode = 0;
592 size_t nely = static_cast<size_t>(m_nhy.size());
593
594 // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
595 // - allocate
596 array_type::tensor<double, 1> y = xt::empty<double>({nely + 1});
597 // - initialize
598 y(0) = 0.0;
599 // - compute
600 for (size_t iy = 1; iy < nely + 1; ++iy) {
601 y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h;
602 }
603
604 // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of
605 // nodes
606
607 for (size_t iy = 0;; ++iy) {
608 // get positions along the x- and z-axis
609 array_type::tensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
610
611 // add nodes of the bottom layer of this element
612 for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
613 ret(inode, 0) = x(ix);
614 ret(inode, 1) = y(iy);
615 ++inode;
616 }
617
618 // stop at middle layer
619 if (iy == (nely - 1) / 2) {
620 break;
621 }
622
623 // add extra nodes of the intermediate layer, for refinement in x-direction
624 if (m_refine(iy) == 0) {
625 // - get position offset in x- and y-direction
626 double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
627 double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
628 // - add nodes of the intermediate layer
629 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
630 for (size_t j = 0; j < 2; ++j) {
631 ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
632 ret(inode, 1) = y(iy) + dy;
633 ++inode;
634 }
635 }
636 }
637 }
638
639 // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes
640
641 for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
642 // get positions along the x- and z-axis
643 array_type::tensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
644
645 // add extra nodes of the intermediate layer, for refinement in x-direction
646 if (m_refine(iy) == 0) {
647 // - get position offset in x- and y-direction
648 double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
649 double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
650 // - add nodes of the intermediate layer
651 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
652 for (size_t j = 0; j < 2; ++j) {
653 ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
654 ret(inode, 1) = y(iy) + dy;
655 ++inode;
656 }
657 }
658 }
659
660 // add nodes of the top layer of this element
661 for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
662 ret(inode, 0) = x(ix);
663 ret(inode, 1) = y(iy + 1);
664 ++inode;
665 }
666 }
667
668 return ret;
669 }
670
671 array_type::tensor<size_t, 2> conn_impl() const
672 {
673 // allocate output
674 array_type::tensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
675
676 // current element, number of element layers, starting nodes of each node layer
677 size_t ielem = 0;
678 size_t nely = static_cast<size_t>(m_nhy.size());
679 size_t bot, mid, top;
680
681 // loop over all element layers
682 for (size_t iy = 0; iy < nely; ++iy) {
683 // - get: starting nodes of bottom(, middle) and top layer
684 bot = m_startNode(iy);
685 mid = m_startNode(iy) + m_nnd(iy);
686 top = m_startNode(iy + 1);
687
688 // - define connectivity: no coarsening/refinement
689 if (m_refine(iy) == -1) {
690 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
691 ret(ielem, 0) = bot + (ix);
692 ret(ielem, 1) = bot + (ix + 1);
693 ret(ielem, 2) = top + (ix + 1);
694 ret(ielem, 3) = top + (ix);
695 ielem++;
696 }
697 }
698
699 // - define connectivity: refinement along the x-direction (below the middle layer)
700 else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
701 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
702 // -- bottom element
703 ret(ielem, 0) = bot + (ix);
704 ret(ielem, 1) = bot + (ix + 1);
705 ret(ielem, 2) = mid + (2 * ix + 1);
706 ret(ielem, 3) = mid + (2 * ix);
707 ielem++;
708 // -- top-right element
709 ret(ielem, 0) = bot + (ix + 1);
710 ret(ielem, 1) = top + (3 * ix + 3);
711 ret(ielem, 2) = top + (3 * ix + 2);
712 ret(ielem, 3) = mid + (2 * ix + 1);
713 ielem++;
714 // -- top-center element
715 ret(ielem, 0) = mid + (2 * ix);
716 ret(ielem, 1) = mid + (2 * ix + 1);
717 ret(ielem, 2) = top + (3 * ix + 2);
718 ret(ielem, 3) = top + (3 * ix + 1);
719 ielem++;
720 // -- top-left element
721 ret(ielem, 0) = bot + (ix);
722 ret(ielem, 1) = mid + (2 * ix);
723 ret(ielem, 2) = top + (3 * ix + 1);
724 ret(ielem, 3) = top + (3 * ix);
725 ielem++;
726 }
727 }
728
729 // - define connectivity: coarsening along the x-direction (above the middle layer)
730 else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) {
731 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
732 // -- lower-left element
733 ret(ielem, 0) = bot + (3 * ix);
734 ret(ielem, 1) = bot + (3 * ix + 1);
735 ret(ielem, 2) = mid + (2 * ix);
736 ret(ielem, 3) = top + (ix);
737 ielem++;
738 // -- lower-center element
739 ret(ielem, 0) = bot + (3 * ix + 1);
740 ret(ielem, 1) = bot + (3 * ix + 2);
741 ret(ielem, 2) = mid + (2 * ix + 1);
742 ret(ielem, 3) = mid + (2 * ix);
743 ielem++;
744 // -- lower-right element
745 ret(ielem, 0) = bot + (3 * ix + 2);
746 ret(ielem, 1) = bot + (3 * ix + 3);
747 ret(ielem, 2) = top + (ix + 1);
748 ret(ielem, 3) = mid + (2 * ix + 1);
749 ielem++;
750 // -- upper element
751 ret(ielem, 0) = mid + (2 * ix);
752 ret(ielem, 1) = mid + (2 * ix + 1);
753 ret(ielem, 2) = top + (ix + 1);
754 ret(ielem, 3) = top + (ix);
755 ielem++;
756 }
757 }
758 }
759
760 return ret;
761 }
762
763 array_type::tensor<size_t, 1> nodesBottomEdge_impl() const
764 {
765 return m_startNode(0) + xt::arange<size_t>(m_layer_nelx(0) + 1);
766 }
767
768 array_type::tensor<size_t, 1> nodesTopEdge_impl() const
769 {
770 size_t nely = m_nhy.size();
771 return m_startNode(nely) + xt::arange<size_t>(m_layer_nelx(nely - 1) + 1);
772 }
773
774 array_type::tensor<size_t, 1> nodesLeftEdge_impl() const
775 {
776 size_t nely = m_nhy.size();
777
778 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
779
780 size_t i = 0;
781 size_t j = (nely + 1) / 2;
782 size_t k = (nely - 1) / 2;
783 size_t l = nely;
784
785 xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j));
786 xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1));
787
788 return ret;
789 }
790
791 array_type::tensor<size_t, 1> nodesRightEdge_impl() const
792 {
793 size_t nely = m_nhy.size();
794
795 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
796
797 size_t i = 0;
798 size_t j = (nely + 1) / 2;
799 size_t k = (nely - 1) / 2;
800 size_t l = nely;
801
802 xt::view(ret, xt::range(i, j)) =
803 xt::view(m_startNode, xt::range(i, j)) + xt::view(m_layer_nelx, xt::range(i, j));
804
805 xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1)) +
806 xt::view(m_layer_nelx, xt::range(k, l));
807
808 return ret;
809 }
810
811 array_type::tensor<size_t, 1> nodesBottomOpenEdge_impl() const
812 {
813 return m_startNode(0) + xt::arange<size_t>(1, m_layer_nelx(0));
814 }
815
816 array_type::tensor<size_t, 1> nodesTopOpenEdge_impl() const
817 {
818 size_t nely = m_nhy.size();
819
820 return m_startNode(nely) + xt::arange<size_t>(1, m_layer_nelx(nely - 1));
821 }
822
823 array_type::tensor<size_t, 1> nodesLeftOpenEdge_impl() const
824 {
825 size_t nely = m_nhy.size();
826
827 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
828
829 size_t i = 0;
830 size_t j = (nely + 1) / 2;
831 size_t k = (nely - 1) / 2;
832 size_t l = nely;
833
834 xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j));
835 xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l));
836
837 return ret;
838 }
839
840 array_type::tensor<size_t, 1> nodesRightOpenEdge_impl() const
841 {
842 size_t nely = m_nhy.size();
843
844 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
845
846 size_t i = 0;
847 size_t j = (nely + 1) / 2;
848 size_t k = (nely - 1) / 2;
849 size_t l = nely;
850
851 xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j)) +
852 xt::view(m_layer_nelx, xt::range(i + 1, j));
853
854 xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l)) +
855 xt::view(m_layer_nelx, xt::range(k, l - 1));
856
857 return ret;
858 }
859
860 size_t nodesBottomLeftCorner_impl() const
861 {
862 return m_startNode(0);
863 }
864
865 size_t nodesBottomRightCorner_impl() const
866 {
867 return m_startNode(0) + m_layer_nelx(0);
868 }
869
870 size_t nodesTopLeftCorner_impl() const
871 {
872 size_t nely = m_nhy.size();
873 return m_startNode(nely);
874 }
875
876 size_t nodesTopRightCorner_impl() const
877 {
878 size_t nely = m_nhy.size();
879 return m_startNode(nely) + m_layer_nelx(nely - 1);
880 }
881
882 double m_h;
883 size_t m_nelem;
884 size_t m_nnode;
885 size_t m_nne;
886 size_t m_ndim;
887 double m_Lx;
888 array_type::tensor<size_t, 1> m_layer_nelx;
889 array_type::tensor<size_t, 1> m_nhx;
890 array_type::tensor<size_t, 1> m_nhy;
891 array_type::tensor<size_t, 1> m_nnd;
892 array_type::tensor<int, 1> m_refine;
893 array_type::tensor<size_t, 1> m_startElem;
894 array_type::tensor<size_t, 1> m_startNode;
895
899 void init(size_t nelx, size_t nely, double h, size_t nfine = 1)
900 {
901 GOOSEFEM_ASSERT(nelx >= 1ul);
902 GOOSEFEM_ASSERT(nely >= 1ul);
903
904 m_h = h;
905 m_ndim = 2;
906 m_nne = 4;
907 m_Lx = m_h * static_cast<double>(nelx);
908
909 // compute element size in y-direction (use symmetry, compute upper half)
910
911 // temporary variables
912 size_t nmin, ntot;
913 array_type::tensor<size_t, 1> nhx = xt::ones<size_t>({nely});
914 array_type::tensor<size_t, 1> nhy = xt::ones<size_t>({nely});
915 array_type::tensor<int, 1> refine = -1 * xt::ones<int>({nely});
916
917 // minimum height in y-direction (half of the height because of symmetry)
918 if (nely % 2 == 0) {
919 nmin = nely / 2;
920 }
921 else {
922 nmin = (nely + 1) / 2;
923 }
924
925 // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half)
926 if (nfine % 2 == 0) {
927 nfine = nfine / 2 + 1;
928 }
929 else {
930 nfine = (nfine + 1) / 2;
931 }
932
933 if (nfine < 1) {
934 nfine = 1;
935 }
936
937 if (nfine > nmin) {
938 nfine = nmin;
939 }
940
941 // loop over element layers in y-direction, try to coarsen using these rules:
942 // (1) element size in y-direction <= distance to origin in y-direction
943 // (2) element size in x-direction should fit the total number of elements in x-direction
944 // (3) a certain number of layers have the minimum size "1" (are fine)
945 for (size_t iy = nfine;;) {
946 // initialize current size in y-direction
947 if (iy == nfine) {
948 ntot = nfine;
949 }
950 // check to stop
951 if (iy >= nely || ntot >= nmin) {
952 nely = iy;
953 break;
954 }
955
956 // rules (1,2) satisfied: coarsen in x-direction
957 if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
958 refine(iy) = 0;
959 nhy(iy) *= 2;
960 auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
961 auto vnhx = xt::view(nhx, xt::range(iy, _));
962 vnhy *= 3;
963 vnhx *= 3;
964 }
965
966 // update the number of elements in y-direction
967 ntot += nhy(iy);
968 // proceed to next element layer in y-direction
969 ++iy;
970 // check to stop
971 if (iy >= nely || ntot >= nmin) {
972 nely = iy;
973 break;
974 }
975 }
976
977 // symmetrize, compute full information
978
979 // allocate mesh constructor parameters
980 m_nhx = xt::empty<size_t>({nely * 2 - 1});
981 m_nhy = xt::empty<size_t>({nely * 2 - 1});
982 m_refine = xt::empty<int>({nely * 2 - 1});
983 m_layer_nelx = xt::empty<size_t>({nely * 2 - 1});
984 m_nnd = xt::empty<size_t>({nely * 2});
985 m_startElem = xt::empty<size_t>({nely * 2 - 1});
986 m_startNode = xt::empty<size_t>({nely * 2});
987
988 // fill
989 // - lower half
990 for (size_t iy = 0; iy < nely; ++iy) {
991 m_nhx(iy) = nhx(nely - iy - 1);
992 m_nhy(iy) = nhy(nely - iy - 1);
993 m_refine(iy) = refine(nely - iy - 1);
994 }
995 // - upper half
996 for (size_t iy = 0; iy < nely - 1; ++iy) {
997 m_nhx(iy + nely) = nhx(iy + 1);
998 m_nhy(iy + nely) = nhy(iy + 1);
999 m_refine(iy + nely) = refine(iy + 1);
1000 }
1001
1002 // update size
1003 nely = m_nhx.size();
1004
1005 // compute the number of elements per element layer in y-direction
1006 for (size_t iy = 0; iy < nely; ++iy) {
1007 m_layer_nelx(iy) = nelx / m_nhx(iy);
1008 }
1009
1010 // compute the number of nodes per node layer in y-direction
1011 for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
1012 m_nnd(iy) = m_layer_nelx(iy) + 1;
1013 }
1014 for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
1015 m_nnd(iy + 1) = m_layer_nelx(iy) + 1;
1016 }
1017
1018 // compute mesh dimensions
1019
1020 // initialize
1021 m_nnode = 0;
1022 m_nelem = 0;
1023 m_startNode(0) = 0;
1024
1025 // loop over element layers (bottom -> middle, elements become finer)
1026 for (size_t i = 0; i < (nely - 1) / 2; ++i) {
1027 // - store the first element of the layer
1028 m_startElem(i) = m_nelem;
1029 // - add the nodes of this layer
1030 if (m_refine(i) == 0) {
1031 m_nnode += (3 * m_layer_nelx(i) + 1);
1032 }
1033 else {
1034 m_nnode += (m_layer_nelx(i) + 1);
1035 }
1036 // - add the elements of this layer
1037 if (m_refine(i) == 0) {
1038 m_nelem += (4 * m_layer_nelx(i));
1039 }
1040 else {
1041 m_nelem += (m_layer_nelx(i));
1042 }
1043 // - store the starting node of the next layer
1044 m_startNode(i + 1) = m_nnode;
1045 }
1046
1047 // loop over element layers (middle -> top, elements become coarser)
1048 for (size_t i = (nely - 1) / 2; i < nely; ++i) {
1049 // - store the first element of the layer
1050 m_startElem(i) = m_nelem;
1051 // - add the nodes of this layer
1052 if (m_refine(i) == 0) {
1053 m_nnode += (5 * m_layer_nelx(i) + 1);
1054 }
1055 else {
1056 m_nnode += (m_layer_nelx(i) + 1);
1057 }
1058 // - add the elements of this layer
1059 if (m_refine(i) == 0) {
1060 m_nelem += (4 * m_layer_nelx(i));
1061 }
1062 else {
1063 m_nelem += (m_layer_nelx(i));
1064 }
1065 // - store the starting node of the next layer
1066 m_startNode(i + 1) = m_nnode;
1067 }
1068 // - add the top row of nodes
1069 m_nnode += m_layer_nelx(nely - 1) + 1;
1070 }
1071
1075 template <class C, class E>
1076 void init_by_mapping(const C& coor, const E& conn)
1077 {
1078 GOOSEFEM_ASSERT(coor.dimension() == 2);
1079 GOOSEFEM_ASSERT(conn.dimension() == 2);
1080 GOOSEFEM_ASSERT(coor.shape(1) == 2);
1081 GOOSEFEM_ASSERT(conn.shape(1) == 4);
1082 GOOSEFEM_ASSERT(conn.shape(0) > 0);
1083 GOOSEFEM_ASSERT(coor.shape(0) >= 4);
1084 GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
1085
1086 if (conn.shape(0) == 1) {
1087 this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0));
1088 GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
1089 GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
1090 return;
1091 }
1092
1093 // Identify the middle layer
1094
1095 size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2;
1096 size_t eleft = emid;
1097 size_t eright = emid;
1098
1099 while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) {
1100 eleft--;
1101 }
1102
1103 while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) {
1104 eright++;
1105 }
1106
1107 GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0));
1108
1109 // Get element sizes along the middle layer
1110
1111 auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0);
1112 auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1);
1113 auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2);
1114 auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0);
1115 auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1);
1116 auto hx = xt::amin(dx)();
1117 auto hy = xt::amin(dy)();
1118
1119 GOOSEFEM_CHECK(xt::allclose(hx, hy));
1120 GOOSEFEM_CHECK(xt::allclose(dx, hx));
1121 GOOSEFEM_CHECK(xt::allclose(dy, hy));
1122
1123 // Extract shape and initialise
1124
1125 size_t nelx = eright - eleft + 1;
1126 size_t nely = static_cast<size_t>((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx);
1127 this->init(nelx, nely, hx);
1128 GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
1129 GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
1131 xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange<size_t>(nelx))));
1132 }
1133};
1134
1138namespace Map {
1139
1144public:
1145 RefineRegular() = default;
1146
1155 : m_coarse(mesh), m_nx(nx), m_ny(ny)
1156 {
1157 m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h());
1158
1159 array_type::tensor<size_t, 2> elmat_coarse = m_coarse.elementgrid();
1160 array_type::tensor<size_t, 2> elmat_fine = m_fine.elementgrid();
1161
1162 m_coarse2fine = xt::empty<size_t>({m_coarse.nelem(), nx * ny});
1163
1164 for (size_t i = 0; i < elmat_coarse.shape(0); ++i) {
1165 for (size_t j = 0; j < elmat_coarse.shape(1); ++j) {
1166 xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view(
1167 elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx)));
1168 }
1169 }
1170 }
1171
1177 size_t nx() const
1178 {
1179 return m_nx;
1180 }
1181
1187 size_t ny() const
1188 {
1189 return m_ny;
1190 }
1191
1197 {
1198 return m_coarse;
1199 }
1200
1206 {
1207 return m_fine;
1208 }
1209
1215 {
1216 return m_coarse2fine;
1217 }
1218
1224 {
1225 return m_coarse;
1226 }
1227
1233 {
1234 return m_fine;
1235 }
1236
1241 [[deprecated]] const array_type::tensor<size_t, 2>& getMap() const
1242 {
1243 return m_coarse2fine;
1244 }
1245
1254 template <class T, size_t rank>
1256 {
1257 GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
1258
1259 std::array<size_t, rank> shape;
1260 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1261 shape[0] = m_coarse2fine.shape(0);
1262
1263 array_type::tensor<T, rank> ret = xt::empty<T>(shape);
1264
1265 for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
1266 auto e = xt::view(m_coarse2fine, i, xt::all());
1267 auto d = xt::view(data, xt::keep(e));
1268 xt::view(ret, i) = xt::mean(d, 0);
1269 }
1270
1271 return ret;
1272 }
1273
1284 template <class T, size_t rank, class S>
1286 const array_type::tensor<T, rank>& data,
1287 const array_type::tensor<S, rank>& weights) const
1288 {
1289 GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
1290
1291 std::array<size_t, rank> shape;
1292 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1293 shape[0] = m_coarse2fine.shape(0);
1294
1295 array_type::tensor<T, rank> ret = xt::empty<T>(shape);
1296
1297 for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
1298 auto e = xt::view(m_coarse2fine, i, xt::all());
1299 array_type::tensor<T, rank> d = xt::view(data, xt::keep(e));
1300 array_type::tensor<T, rank> w = xt::view(weights, xt::keep(e));
1301 xt::view(ret, i) = xt::average(d, w, {0});
1302 }
1303
1304 return ret;
1305 }
1306
1319 template <class T, size_t rank>
1321 {
1322 GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
1323
1324 std::array<size_t, rank> shape;
1325 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1326 shape[0] = m_coarse2fine.size();
1327
1328 array_type::tensor<T, rank> ret = xt::empty<T>(shape);
1329
1330 for (size_t e = 0; e < m_coarse2fine.shape(0); ++e) {
1331 for (size_t i = 0; i < m_coarse2fine.shape(1); ++i) {
1332 xt::view(ret, m_coarse2fine(e, i)) = xt::view(data, e);
1333 }
1334 }
1335
1336 return ret;
1337 }
1338
1339private:
1342 size_t m_nx;
1343 size_t m_ny;
1344 array_type::tensor<size_t, 2> m_coarse2fine;
1345};
1346
1353public:
1354 FineLayer2Regular() = default;
1355
1362 {
1363 // ------------
1364 // Regular-mesh
1365 // ------------
1366
1368 xt::amax(m_finelayer.m_layer_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);
1369
1370 // -------
1371 // mapping
1372 // -------
1373
1374 // allocate mapping
1375 m_elem_regular.resize(m_finelayer.m_nelem);
1376 m_frac_regular.resize(m_finelayer.m_nelem);
1377
1378 // alias
1379 array_type::tensor<size_t, 1> nhx = m_finelayer.m_nhx;
1380 array_type::tensor<size_t, 1> nhy = m_finelayer.m_nhy;
1381 array_type::tensor<size_t, 1> nelx = m_finelayer.m_layer_nelx;
1382 array_type::tensor<size_t, 1> start = m_finelayer.m_startElem;
1383
1384 // 'matrix' of element numbers of the Regular-mesh
1385 array_type::tensor<size_t, 2> elementgrid = m_regular.elementgrid();
1386
1387 // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh
1389 xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));
1390
1391 // number of element layers in y-direction of the FineLayer-mesh
1392 size_t nely = nhy.size();
1393
1394 // loop over layers of the FineLayer-mesh
1395 for (size_t iy = 0; iy < nely; ++iy) {
1396 // element numbers of the Regular-mesh along this layer of the FineLayer-mesh
1397 auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());
1398
1399 // no coarsening/refinement
1400 // ------------------------
1401
1402 if (m_finelayer.m_refine(iy) == -1) {
1403 // element numbers of the FineLayer-mesh along this layer
1404 array_type::tensor<size_t, 1> el_old = start(iy) + xt::arange<size_t>(nelx(iy));
1405
1406 // loop along this layer of the FineLayer-mesh
1407 for (size_t ix = 0; ix < nelx(iy); ++ix) {
1408 // get the element numbers of the Regular-mesh for this element of the
1409 // FineLayer-mesh
1410 auto block =
1411 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
1412
1413 // write to mapping
1414 for (auto& i : block) {
1415 m_elem_regular[el_old(ix)].push_back(i);
1416 m_frac_regular[el_old(ix)].push_back(1.0);
1417 }
1418 }
1419 }
1420
1421 // refinement along the x-direction (below the middle layer)
1422 // ---------------------------------------------------------
1423
1424 else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
1425 // element numbers of the FineLayer-mesh along this layer
1426 // rows: coarse block, columns element numbers per block
1428 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
1429
1430 // loop along this layer of the FineLayer-mesh
1431 for (size_t ix = 0; ix < nelx(iy); ++ix) {
1432 // get the element numbers of the Regular-mesh for this block of the
1433 // FineLayer-mesh
1434 auto block =
1435 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
1436
1437 // bottom: wide-to-narrow
1438 {
1439 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1440 auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));
1441
1442 m_elem_regular[el_old(ix, 0)].push_back(e(0));
1443 m_frac_regular[el_old(ix, 0)].push_back(0.5);
1444
1445 for (size_t k = 1; k < e.size() - 1; ++k) {
1446 m_elem_regular[el_old(ix, 0)].push_back(e(k));
1447 m_frac_regular[el_old(ix, 0)].push_back(1.0);
1448 }
1449
1450 m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
1451 m_frac_regular[el_old(ix, 0)].push_back(0.5);
1452 }
1453 }
1454
1455 // top: regular small element
1456 {
1457 auto e = xt::view(
1458 block,
1459 xt::range(nhy(iy) / 2, nhy(iy)),
1460 xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
1461
1462 for (auto& i : e) {
1463 m_elem_regular[el_old(ix, 2)].push_back(i);
1464 m_frac_regular[el_old(ix, 2)].push_back(1.0);
1465 }
1466 }
1467
1468 // left
1469 {
1470 // left-bottom: narrow-to-wide
1471 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1472 auto e = xt::view(block, j, xt::range(0, j + 1));
1473
1474 for (size_t k = 0; k < e.size() - 1; ++k) {
1475 m_elem_regular[el_old(ix, 3)].push_back(e(k));
1476 m_frac_regular[el_old(ix, 3)].push_back(1.0);
1477 }
1478
1479 m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
1480 m_frac_regular[el_old(ix, 3)].push_back(0.5);
1481 }
1482
1483 // left-top: regular
1484 {
1485 auto e = xt::view(
1486 block,
1487 xt::range(nhy(iy) / 2, nhy(iy)),
1488 xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
1489
1490 for (auto& i : e) {
1491 m_elem_regular[el_old(ix, 3)].push_back(i);
1492 m_frac_regular[el_old(ix, 3)].push_back(1.0);
1493 }
1494 }
1495 }
1496
1497 // right
1498 {
1499 // right-bottom: narrow-to-wide
1500 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1501 auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy)));
1502
1503 m_elem_regular[el_old(ix, 1)].push_back(e(0));
1504 m_frac_regular[el_old(ix, 1)].push_back(0.5);
1505
1506 for (size_t k = 1; k < e.size(); ++k) {
1507 m_elem_regular[el_old(ix, 1)].push_back(e(k));
1508 m_frac_regular[el_old(ix, 1)].push_back(1.0);
1509 }
1510 }
1511
1512 // right-top: regular
1513 {
1514 auto e = xt::view(
1515 block,
1516 xt::range(nhy(iy) / 2, nhy(iy)),
1517 xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
1518
1519 for (auto& i : e) {
1520 m_elem_regular[el_old(ix, 1)].push_back(i);
1521 m_frac_regular[el_old(ix, 1)].push_back(1.0);
1522 }
1523 }
1524 }
1525 }
1526 }
1527
1528 // coarsening along the x-direction (above the middle layer)
1529 else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
1530 // element numbers of the FineLayer-mesh along this layer
1531 // rows: coarse block, columns element numbers per block
1533 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
1534
1535 // loop along this layer of the FineLayer-mesh
1536 for (size_t ix = 0; ix < nelx(iy); ++ix) {
1537 // get the element numbers of the Regular-mesh for this block of the
1538 // FineLayer-mesh
1539 auto block =
1540 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
1541
1542 // top: narrow-to-wide
1543 {
1544 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1545 auto e = xt::view(
1546 block,
1547 nhy(iy) / 2 + j,
1548 xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));
1549
1550 m_elem_regular[el_old(ix, 3)].push_back(e(0));
1551 m_frac_regular[el_old(ix, 3)].push_back(0.5);
1552
1553 for (size_t k = 1; k < e.size() - 1; ++k) {
1554 m_elem_regular[el_old(ix, 3)].push_back(e(k));
1555 m_frac_regular[el_old(ix, 3)].push_back(1.0);
1556 }
1557
1558 m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
1559 m_frac_regular[el_old(ix, 3)].push_back(0.5);
1560 }
1561 }
1562
1563 // bottom: regular small element
1564 {
1565 auto e = xt::view(
1566 block,
1567 xt::range(0, nhy(iy) / 2),
1568 xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
1569
1570 for (auto& i : e) {
1571 m_elem_regular[el_old(ix, 1)].push_back(i);
1572 m_frac_regular[el_old(ix, 1)].push_back(1.0);
1573 }
1574 }
1575
1576 // left
1577 {
1578 // left-bottom: regular
1579 {
1580 auto e = xt::view(
1581 block,
1582 xt::range(0, nhy(iy) / 2),
1583 xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
1584
1585 for (auto& i : e) {
1586 m_elem_regular[el_old(ix, 0)].push_back(i);
1587 m_frac_regular[el_old(ix, 0)].push_back(1.0);
1588 }
1589 }
1590
1591 // left-top: narrow-to-wide
1592 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1593 auto e =
1594 xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));
1595
1596 for (size_t k = 0; k < e.size() - 1; ++k) {
1597 m_elem_regular[el_old(ix, 0)].push_back(e(k));
1598 m_frac_regular[el_old(ix, 0)].push_back(1.0);
1599 }
1600
1601 m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
1602 m_frac_regular[el_old(ix, 0)].push_back(0.5);
1603 }
1604 }
1605
1606 // right
1607 {
1608 // right-bottom: regular
1609 {
1610 auto e = xt::view(
1611 block,
1612 xt::range(0, nhy(iy) / 2),
1613 xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
1614
1615 for (auto& i : e) {
1616 m_elem_regular[el_old(ix, 2)].push_back(i);
1617 m_frac_regular[el_old(ix, 2)].push_back(1.0);
1618 }
1619 }
1620
1621 // right-top: narrow-to-wide
1622 for (size_t j = 0; j < nhy(iy) / 2; ++j) {
1623 auto e = xt::view(
1624 block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));
1625
1626 m_elem_regular[el_old(ix, 2)].push_back(e(0));
1627 m_frac_regular[el_old(ix, 2)].push_back(0.5);
1628
1629 for (size_t k = 1; k < e.size(); ++k) {
1630 m_elem_regular[el_old(ix, 2)].push_back(e(k));
1631 m_frac_regular[el_old(ix, 2)].push_back(1.0);
1632 }
1633 }
1634 }
1635 }
1636 }
1637 }
1638 }
1639
1646 {
1647 return m_regular;
1648 }
1649
1656 {
1657 return m_finelayer;
1658 }
1659
1660 // elements of the Regular mesh per element of the FineLayer mesh
1661 // and the fraction by which the overlap is
1662
1669 std::vector<std::vector<size_t>> map() const
1670 {
1671 return m_elem_regular;
1672 }
1673
1679 std::vector<std::vector<double>> mapFraction() const
1680 {
1681 return m_frac_regular;
1682 }
1683
1690 {
1691 return m_regular;
1692 }
1693
1700 {
1701 return m_finelayer;
1702 }
1703
1704 // elements of the Regular mesh per element of the FineLayer mesh
1705 // and the fraction by which the overlap is
1706
1713 [[deprecated]] std::vector<std::vector<size_t>> getMap() const
1714 {
1715 return m_elem_regular;
1716 }
1717
1723 [[deprecated]] std::vector<std::vector<double>> getMapFraction() const
1724 {
1725 return m_frac_regular;
1726 }
1727
1741 template <class T, size_t rank>
1743 {
1744 GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
1745
1746 std::array<size_t, rank> shape;
1747 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1748 shape[0] = m_regular.nelem();
1749
1750 array_type::tensor<T, rank> ret = xt::zeros<T>(shape);
1751
1752 for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
1753 for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
1754 xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
1755 }
1756 }
1757
1758 return ret;
1759 }
1760
1761private:
1764 std::vector<std::vector<size_t>> m_elem_regular;
1765 std::vector<std::vector<double>> m_frac_regular;
1766};
1767
1768} // namespace Map
1769
1770} // namespace Quad4
1771} // namespace Mesh
1772} // namespace GooseFEM
1773
1774#endif
Generic mesh operations.
Mesh with fine middle layer, and coarser elements towards the top and bottom.
Definition: MeshQuad4.h:200
const array_type::tensor< size_t, 1 > & elemrow_nhy() const
Edge size in y-direction of a block, in units of h, per row of blocks.
Definition: MeshQuad4.h:253
array_type::tensor< size_t, 1 > roll(size_t n)
Mapping to 'roll' periodically in the x-direction,.
Definition: MeshQuad4.h:538
array_type::tensor< size_t, 1 > elementsLayer(size_t layer) const
Elements along a layer.
Definition: MeshQuad4.h:299
array_type::tensor< size_t, 1 > elementgrid_leftright(size_t e, size_t left, size_t right, bool periodic=true)
Select region of elements from 'matrix' of element numbers around an element: left/right from element...
Definition: MeshQuad4.h:482
array_type::tensor< size_t, 1 > elementgrid_around_ravel(size_t e, size_t size, bool periodic=true)
Select region of elements from 'matrix' of element numbers around an element: square box with edge-si...
Definition: MeshQuad4.h:409
FineLayer(const C &coor, const E &conn)
Reconstruct class for given coordinates / connectivity.
Definition: MeshQuad4.h:229
FineLayer(size_t nelx, size_t nely, double h=1.0, size_t nfine=1)
Constructor.
Definition: MeshQuad4.h:215
const array_type::tensor< int, 1 > & elemrow_type() const
Per row of blocks: -1: normal layer 0: transition layer to match coarse and finer element on the prev...
Definition: MeshQuad4.h:265
const array_type::tensor< size_t, 1 > & elemrow_nhx() const
Edge size in x-direction of a block, in units of h, per row of blocks.
Definition: MeshQuad4.h:241
array_type::tensor< size_t, 1 > elementgrid_ravel(std::vector< size_t > start_stop_rows, std::vector< size_t > start_stop_cols) const
Select region of elements from 'matrix' of element numbers.
Definition: MeshQuad4.h:314
const array_type::tensor< size_t, 1 > & elemrow_nelem() const
Number of elements per row of blocks.
Definition: MeshQuad4.h:277
array_type::tensor< size_t, 1 > elementsMiddleLayer() const
Elements in the middle (fine) layer.
Definition: MeshQuad4.h:287
Map a FineLayer mesh to a Regular mesh.
Definition: MeshQuad4.h:1352
array_type::tensor< T, rank > mapToRegular(const array_type::tensor< T, rank > &data) const
Map element quantities to Regular.
Definition: MeshQuad4.h:1742
std::vector< std::vector< double > > getMapFraction() const
To overlap fraction for each item in the mapping in getMap().
Definition: MeshQuad4.h:1723
GooseFEM::Mesh::Quad4::FineLayer fineLayerMesh() const
Obtain the FineLayer mesh (copy of the mesh passed to the constructor).
Definition: MeshQuad4.h:1655
std::vector< std::vector< size_t > > map() const
Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh.
Definition: MeshQuad4.h:1669
GooseFEM::Mesh::Quad4::Regular regularMesh() const
Obtain the Regular mesh.
Definition: MeshQuad4.h:1645
GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const
Obtain the FineLayer mesh (copy of the mesh passed to the constructor).
Definition: MeshQuad4.h:1699
GooseFEM::Mesh::Quad4::Regular getRegularMesh() const
Obtain the Regular mesh.
Definition: MeshQuad4.h:1689
std::vector< std::vector< double > > mapFraction() const
To overlap fraction for each item in the mapping in map().
Definition: MeshQuad4.h:1679
std::vector< std::vector< size_t > > getMap() const
Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh.
Definition: MeshQuad4.h:1713
FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer &mesh)
Constructors.
Definition: MeshQuad4.h:1361
Refine a Regular mesh: subdivide elements in several smaller elements.
Definition: MeshQuad4.h:1143
GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const
Obtain the coarse mesh (copy of the mesh passed to the constructor).
Definition: MeshQuad4.h:1223
const array_type::tensor< size_t, 2 > & getMap() const
Get element-mapping: elements of the fine mesh per element of the coarse mesh.
Definition: MeshQuad4.h:1241
array_type::tensor< T, rank > mapToFine(const array_type::tensor< T, rank > &data) const
Map element quantities to the fine mesh.
Definition: MeshQuad4.h:1320
GooseFEM::Mesh::Quad4::Regular fineMesh() const
Obtain the fine mesh.
Definition: MeshQuad4.h:1205
GooseFEM::Mesh::Quad4::Regular coarseMesh() const
Obtain the coarse mesh (copy of the mesh passed to the constructor).
Definition: MeshQuad4.h:1196
size_t ny() const
For each coarse element: number of fine elements in y-direction.
Definition: MeshQuad4.h:1187
array_type::tensor< T, rank > meanToCoarse(const array_type::tensor< T, rank > &data) const
Compute the mean of the quantity define on the fine mesh when mapped on the coarse mesh.
Definition: MeshQuad4.h:1255
RefineRegular(const GooseFEM::Mesh::Quad4::Regular &mesh, size_t nx, size_t ny)
Constructor.
Definition: MeshQuad4.h:1154
size_t nx() const
For each coarse element: number of fine elements in x-direction.
Definition: MeshQuad4.h:1177
array_type::tensor< T, rank > averageToCoarse(const array_type::tensor< T, rank > &data, const array_type::tensor< S, rank > &weights) const
Compute the average of the quantity define on the fine mesh when mapped on the coarse mesh.
Definition: MeshQuad4.h:1285
GooseFEM::Mesh::Quad4::Regular getFineMesh() const
Obtain the fine mesh.
Definition: MeshQuad4.h:1232
const array_type::tensor< size_t, 2 > & map() const
Get element-mapping: elements of the fine mesh per element of the coarse mesh.
Definition: MeshQuad4.h:1214
Regular mesh: equi-sized elements.
Definition: MeshQuad4.h:32
Regular(size_t nelx, size_t nely, double h=1.0)
Constructor.
Definition: MeshQuad4.h:43
array_type::tensor< size_t, 2 > elementgrid() const
Element numbers as 'matrix'.
Definition: MeshQuad4.h:63
CRTP base class for regular meshes in 2d.
Definition: Mesh.h:339
CRTP base class for regular meshes.
Definition: Mesh.h:181
auto nely() const
Number of elements in y-direction == height of the mesh, in units of h,.
Definition: Mesh.h:237
auto coor() const
Nodal coordinates [nnode, ndim].
Definition: Mesh.h:264
auto h() const
Linear edge size of one 'block'.
Definition: Mesh.h:246
auto nelx() const
Number of elements in x-direction == width of the mesh in units of h.
Definition: Mesh.h:228
auto nelem() const
Number of elements.
Definition: Mesh.h:192
auto conn() const
Connectivity [nelem, nne].
Definition: Mesh.h:273
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
Definition: config.h:96
#define GOOSEFEM_WIP_ASSERT(expr)
Assertion that concerns temporary implementation limitations.
Definition: config.h:116
#define GOOSEFEM_CHECK(expr)
Assertion that cannot be switched off.
Definition: config.h:106
ElementType
Enumerator for element-types.
Definition: Mesh.h:31
@ Quad4
Quadrilateral: 4-noded element in 2-d.
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:176
Toolbox to perform finite element computations.
Definition: Allocate.h:14