10#ifndef GOOSEFEM_MESHQUAD4_H
11#define GOOSEFEM_MESHQUAD4_H
54 m_nnode = (m_nelx + 1) * (m_nely + 1);
55 m_nelem = m_nelx * m_nely;
65 return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
72 size_t nelx_impl()
const
77 size_t nely_impl()
const
87 array_type::tensor<double, 2> coor_impl()
const
89 array_type::tensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
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);
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);
109 array_type::tensor<size_t, 2> conn_impl()
const
111 array_type::tensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
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);
128 array_type::tensor<size_t, 1> nodesBottomEdge_impl()
const
130 return xt::arange<size_t>(m_nelx + 1);
133 array_type::tensor<size_t, 1> nodesTopEdge_impl()
const
135 return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
138 array_type::tensor<size_t, 1> nodesLeftEdge_impl()
const
140 return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
143 array_type::tensor<size_t, 1> nodesRightEdge_impl()
const
145 return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
148 array_type::tensor<size_t, 1> nodesBottomOpenEdge_impl()
const
150 return xt::arange<size_t>(1, m_nelx);
153 array_type::tensor<size_t, 1> nodesTopOpenEdge_impl()
const
155 return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
158 array_type::tensor<size_t, 1> nodesLeftOpenEdge_impl()
const
160 return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
163 array_type::tensor<size_t, 1> nodesRightOpenEdge_impl()
const
165 return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
168 size_t nodesBottomLeftCorner_impl()
const
173 size_t nodesBottomRightCorner_impl()
const
178 size_t nodesTopLeftCorner_impl()
const
180 return m_nely * (m_nelx + 1);
183 size_t nodesTopRightCorner_impl()
const
185 return m_nely * (m_nelx + 1) + m_nelx;
228 template <class C, class E, std::enable_if_t<xt::is_xexpression<C>::value,
bool> =
true>
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));
302 size_t n = m_layer_nelx(layer);
303 if (m_refine(layer) != -1) {
306 return m_startElem(layer) + xt::arange<size_t>(n);
315 std::vector<size_t> start_stop_rows,
316 std::vector<size_t> start_stop_cols)
const
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);
321 std::array<size_t, 2> rows;
322 std::array<size_t, 2> cols;
324 if (start_stop_rows.size() == 2) {
325 std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin());
331 rows[1] = this->
nely();
334 if (start_stop_cols.size() == 2) {
335 std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin());
341 cols[1] = this->
nelx();
344 if (rows[0] == rows[1] || cols[0] == cols[1]) {
351 auto H = xt::cumsum(m_nhy);
354 yl = xt::argmax(H > rows[0])();
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;
365 for (
size_t i = yl; i <= yu; ++i) {
367 size_t n = (xu - xl) * hx / m_nhx(i);
369 if (m_refine(i) != -1) {
381 for (
size_t i = yl; i <= yu; ++i) {
383 size_t n = (xu - xl) * hx / m_nhx(i);
386 if (m_refine(i) != -1) {
391 m_startElem(i) + xl *
h / m_nhx(i) + xt::arange<size_t>(n);
392 xt::view(ret, xt::range(N, N + n)) = e;
413 size_t iy = xt::argmin(m_startElem <= e)() - 1;
414 size_t nel = m_layer_nelx(iy);
418 auto H = xt::cumsum(m_nhy);
420 if (2 * size >= H(H.size() - 1)) {
421 return xt::arange<size_t>(this->
nelem());
425 size_t l = xt::argmax(H > (hy - size - 1))();
426 size_t u = xt::argmax(H >= (hy + size))();
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);
443 if (mid + size < nel) {
448 if (mid - xl < size) {
456 if (xu - mid < size) {
466 auto map = this->
roll(nroll);
467 return xt::view(map, xt::keep(ret));
486 size_t iy = xt::argmin(m_startElem <= e)() - 1;
487 size_t nel = m_layer_nelx(iy);
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);
501 if (mid + right < nel) {
502 xu = mid + right + 1;
506 if (mid - xl < left) {
514 if (xu - mid < right) {
523 auto H = xt::cumsum(m_nhy);
529 auto map = this->
roll(nroll);
530 return xt::view(map, xt::keep(ret));
541 size_t nely =
static_cast<size_t>(m_nhy.size());
545 for (
size_t iy = 0; iy <
nely; ++iy) {
548 size_t shift = n * (m_layer_nelx(iy) / m_layer_nelx(0));
549 size_t nel = m_layer_nelx(iy);
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;
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);
570 size_t nelx_impl()
const
572 return xt::amax(m_layer_nelx)();
575 size_t nely_impl()
const
577 return xt::sum(m_nhy)();
585 array_type::tensor<double, 2> coor_impl()
const
588 array_type::tensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
592 size_t nely =
static_cast<size_t>(m_nhy.size());
596 array_type::tensor<double, 1> y = xt::empty<double>({
nely + 1});
600 for (
size_t iy = 1; iy <
nely + 1; ++iy) {
601 y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h;
607 for (
size_t iy = 0;; ++iy) {
609 array_type::tensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
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);
619 if (iy == (
nely - 1) / 2) {
624 if (m_refine(iy) == 0) {
626 double dx = m_h *
static_cast<double>(m_nhx(iy) / 3);
627 double dy = m_h *
static_cast<double>(m_nhy(iy) / 2);
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;
641 for (
size_t iy = (
nely - 1) / 2; iy <
nely; ++iy) {
643 array_type::tensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
646 if (m_refine(iy) == 0) {
648 double dx = m_h *
static_cast<double>(m_nhx(iy) / 3);
649 double dy = m_h *
static_cast<double>(m_nhy(iy) / 2);
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;
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);
671 array_type::tensor<size_t, 2> conn_impl()
const
674 array_type::tensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
678 size_t nely =
static_cast<size_t>(m_nhy.size());
679 size_t bot, mid, top;
682 for (
size_t iy = 0; iy <
nely; ++iy) {
684 bot = m_startNode(iy);
685 mid = m_startNode(iy) + m_nnd(iy);
686 top = m_startNode(iy + 1);
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);
700 else if (m_refine(iy) == 0 && iy <= (
nely - 1) / 2) {
701 for (
size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
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);
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);
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);
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);
730 else if (m_refine(iy) == 0 && iy > (
nely - 1) / 2) {
731 for (
size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
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);
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);
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);
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);
763 array_type::tensor<size_t, 1> nodesBottomEdge_impl()
const
765 return m_startNode(0) + xt::arange<size_t>(m_layer_nelx(0) + 1);
768 array_type::tensor<size_t, 1> nodesTopEdge_impl()
const
770 size_t nely = m_nhy.size();
771 return m_startNode(
nely) + xt::arange<size_t>(m_layer_nelx(
nely - 1) + 1);
774 array_type::tensor<size_t, 1> nodesLeftEdge_impl()
const
776 size_t nely = m_nhy.size();
778 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({
nely + 1});
781 size_t j = (
nely + 1) / 2;
782 size_t k = (
nely - 1) / 2;
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));
791 array_type::tensor<size_t, 1> nodesRightEdge_impl()
const
793 size_t nely = m_nhy.size();
795 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({
nely + 1});
798 size_t j = (
nely + 1) / 2;
799 size_t k = (
nely - 1) / 2;
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));
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));
811 array_type::tensor<size_t, 1> nodesBottomOpenEdge_impl()
const
813 return m_startNode(0) + xt::arange<size_t>(1, m_layer_nelx(0));
816 array_type::tensor<size_t, 1> nodesTopOpenEdge_impl()
const
818 size_t nely = m_nhy.size();
820 return m_startNode(
nely) + xt::arange<size_t>(1, m_layer_nelx(
nely - 1));
823 array_type::tensor<size_t, 1> nodesLeftOpenEdge_impl()
const
825 size_t nely = m_nhy.size();
827 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({
nely - 1});
830 size_t j = (
nely + 1) / 2;
831 size_t k = (
nely - 1) / 2;
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));
840 array_type::tensor<size_t, 1> nodesRightOpenEdge_impl()
const
842 size_t nely = m_nhy.size();
844 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({
nely - 1});
847 size_t j = (
nely + 1) / 2;
848 size_t k = (
nely - 1) / 2;
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));
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));
860 size_t nodesBottomLeftCorner_impl()
const
862 return m_startNode(0);
865 size_t nodesBottomRightCorner_impl()
const
867 return m_startNode(0) + m_layer_nelx(0);
870 size_t nodesTopLeftCorner_impl()
const
872 size_t nely = m_nhy.size();
873 return m_startNode(
nely);
876 size_t nodesTopRightCorner_impl()
const
878 size_t nely = m_nhy.size();
879 return m_startNode(
nely) + m_layer_nelx(
nely - 1);
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;
899 void init(
size_t nelx,
size_t nely,
double h,
size_t nfine = 1)
907 m_Lx = m_h *
static_cast<double>(
nelx);
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});
922 nmin = (
nely + 1) / 2;
926 if (nfine % 2 == 0) {
927 nfine = nfine / 2 + 1;
930 nfine = (nfine + 1) / 2;
945 for (
size_t iy = nfine;;) {
951 if (iy >=
nely || ntot >= nmin) {
957 if (3 * nhy(iy) <= ntot &&
nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
960 auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
961 auto vnhx = xt::view(nhx, xt::range(iy, _));
971 if (iy >=
nely || ntot >= nmin) {
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});
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);
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);
1003 nely = m_nhx.size();
1006 for (
size_t iy = 0; iy <
nely; ++iy) {
1007 m_layer_nelx(iy) =
nelx / m_nhx(iy);
1011 for (
size_t iy = 0; iy < (
nely + 1) / 2; ++iy) {
1012 m_nnd(iy) = m_layer_nelx(iy) + 1;
1014 for (
size_t iy = (
nely - 1) / 2; iy <
nely; ++iy) {
1015 m_nnd(iy + 1) = m_layer_nelx(iy) + 1;
1026 for (
size_t i = 0; i < (
nely - 1) / 2; ++i) {
1028 m_startElem(i) = m_nelem;
1030 if (m_refine(i) == 0) {
1031 m_nnode += (3 * m_layer_nelx(i) + 1);
1034 m_nnode += (m_layer_nelx(i) + 1);
1037 if (m_refine(i) == 0) {
1038 m_nelem += (4 * m_layer_nelx(i));
1041 m_nelem += (m_layer_nelx(i));
1044 m_startNode(i + 1) = m_nnode;
1048 for (
size_t i = (
nely - 1) / 2; i <
nely; ++i) {
1050 m_startElem(i) = m_nelem;
1052 if (m_refine(i) == 0) {
1053 m_nnode += (5 * m_layer_nelx(i) + 1);
1056 m_nnode += (m_layer_nelx(i) + 1);
1059 if (m_refine(i) == 0) {
1060 m_nelem += (4 * m_layer_nelx(i));
1063 m_nelem += (m_layer_nelx(i));
1066 m_startNode(i + 1) = m_nnode;
1069 m_nnode += m_layer_nelx(
nely - 1) + 1;
1075 template <
class C,
class E>
1076 void init_by_mapping(
const C&
coor,
const E&
conn)
1086 if (
conn.shape(0) == 1) {
1095 size_t emid = (
conn.shape(0) -
conn.shape(0) % 2) / 2;
1096 size_t eleft = emid;
1097 size_t eright = emid;
1099 while (
conn(eleft, 0) ==
conn(eleft - 1, 1) && eleft > 0) {
1103 while (
conn(eright, 1) ==
conn(eright + 1, 0) && eright <
conn.shape(0) - 1) {
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)();
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);
1155 : m_coarse(mesh), m_nx(
nx), m_ny(
ny)
1162 m_coarse2fine = xt::empty<size_t>({m_coarse.
nelem(),
nx *
ny});
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)));
1216 return m_coarse2fine;
1243 return m_coarse2fine;
1254 template <
class T,
size_t rank>
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);
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);
1284 template <
class T,
size_t rank,
class S>
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);
1297 for (
size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
1298 auto e = xt::view(m_coarse2fine, i, xt::all());
1301 xt::view(ret, i) = xt::average(d, w, {0});
1319 template <
class T,
size_t rank>
1324 std::array<size_t, rank> shape;
1325 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1326 shape[0] = m_coarse2fine.size();
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);
1368 xt::amax(m_finelayer.m_layer_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);
1375 m_elem_regular.resize(m_finelayer.m_nelem);
1376 m_frac_regular.resize(m_finelayer.m_nelem);
1389 xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));
1392 size_t nely = nhy.size();
1395 for (
size_t iy = 0; iy < nely; ++iy) {
1397 auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());
1402 if (m_finelayer.m_refine(iy) == -1) {
1407 for (
size_t ix = 0; ix < nelx(iy); ++ix) {
1411 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
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);
1424 else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
1428 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
1431 for (
size_t ix = 0; ix < nelx(iy); ++ix) {
1435 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
1439 for (
size_t j = 0; j < nhy(iy) / 2; ++j) {
1440 auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));
1442 m_elem_regular[el_old(ix, 0)].push_back(e(0));
1443 m_frac_regular[el_old(ix, 0)].push_back(0.5);
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);
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);
1459 xt::range(nhy(iy) / 2, nhy(iy)),
1460 xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
1463 m_elem_regular[el_old(ix, 2)].push_back(i);
1464 m_frac_regular[el_old(ix, 2)].push_back(1.0);
1471 for (
size_t j = 0; j < nhy(iy) / 2; ++j) {
1472 auto e = xt::view(block, j, xt::range(0, j + 1));
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);
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);
1487 xt::range(nhy(iy) / 2, nhy(iy)),
1488 xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
1491 m_elem_regular[el_old(ix, 3)].push_back(i);
1492 m_frac_regular[el_old(ix, 3)].push_back(1.0);
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)));
1503 m_elem_regular[el_old(ix, 1)].push_back(e(0));
1504 m_frac_regular[el_old(ix, 1)].push_back(0.5);
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);
1516 xt::range(nhy(iy) / 2, nhy(iy)),
1517 xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
1520 m_elem_regular[el_old(ix, 1)].push_back(i);
1521 m_frac_regular[el_old(ix, 1)].push_back(1.0);
1529 else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
1533 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
1536 for (
size_t ix = 0; ix < nelx(iy); ++ix) {
1540 xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
1544 for (
size_t j = 0; j < nhy(iy) / 2; ++j) {
1548 xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));
1550 m_elem_regular[el_old(ix, 3)].push_back(e(0));
1551 m_frac_regular[el_old(ix, 3)].push_back(0.5);
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);
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);
1567 xt::range(0, nhy(iy) / 2),
1568 xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
1571 m_elem_regular[el_old(ix, 1)].push_back(i);
1572 m_frac_regular[el_old(ix, 1)].push_back(1.0);
1582 xt::range(0, nhy(iy) / 2),
1583 xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
1586 m_elem_regular[el_old(ix, 0)].push_back(i);
1587 m_frac_regular[el_old(ix, 0)].push_back(1.0);
1592 for (
size_t j = 0; j < nhy(iy) / 2; ++j) {
1594 xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));
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);
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);
1612 xt::range(0, nhy(iy) / 2),
1613 xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
1616 m_elem_regular[el_old(ix, 2)].push_back(i);
1617 m_frac_regular[el_old(ix, 2)].push_back(1.0);
1622 for (
size_t j = 0; j < nhy(iy) / 2; ++j) {
1624 block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));
1626 m_elem_regular[el_old(ix, 2)].push_back(e(0));
1627 m_frac_regular[el_old(ix, 2)].push_back(0.5);
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);
1669 std::vector<std::vector<size_t>>
map()
const
1671 return m_elem_regular;
1681 return m_frac_regular;
1713 [[deprecated]] std::vector<std::vector<size_t>>
getMap()
const
1715 return m_elem_regular;
1725 return m_frac_regular;
1741 template <
class T,
size_t rank>
1746 std::array<size_t, rank> shape;
1747 std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
1748 shape[0] = m_regular.
nelem();
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);
1764 std::vector<std::vector<size_t>> m_elem_regular;
1765 std::vector<std::vector<double>> m_frac_regular;
Mesh with fine middle layer, and coarser elements towards the top and bottom.
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.
array_type::tensor< size_t, 1 > roll(size_t n)
Mapping to 'roll' periodically in the x-direction,.
array_type::tensor< size_t, 1 > elementsLayer(size_t layer) const
Elements along a layer.
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...
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...
FineLayer(const C &coor, const E &conn)
Reconstruct class for given coordinates / connectivity.
FineLayer(size_t nelx, size_t nely, double h=1.0, size_t nfine=1)
Constructor.
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...
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.
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.
const array_type::tensor< size_t, 1 > & elemrow_nelem() const
Number of elements per row of blocks.
array_type::tensor< size_t, 1 > elementsMiddleLayer() const
Elements in the middle (fine) layer.
Map a FineLayer mesh to a Regular mesh.
array_type::tensor< T, rank > mapToRegular(const array_type::tensor< T, rank > &data) const
Map element quantities to Regular.
std::vector< std::vector< double > > getMapFraction() const
To overlap fraction for each item in the mapping in getMap().
GooseFEM::Mesh::Quad4::FineLayer fineLayerMesh() const
Obtain the FineLayer mesh (copy of the mesh passed to the constructor).
std::vector< std::vector< size_t > > map() const
Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh.
GooseFEM::Mesh::Quad4::Regular regularMesh() const
Obtain the Regular mesh.
GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const
Obtain the FineLayer mesh (copy of the mesh passed to the constructor).
GooseFEM::Mesh::Quad4::Regular getRegularMesh() const
Obtain the Regular mesh.
std::vector< std::vector< double > > mapFraction() const
To overlap fraction for each item in the mapping in map().
std::vector< std::vector< size_t > > getMap() const
Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh.
FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer &mesh)
Constructors.
Refine a Regular mesh: subdivide elements in several smaller elements.
GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const
Obtain the coarse mesh (copy of the mesh passed to the constructor).
const array_type::tensor< size_t, 2 > & getMap() const
Get element-mapping: elements of the fine mesh per element of the coarse mesh.
array_type::tensor< T, rank > mapToFine(const array_type::tensor< T, rank > &data) const
Map element quantities to the fine mesh.
GooseFEM::Mesh::Quad4::Regular fineMesh() const
Obtain the fine mesh.
GooseFEM::Mesh::Quad4::Regular coarseMesh() const
Obtain the coarse mesh (copy of the mesh passed to the constructor).
size_t ny() const
For each coarse element: number of fine elements in y-direction.
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.
RefineRegular(const GooseFEM::Mesh::Quad4::Regular &mesh, size_t nx, size_t ny)
Constructor.
size_t nx() const
For each coarse element: number of fine elements in x-direction.
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.
GooseFEM::Mesh::Quad4::Regular getFineMesh() const
Obtain the fine mesh.
const array_type::tensor< size_t, 2 > & map() const
Get element-mapping: elements of the fine mesh per element of the coarse mesh.
Regular mesh: equi-sized elements.
Regular(size_t nelx, size_t nely, double h=1.0)
Constructor.
array_type::tensor< size_t, 2 > elementgrid() const
Element numbers as 'matrix'.
CRTP base class for regular meshes in 2d.
CRTP base class for regular meshes.
auto nely() const
Number of elements in y-direction == height of the mesh, in units of h,.
auto coor() const
Nodal coordinates [nnode, ndim].
auto h() const
Linear edge size of one 'block'.
auto nelx() const
Number of elements in x-direction == width of the mesh in units of h.
auto nelem() const
Number of elements.
auto conn() const
Connectivity [nelem, nne].
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
#define GOOSEFEM_WIP_ASSERT(expr)
Assertion that concerns temporary implementation limitations.
#define GOOSEFEM_CHECK(expr)
Assertion that cannot be switched off.
ElementType
Enumerator for element-types.
@ Quad4
Quadrilateral: 4-noded element in 2-d.
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.