10#define GOOSEFEM_MESH_H
45template <
class S,
class T>
51 if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
54 if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
57 if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
61 throw std::runtime_error(
"Element-type not implemented");
66template <
class T,
class R>
67inline T renum(
const T& arg,
const R& mapping)
69 T ret = T::from_shape(arg.shape());
71 auto jt = ret.begin();
73 for (
auto it = arg.begin(); it != arg.end(); ++it, ++jt) {
95 return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
124 size_t n = xt::amax(
dofs)() + 1;
129 m_renum = xt::empty<size_t>({n});
131 for (
auto& j : unique) {
146 return detail::renum(list, m_renum);
194 return derived_cast().m_nelem;
203 return derived_cast().m_nnode;
212 return derived_cast().m_nne;
221 return derived_cast().m_ndim;
230 return derived_cast().nelx_impl();
239 return derived_cast().nely_impl();
248 return derived_cast().m_h;
257 return derived_cast().getElementType_impl();
266 return derived_cast().coor_impl();
275 return derived_cast().conn_impl();
298 for (
size_t j = 0; j < this->
ndim(); ++j) {
299 xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
311 return derived_cast().nodesPeriodic_impl();
320 return derived_cast().nodesOrigin_impl();
352 return derived_cast().nodesBottomEdge_impl();
361 return derived_cast().nodesTopEdge_impl();
370 return derived_cast().nodesLeftEdge_impl();
379 return derived_cast().nodesRightEdge_impl();
389 return derived_cast().nodesBottomOpenEdge_impl();
399 return derived_cast().nodesTopOpenEdge_impl();
409 return derived_cast().nodesLeftOpenEdge_impl();
419 return derived_cast().nodesRightOpenEdge_impl();
429 return derived_cast().nodesBottomLeftCorner_impl();
439 return derived_cast().nodesBottomRightCorner_impl();
449 return derived_cast().nodesTopLeftCorner_impl();
459 return derived_cast().nodesTopRightCorner_impl();
468 return derived_cast().nodesBottomLeftCorner_impl();
477 return derived_cast().nodesTopLeftCorner_impl();
486 return derived_cast().nodesBottomRightCorner_impl();
495 return derived_cast().nodesTopRightCorner_impl();
509 friend class RegularBase<D>;
511 array_type::tensor<size_t, 2> nodesPeriodic_impl()
const
513 array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomOpenEdge_impl();
514 array_type::tensor<size_t, 1> top = derived_cast().nodesTopOpenEdge_impl();
515 array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftOpenEdge_impl();
516 array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightOpenEdge_impl();
517 std::array<size_t, 2> shape = {bot.size() + lft.size() + size_t(3), size_t(2)};
518 auto ret = array_type::tensor<size_t, 2>::from_shape(shape);
520 ret(0, 0) = derived_cast().nodesBottomLeftCorner_impl();
521 ret(0, 1) = derived_cast().nodesBottomRightCorner_impl();
523 ret(1, 0) = derived_cast().nodesBottomLeftCorner_impl();
524 ret(1, 1) = derived_cast().nodesTopRightCorner_impl();
526 ret(2, 0) = derived_cast().nodesBottomLeftCorner_impl();
527 ret(2, 1) = derived_cast().nodesTopLeftCorner_impl();
531 xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
532 xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
536 xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
537 xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
542 auto nodesOrigin_impl()
const
544 return derived_cast().nodesBottomLeftCorner_impl();
565 return derived_cast().nelz_impl();
574 return derived_cast().nodesBottom_impl();
583 return derived_cast().nodesTop_impl();
592 return derived_cast().nodesLeft_impl();
601 return derived_cast().nodesRight_impl();
610 return derived_cast().nodesFront_impl();
619 return derived_cast().nodesBack_impl();
629 return derived_cast().nodesFrontBottomEdge_impl();
639 return derived_cast().nodesFrontTopEdge_impl();
649 return derived_cast().nodesFrontLeftEdge_impl();
659 return derived_cast().nodesFrontRightEdge_impl();
669 return derived_cast().nodesBackBottomEdge_impl();
679 return derived_cast().nodesBackTopEdge_impl();
689 return derived_cast().nodesBackLeftEdge_impl();
699 return derived_cast().nodesBackRightEdge_impl();
709 return derived_cast().nodesBottomLeftEdge_impl();
719 return derived_cast().nodesBottomRightEdge_impl();
729 return derived_cast().nodesTopLeftEdge_impl();
739 return derived_cast().nodesTopRightEdge_impl();
748 return derived_cast().nodesFrontBottomEdge_impl();
757 return derived_cast().nodesBackBottomEdge_impl();
766 return derived_cast().nodesFrontTopEdge_impl();
775 return derived_cast().nodesBackTopEdge_impl();
784 return derived_cast().nodesBottomLeftEdge_impl();
793 return derived_cast().nodesFrontLeftEdge_impl();
802 return derived_cast().nodesBackLeftEdge_impl();
811 return derived_cast().nodesTopLeftEdge_impl();
820 return derived_cast().nodesBottomRightEdge_impl();
829 return derived_cast().nodesTopRightEdge_impl();
838 return derived_cast().nodesFrontRightEdge_impl();
847 return derived_cast().nodesBackRightEdge_impl();
858 return derived_cast().nodesFrontFace_impl();
869 return derived_cast().nodesBackFace_impl();
880 return derived_cast().nodesLeftFace_impl();
891 return derived_cast().nodesRightFace_impl();
902 return derived_cast().nodesBottomFace_impl();
913 return derived_cast().nodesTopFace_impl();
922 return derived_cast().nodesFrontBottomOpenEdge_impl();
931 return derived_cast().nodesFrontTopOpenEdge_impl();
940 return derived_cast().nodesFrontLeftOpenEdge_impl();
949 return derived_cast().nodesFrontRightOpenEdge_impl();
958 return derived_cast().nodesBackBottomOpenEdge_impl();
967 return derived_cast().nodesBackTopOpenEdge_impl();
976 return derived_cast().nodesBackLeftOpenEdge_impl();
985 return derived_cast().nodesBackRightOpenEdge_impl();
994 return derived_cast().nodesBottomLeftOpenEdge_impl();
1003 return derived_cast().nodesBottomRightOpenEdge_impl();
1012 return derived_cast().nodesTopLeftOpenEdge_impl();
1021 return derived_cast().nodesTopRightOpenEdge_impl();
1030 return derived_cast().nodesFrontBottomOpenEdge_impl();
1039 return derived_cast().nodesBackBottomOpenEdge_impl();
1048 return derived_cast().nodesFrontTopOpenEdge_impl();
1057 return derived_cast().nodesBackTopOpenEdge_impl();
1066 return derived_cast().nodesBottomLeftOpenEdge_impl();
1075 return derived_cast().nodesFrontLeftOpenEdge_impl();
1084 return derived_cast().nodesBackLeftOpenEdge_impl();
1093 return derived_cast().nodesTopLeftOpenEdge_impl();
1102 return derived_cast().nodesBottomRightOpenEdge_impl();
1111 return derived_cast().nodesTopRightOpenEdge_impl();
1120 return derived_cast().nodesFrontRightOpenEdge_impl();
1129 return derived_cast().nodesBackRightOpenEdge_impl();
1138 return derived_cast().nodesFrontBottomLeftCorner_impl();
1147 return derived_cast().nodesFrontBottomRightCorner_impl();
1156 return derived_cast().nodesFrontTopLeftCorner_impl();
1165 return derived_cast().nodesFrontTopRightCorner_impl();
1174 return derived_cast().nodesBackBottomLeftCorner_impl();
1183 return derived_cast().nodesBackBottomRightCorner_impl();
1192 return derived_cast().nodesBackTopLeftCorner_impl();
1201 return derived_cast().nodesBackTopRightCorner_impl();
1210 return derived_cast().nodesFrontBottomLeftCorner_impl();
1219 return derived_cast().nodesFrontBottomLeftCorner_impl();
1228 return derived_cast().nodesFrontBottomLeftCorner_impl();
1237 return derived_cast().nodesFrontBottomLeftCorner_impl();
1246 return derived_cast().nodesFrontBottomLeftCorner_impl();
1255 return derived_cast().nodesFrontBottomRightCorner_impl();
1264 return derived_cast().nodesFrontBottomRightCorner_impl();
1273 return derived_cast().nodesFrontBottomRightCorner_impl();
1282 return derived_cast().nodesFrontBottomRightCorner_impl();
1291 return derived_cast().nodesFrontBottomRightCorner_impl();
1300 return derived_cast().nodesFrontTopLeftCorner_impl();
1309 return derived_cast().nodesFrontTopLeftCorner_impl();
1318 return derived_cast().nodesFrontTopLeftCorner_impl();
1327 return derived_cast().nodesFrontTopLeftCorner_impl();
1336 return derived_cast().nodesFrontTopLeftCorner_impl();
1345 return derived_cast().nodesFrontTopRightCorner_impl();
1354 return derived_cast().nodesFrontTopRightCorner_impl();
1363 return derived_cast().nodesFrontTopRightCorner_impl();
1372 return derived_cast().nodesFrontTopRightCorner_impl();
1381 return derived_cast().nodesFrontTopRightCorner_impl();
1390 return derived_cast().nodesBackBottomLeftCorner_impl();
1399 return derived_cast().nodesBackBottomLeftCorner_impl();
1408 return derived_cast().nodesBackBottomLeftCorner_impl();
1417 return derived_cast().nodesBackBottomLeftCorner_impl();
1426 return derived_cast().nodesBackBottomLeftCorner_impl();
1435 return derived_cast().nodesBackBottomRightCorner_impl();
1444 return derived_cast().nodesBackBottomRightCorner_impl();
1453 return derived_cast().nodesBackBottomRightCorner_impl();
1462 return derived_cast().nodesBackBottomRightCorner_impl();
1471 return derived_cast().nodesBackBottomRightCorner_impl();
1480 return derived_cast().nodesBackTopLeftCorner_impl();
1489 return derived_cast().nodesBackTopLeftCorner_impl();
1498 return derived_cast().nodesBackTopLeftCorner_impl();
1507 return derived_cast().nodesBackTopLeftCorner_impl();
1516 return derived_cast().nodesBackTopLeftCorner_impl();
1525 return derived_cast().nodesBackTopRightCorner_impl();
1534 return derived_cast().nodesBackTopRightCorner_impl();
1543 return derived_cast().nodesBackTopRightCorner_impl();
1552 return derived_cast().nodesBackTopRightCorner_impl();
1561 return derived_cast().nodesBackTopRightCorner_impl();
1575 friend class RegularBase<D>;
1577 array_type::tensor<size_t, 2> nodesPeriodic_impl()
const
1579 array_type::tensor<size_t, 1> fro = derived_cast().nodesFrontFace_impl();
1580 array_type::tensor<size_t, 1> bck = derived_cast().nodesBackFace_impl();
1581 array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftFace_impl();
1582 array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightFace_impl();
1583 array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomFace_impl();
1584 array_type::tensor<size_t, 1> top = derived_cast().nodesTopFace_impl();
1586 array_type::tensor<size_t, 1> froBot = derived_cast().nodesFrontBottomOpenEdge_impl();
1587 array_type::tensor<size_t, 1> froTop = derived_cast().nodesFrontTopOpenEdge_impl();
1588 array_type::tensor<size_t, 1> froLft = derived_cast().nodesFrontLeftOpenEdge_impl();
1589 array_type::tensor<size_t, 1> froRgt = derived_cast().nodesFrontRightOpenEdge_impl();
1590 array_type::tensor<size_t, 1> bckBot = derived_cast().nodesBackBottomOpenEdge_impl();
1591 array_type::tensor<size_t, 1> bckTop = derived_cast().nodesBackTopOpenEdge_impl();
1592 array_type::tensor<size_t, 1> bckLft = derived_cast().nodesBackLeftOpenEdge_impl();
1593 array_type::tensor<size_t, 1> bckRgt = derived_cast().nodesBackRightOpenEdge_impl();
1594 array_type::tensor<size_t, 1> botLft = derived_cast().nodesBottomLeftOpenEdge_impl();
1595 array_type::tensor<size_t, 1> botRgt = derived_cast().nodesBottomRightOpenEdge_impl();
1596 array_type::tensor<size_t, 1> topLft = derived_cast().nodesTopLeftOpenEdge_impl();
1597 array_type::tensor<size_t, 1> topRgt = derived_cast().nodesTopRightOpenEdge_impl();
1599 size_t tface = fro.size() + lft.size() + bot.size();
1600 size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size();
1602 array_type::tensor<size_t, 2> ret =
1603 xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
1607 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1608 ret(i, 1) = derived_cast().nodesFrontBottomRightCorner_impl();
1610 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1611 ret(i, 1) = derived_cast().nodesBackBottomRightCorner_impl();
1613 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1614 ret(i, 1) = derived_cast().nodesBackBottomLeftCorner_impl();
1616 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1617 ret(i, 1) = derived_cast().nodesFrontTopLeftCorner_impl();
1619 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1620 ret(i, 1) = derived_cast().nodesFrontTopRightCorner_impl();
1622 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1623 ret(i, 1) = derived_cast().nodesBackTopRightCorner_impl();
1625 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1626 ret(i, 1) = derived_cast().nodesBackTopLeftCorner_impl();
1629 for (
size_t j = 0; j < froBot.size(); ++j) {
1630 ret(i, 0) = froBot(j);
1631 ret(i, 1) = bckBot(j);
1634 for (
size_t j = 0; j < froBot.size(); ++j) {
1635 ret(i, 0) = froBot(j);
1636 ret(i, 1) = bckTop(j);
1639 for (
size_t j = 0; j < froBot.size(); ++j) {
1640 ret(i, 0) = froBot(j);
1641 ret(i, 1) = froTop(j);
1644 for (
size_t j = 0; j < botLft.size(); ++j) {
1645 ret(i, 0) = botLft(j);
1646 ret(i, 1) = botRgt(j);
1649 for (
size_t j = 0; j < botLft.size(); ++j) {
1650 ret(i, 0) = botLft(j);
1651 ret(i, 1) = topRgt(j);
1654 for (
size_t j = 0; j < botLft.size(); ++j) {
1655 ret(i, 0) = botLft(j);
1656 ret(i, 1) = topLft(j);
1659 for (
size_t j = 0; j < froLft.size(); ++j) {
1660 ret(i, 0) = froLft(j);
1661 ret(i, 1) = froRgt(j);
1664 for (
size_t j = 0; j < froLft.size(); ++j) {
1665 ret(i, 0) = froLft(j);
1666 ret(i, 1) = bckRgt(j);
1669 for (
size_t j = 0; j < froLft.size(); ++j) {
1670 ret(i, 0) = froLft(j);
1671 ret(i, 1) = bckLft(j);
1675 for (
size_t j = 0; j < fro.size(); ++j) {
1680 for (
size_t j = 0; j < lft.size(); ++j) {
1685 for (
size_t j = 0; j < bot.size(); ++j) {
1694 auto nodesOrigin_impl()
const
1696 return derived_cast().nodesFrontBottomLeftCorner_impl();
1714template <
class S,
class T>
1715inline array_type::tensor<size_t, 2>
1716overlapping(
const S& coor_a,
const T& coor_b,
double rtol = 1e-5,
double atol = 1e-8)
1722 std::vector<size_t> ret_a;
1723 std::vector<size_t> ret_b;
1725 for (
size_t i = 0; i < coor_a.shape(0); ++i) {
1727 auto idx = xt::flatten_indices(xt::argwhere(
1728 xt::prod(xt::isclose(coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1)));
1730 for (
auto& j : idx) {
1737 for (
size_t i = 0; i < ret_a.size(); ++i) {
1738 ret(0, i) = ret_a[i];
1739 ret(1, i) = ret_b[i];
1763 template <
class CA,
class EA,
class NA,
class CB,
class EB,
class NB>
1767 const NA& overlapping_nodes_a,
1770 const NB& overlapping_nodes_b,
1771 bool check_position =
true,
1784 GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape()));
1788 if (check_position) {
1790 xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()),
1791 xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()),
1796 size_t nnda = coor_a.shape(0);
1797 size_t nndb = coor_b.shape(0);
1798 size_t ndim = coor_a.shape(1);
1799 size_t nelim = overlapping_nodes_a.size();
1801 size_t nela = conn_a.shape(0);
1802 size_t nelb = conn_b.shape(0);
1803 size_t nne = conn_a.shape(1);
1810 xt::setdiff1d(xt::arange<size_t>(nndb), overlapping_nodes_b);
1812 m_map_b = xt::empty<size_t>({nndb});
1813 xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a;
1814 xt::view(m_map_b, xt::keep(keep_b)) = xt::arange<size_t>(keep_b.size()) + nnda;
1816 m_conn = xt::empty<size_t>({nela + nelb,
nne});
1817 xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a;
1818 xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b);
1820 m_coor = xt::empty<size_t>({nnda + nndb - nelim,
ndim});
1821 xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a;
1822 xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) =
1823 xt::view(coor_b, xt::keep(keep_b), xt::all());
1841 return m_conn.shape(0);
1850 return m_coor.shape(0);
1859 return m_conn.shape(1);
1868 return m_coor.shape(1);
1904 std::vector<array_type::tensor<size_t, 1>>
nodemap()
const
1906 std::vector<array_type::tensor<size_t, 1>> ret(this->
nmesh());
1907 for (
size_t i = 0; i < this->
nmesh(); ++i) {
1917 std::vector<array_type::tensor<size_t, 1>>
elemmap()
const
1919 std::vector<array_type::tensor<size_t, 1>> ret(this->
nmesh());
1920 for (
size_t i = 0; i < this->
nmesh(); ++i) {
1934 if (mesh_index == 0) {
1935 return xt::arange<size_t>(m_nnd_a);
1949 if (mesh_index == 0) {
1950 return xt::arange<size_t>(m_nel_a);
1953 return xt::arange<size_t>(m_nel_b) + m_nel_a;
1968 if (mesh_index == 0) {
1974 return detail::renum(set, m_map_b);
1989 if (mesh_index == 0) {
1995 return set + m_nel_a;
2016 Stitch(
double rtol = 1e-5,
double atol = 1e-8)
2028 template <
class C,
class E>
2034 if (
m_map.size() == 0) {
2037 m_map.push_back(xt::eval(xt::arange<size_t>(
coor.shape(0))));
2044 size_t index =
m_map.size();
2049 xt::eval(xt::view(overlap, 0, xt::all())),
2052 xt::eval(xt::view(overlap, 1, xt::all())),
2068 return m_map.size();
2140 std::vector<array_type::tensor<size_t, 1>>
nodemap()
const
2142 std::vector<array_type::tensor<size_t, 1>> ret(this->
nmesh());
2143 for (
size_t i = 0; i < this->
nmesh(); ++i) {
2153 std::vector<array_type::tensor<size_t, 1>>
elemmap()
const
2155 std::vector<array_type::tensor<size_t, 1>> ret(this->
nmesh());
2156 for (
size_t i = 0; i < this->
nmesh(); ++i) {
2171 return m_map[mesh_index];
2198 return detail::renum(set,
m_map[mesh_index]);
2229 for (
size_t i = 0; i < set.size(); ++i) {
2237 for (
size_t i = 0; i < set.size(); ++i) {
2238 xt::view(ret, xt::range(n, n + set[i].size())) = this->
nodeset(set[i], i);
2242 return xt::unique(ret);
2251 return this->
nodeset(std::vector<T>(set));
2267 for (
size_t i = 0; i < set.size(); ++i) {
2275 for (
size_t i = 0; i < set.size(); ++i) {
2276 xt::view(ret, xt::range(n, n + set[i].size())) = this->
elemset(set[i], i);
2289 return this->
elemset(std::vector<T>(set));
2295 std::vector<array_type::tensor<size_t, 1>>
m_map;
2312 Vstack(
bool check_overlap =
true,
double rtol = 1e-5,
double atol = 1e-8)
2314 m_check_overlap = check_overlap;
2328 template <
class C,
class E,
class N>
2331 if (
m_map.size() == 0) {
2334 m_map.push_back(xt::eval(xt::arange<size_t>(
coor.shape(0))));
2337 m_nodes_bot.push_back(nodes_bot);
2338 m_nodes_top.push_back(nodes_top);
2344 size_t index =
m_map.size();
2346 double shift = xt::amax(xt::view(
m_coor, xt::all(), 1))();
2348 xt::view(x, xt::all(), 1) += shift;
2361 m_nodes_bot.push_back(stitch.
nodeset(nodes_bot, 1));
2362 m_nodes_top.push_back(stitch.
nodeset(nodes_top, 1));
2372 std::vector<array_type::tensor<size_t, 1>>
2374 std::vector<array_type::tensor<size_t, 1>>
2376 bool m_check_overlap;
2400 for (
auto& arg : args) {
2401 if (arg.size() == 0) {
2404 n = std::max(n, xt::amax(arg)() + 1);
2407#ifdef GOOSEFEM_ENABLE_ASSERT
2408 for (
auto& arg : args) {
2413 m_renum = xt::empty<size_t>({n});
2415 for (
auto& arg : args) {
2416 for (
auto& j : arg) {
2432 T ret = T::from_shape(list.shape());
2434 auto jt = ret.begin();
2436 for (
auto it = list.begin(); it != list.end(); ++it, ++jt) {
2470 size_t nnode = xt::amax(conn)() + 1;
2474 for (
auto it = conn.begin(); it != conn.end(); ++it) {
2489inline std::vector<std::vector<size_t>>
elem2node(
const E& conn,
bool sorted =
true)
2492 auto nnode = N.size();
2494 std::vector<std::vector<size_t>> ret(nnode);
2495 for (
size_t i = 0; i < nnode; ++i) {
2496 ret[i].reserve(N(i));
2499 for (
size_t e = 0; e < conn.shape(0); ++e) {
2500 for (
size_t m = 0; m < conn.shape(1); ++m) {
2501 ret[conn(e, m)].push_back(e);
2506 for (
auto& row : ret) {
2507 std::sort(row.begin(), row.end());
2519template <
class E,
class D>
2520inline std::vector<std::vector<size_t>>
elem2node(
const E& conn,
const D&
dofs,
bool sorted =
true)
2522 size_t nnode =
dofs.shape(0);
2526 for (
size_t m = 0; m < nnode; ++m) {
2527 if (nties[m].size() <= 1) {
2530 if (m > nties[m][0]) {
2533 size_t n = ret[m].size();
2534 for (
size_t j = 1; j < nties[m].size(); ++j) {
2535 n += ret[nties[m][j]].size();
2538 for (
size_t j = 1; j < nties[m].size(); ++j) {
2539 ret[m].insert(ret[m].end(), ret[nties[m][j]].begin(), ret[nties[m][j]].end());
2542 std::sort(ret[m].begin(), ret[m].end());
2544 for (
size_t j = 1; j < nties[m].size(); ++j) {
2545 ret[nties[m][j]] = ret[m];
2560inline std::vector<std::vector<size_t>>
node2dof(
const D&
dofs,
bool sorted =
true)
2573template <
class C,
class E>
2596 xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
2597 xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
2598 xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
2599 xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
2603 throw std::runtime_error(
"Element-type not implemented");
2614template <
class C,
class E>
2628template <
class C,
class E>
2639 for (
size_t i = 0; i < 4; ++i) {
2640 auto n = xt::view(conn, xt::all(), i);
2641 ret += xt::view(coor, xt::keep(n), xt::all());
2647 throw std::runtime_error(
"Element-type not implemented");
2658template <
class C,
class E>
2673template <
class T,
class C,
class E>
2682 size_t N = coor.shape(0);
2690 for (
size_t i = 0; i < 4; ++i) {
2692 auto old_nd = xt::view(conn, xt::all(), i);
2693 auto new_nd = xt::view(conn, xt::keep(elem_map), i);
2694 xt::view(t, xt::keep(old_nd)) = new_nd;
2695 ret = xt::where(xt::equal(ret, N), t, ret);
2701 throw std::runtime_error(
"Element-type not implemented");
2713template <
class T,
class C,
class E>
2734template <
class C,
class E>
2737 auto dof =
dofs(coor.shape(0), coor.shape(1));
2750 throw std::runtime_error(
"Element-type not implemented");
2769template <
class C,
class E>
2786 for (
size_t j = 0; j < weights.shape(1); ++j) {
2788 for (
size_t i = 0; i < weights.shape(0); ++i) {
2789 ret(j) += data(i, j) * weights(i, j);
2790 norm += weights(i, j);
2809template <
class C,
class E>
2810inline array_type::tensor<double, 1>
2814 return detail::average_axis_0(coor,
nodal_mass(coor, conn, type));
2827template <
class C,
class E>
2844 size_t nnode =
dofs.shape(0);
2845 size_t ndim =
dofs.shape(1);
2847 std::vector<std::vector<size_t>> ret(nnode);
2849 for (
size_t m = 0; m < nnode; ++m) {
2850 auto r = nodemap[
dofs(m, 0)];
2851 std::sort(r.begin(), r.end());
2853#ifdef GOOSEFEM_ENABLE_ASSERT
2854 for (
size_t i = 1; i < ndim; ++i) {
2855 auto t = nodemap[
dofs(m, i)];
2856 std::sort(t.begin(), t.end());
Quadrature for 4-noded quadrilateral element in 2d (GooseFEM::Mesh::ElementType::Quad4),...
Methods to switch between storage types based on a mesh.
Interpolation and quadrature.
auto Int_N_scalar_NT_dV(const T &qscalar) const -> array_type::tensor< double, 3 >
Element-by-element: integral of the scalar product of the shape function with a scalar.
void assemble(const T &elemmat)
Assemble from "elemmat".
const array_type::tensor< double, 1 > & data() const
Underlying matrix.
Stitch two mesh objects, specifying overlapping nodes by hand.
T elemset(const T &set, size_t mesh_index) const
Convert set of element numbers for an original mesh to the stitched mesh.
array_type::tensor< size_t, 1 > nodemap(size_t mesh_index) const
std::vector< array_type::tensor< size_t, 1 > > elemmap() const
Element-map per sub-mesh.
const array_type::tensor< size_t, 2 > & conn() const
Connectivity [nelem, nne].
array_type::tensor< size_t, 2 > dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
T nodeset(const T &set, size_t mesh_index) const
Convert set of node numbers for an original mesh to the stitched mesh.
size_t nne() const
Number of nodes-per-element.
size_t nmesh() const
Number of sub meshes == 2.
array_type::tensor< size_t, 1 > elemmap(size_t mesh_index) const
std::vector< array_type::tensor< size_t, 1 > > nodemap() const
Node-map per sub-mesh.
const array_type::tensor< double, 2 > & coor() const
Nodal coordinates [nnode, ndim].
size_t nnode() const
Number of nodes.
ManualStitch(const CA &coor_a, const EA &conn_a, const NA &overlapping_nodes_a, const CB &coor_b, const EB &conn_b, const NB &overlapping_nodes_b, bool check_position=true, double rtol=1e-5, double atol=1e-8)
size_t nelem() const
Number of elements.
size_t ndim() const
Number of dimensions.
CRTP base class for regular meshes in 2d.
auto nodesTopRightCorner() const
The top-right corner node (at x = nelx * h, y = nely * h).
auto nodesTopLeftCorner() const
The top-left corner node (at x = 0, y = nely * h).
auto nodesLeftEdge() const
Nodes along the left edge (x = 0), in order of increasing y.
auto nodesBottomEdge() const
Nodes along the bottom edge (y = 0), in order of increasing x.
auto nodesRightOpenEdge() const
Nodes along the right edge (x = nelx * h), without the corners (at y = 0 and y = nely * h).
auto nodesRightBottomCorner() const
Alias of nodesBottomRightCorner().
auto nodesRightEdge() const
Nodes along the right edge (x = nelx * h), in order of increasing y.
auto nodesBottomLeftCorner() const
The bottom-left corner node (at x = 0, y = 0).
auto nodesTopEdge() const
Nodes along the top edge (y = nely * h), in order of increasing x.
auto nodesBottomRightCorner() const
The bottom-right corner node (at x = nelx * h, y = 0).
auto nodesRightTopCorner() const
Alias of nodesTopRightCorner().
auto nodesLeftBottomCorner() const
Alias of nodesBottomLeftCorner().
auto nodesLeftOpenEdge() const
Nodes along the left edge (x = 0), without the corners (at y = 0 and y = nely * h).
auto nodesBottomOpenEdge() const
Nodes along the bottom edge (y = 0), without the corners (at x = 0 and x = nelx * h).
auto nodesLeftTopCorner() const
Alias of nodesTopLeftCorner().
D derived_type
Underlying type.
auto nodesTopOpenEdge() const
Nodes along the top edge (y = nely * h), without the corners (at x = 0 and x = nelx * h).
CRTP base class for regular meshes in 3d.
auto nodesRightBackOpenEdge() const
Alias of nodesBackRightOpenEdge().
auto nodesFrontRightTopCorner() const
Alias of nodesFrontTopRightCorner().
auto nodesRightBottomFrontCorner() const
Alias of nodesFrontBottomRightCorner().
auto nodesFrontRightOpenEdge() const
Same as nodesFrontRightEdge() but without corners.
auto nodesFrontTopLeftCorner() const
Front-Top-Left corner node.
auto nodesFrontBottomRightCorner() const
Front-Bottom-Right corner node.
auto nodesFrontLeftTopCorner() const
Alias of nodesFrontTopLeftCorner().
auto nodesTopBackRightCorner() const
Alias of nodesBackTopRightCorner().
auto nodesFrontTopEdge() const
Nodes along the edge at the intersection of the front and top faces (z = 0 and y = nely * h).
auto nodesFrontLeftOpenEdge() const
Same as nodesFrontLeftEdge() but without corners.
auto nodesTopLeftEdge() const
Nodes along the edge at the intersection of the top and left faces (y = 0 and x = nelx * h).
auto nodesBackLeftEdge() const
Nodes along the edge at the intersection of the back and left faces (z = nelz * h and x = nelx * h).
auto nodesLeftBottomEdge() const
Alias of nodesBottomLeftEdge()
auto nodesFrontLeftBottomCorner() const
Alias of nodesFrontBottomLeftCorner().
auto nodesRightTopEdge() const
Alias of nodesTopRightEdge()
auto nodesBackRightBottomCorner() const
Alias of nodesBackBottomRightCorner().
auto nodesBottomRightOpenEdge() const
Same as nodesBottomRightEdge() but without corners.
auto nodesFrontBottomEdge() const
Nodes along the edge at the intersection of the front and bottom faces (z = 0 and y = 0).
auto nodesTopRightOpenEdge() const
Same as nodesTopRightEdge() but without corners.
auto nodesFrontBottomOpenEdge() const
Same as nodesFrontBottomEdge() but without corners.
auto nodesBack() const
Nodes along the back face (z = nelz * h).
auto nodesBottomRightFrontCorner() const
Alias of nodesFrontBottomRightCorner().
auto nodesTopRightBackCorner() const
Alias of nodesBackTopRightCorner().
auto nelz() const
Number of elements in y-direction == height of the mesh, in units of h,.
auto nodesRightBottomOpenEdge() const
Alias of nodesBottomRightOpenEdge().
auto nodesLeftBottomFrontCorner() const
Alias of nodesFrontBottomLeftCorner().
auto nodesLeftTopEdge() const
Alias of nodesTopLeftEdge()
auto nodesBackRightTopCorner() const
Alias of nodesBackTopRightCorner().
auto nodesTopLeftFrontCorner() const
Alias of nodesFrontTopLeftCorner().
auto nodesLeftFrontEdge() const
Alias of nodesFrontLeftEdge()
auto nodesFrontTopOpenEdge() const
Same as nodesFrontTopEdge() but without corners.
auto nodesTopLeftOpenEdge() const
Same as nodesTopLeftEdge() but without corners.
auto nodesTopLeftBackCorner() const
Alias of nodesBackTopLeftCorner().
auto nodesLeftFrontBottomCorner() const
Alias of nodesFrontBottomLeftCorner().
auto nodesTopFrontLeftCorner() const
Alias of nodesFrontTopLeftCorner().
auto nodesLeftFrontOpenEdge() const
Alias of nodesFrontLeftOpenEdge().
auto nodesBottomBackLeftCorner() const
Alias of nodesBackBottomLeftCorner().
auto nodesTopRightFrontCorner() const
Alias of nodesFrontTopRightCorner().
auto nodesBackLeftBottomCorner() const
Alias of nodesBackBottomLeftCorner().
auto nodesBottomLeftEdge() const
Nodes along the edge at the intersection of the bottom and left faces (y = 0 and x = 0).
auto nodesTopBackOpenEdge() const
Alias of nodesBackTopOpenEdge().
auto nodesBottomBackOpenEdge() const
Alias of nodesBackBottomOpenEdge().
auto nodesBottomLeftOpenEdge() const
Same as nodesBottomLeftEdge() but without corners.
auto nodesBackTopOpenEdge() const
Same as nodesBackTopEdge() but without corners.
auto nodesLeftFace() const
Nodes along the left face excluding edges.
auto nodesRightBottomEdge() const
Alias of nodesBottomRightEdge()
auto nodesBackRightEdge() const
Nodes along the edge at the intersection of the back and right faces (? = nelz * h and ?...
auto nodesRightTopBackCorner() const
Alias of nodesBackTopRightCorner().
auto nodesBottomLeftBackCorner() const
Alias of nodesBackBottomLeftCorner().
auto nodesTopFace() const
Nodes along the top face excluding edges.
auto nodesLeftBackTopCorner() const
Alias of nodesBackTopLeftCorner().
auto nodesLeft() const
Nodes along the left face (x = 0).
auto nodesTopBackLeftCorner() const
Alias of nodesBackTopLeftCorner().
auto nodesBackTopLeftCorner() const
Back-Top-Left corner node.
auto nodesTopFrontOpenEdge() const
Alias of nodesFrontTopOpenEdge().
auto nodesRightTopFrontCorner() const
Alias of nodesFrontTopRightCorner().
auto nodesTopFrontEdge() const
Alias of nodesFrontTopEdge()
auto nodesRightFace() const
Nodes along the right face excluding edges.
auto nodesLeftBottomBackCorner() const
Alias of nodesBackBottomLeftCorner().
auto nodesBottomRightEdge() const
Nodes along the edge at the intersection of the bottom and right faces (y = 0 and x = nelx * h).
auto nodesFrontRightEdge() const
Nodes along the edge at the intersection of the front and right faces (z = 0 and x = nelx * h).
auto nodesTopRightEdge() const
Nodes along the edge at the intersection of the top and right faces (y = nely * h and x = nelx * h).
auto nodesBackTopEdge() const
Nodes along the edge at the intersection of the back and top faces (z = nelz * h and x = 0).
auto nodesLeftBottomOpenEdge() const
Alias of nodesBottomLeftOpenEdge().
auto nodesRightFrontEdge() const
Alias of nodesFrontRightEdge()
auto nodesBackBottomLeftCorner() const
Back-Bottom-Left corner node.
auto nodesFrontFace() const
Nodes along the front face excluding edges.
auto nodesBackLeftTopCorner() const
Alias of nodesBackTopLeftCorner().
auto nodesFrontTopRightCorner() const
Front-Top-Right corner node.
auto nodesRightBackBottomCorner() const
Alias of nodesBackBottomRightCorner().
auto nodesLeftFrontTopCorner() const
Alias of nodesFrontTopLeftCorner().
auto nodesRightBackTopCorner() const
Alias of nodesBackTopRightCorner().
auto nodesLeftTopFrontCorner() const
Alias of nodesFrontTopLeftCorner().
auto nodesRightTopOpenEdge() const
Alias of nodesTopRightOpenEdge().
auto nodesBackRightOpenEdge() const
Same as nodesBackRightEdge() but without corners.
auto nodesBackBottomRightCorner() const
Back-Bottom-Right corner node.
auto nodesBottomFrontRightCorner() const
Alias of nodesFrontBottomRightCorner().
auto nodesBottomRightBackCorner() const
Alias of nodesBackBottomRightCorner().
auto nodesLeftTopOpenEdge() const
Alias of nodesTopLeftOpenEdge().
auto nodesBottomFrontLeftCorner() const
Alias of nodesFrontBottomLeftCorner().
D derived_type
Underlying type.
auto nodesBottomFrontEdge() const
Alias of nodesFrontBottomEdge()
auto nodesBottomLeftFrontCorner() const
Alias of nodesFrontBottomLeftCorner().
auto nodesFront() const
Nodes along the front face (z = 0).
auto nodesBackBottomEdge() const
Nodes along the edge at the intersection of the back and bottom faces (z = nelz * h and y = nely * h)...
auto nodesBottomFrontOpenEdge() const
Alias of nodesFrontBottomOpenEdge().
auto nodesBackBottomOpenEdge() const
Same as nodesBackBottomEdge() but without corners.
auto nodesBackFace() const
Nodes along the back face excluding edges.
auto nodesBottomBackEdge() const
Alias of nodesBackBottomEdge()
auto nodesBackLeftOpenEdge() const
Same as nodesBackLeftEdge() but without corners.
auto nodesBottomBackRightCorner() const
Alias of nodesBackBottomRightCorner().
auto nodesRightBackEdge() const
Alias of nodesBackRightEdge()
auto nodesBottomFace() const
Nodes along the bottom face excluding edges.
auto nodesLeftBackOpenEdge() const
Alias of nodesBackLeftOpenEdge().
auto nodesFrontLeftEdge() const
Nodes along the edge at the intersection of the front and left faces (z = 0 and x = 0).
auto nodesTopFrontRightCorner() const
Alias of nodesFrontTopRightCorner().
auto nodesFrontBottomLeftCorner() const
Front-Bottom-Left corner node.
auto nodesRight() const
Nodes along the right face (x = nelx * h).
auto nodesLeftTopBackCorner() const
Alias of nodesBackTopLeftCorner().
auto nodesRightBottomBackCorner() const
Alias of nodesBackBottomRightCorner().
auto nodesBottom() const
Nodes along the bottom face (y = 0).
auto nodesFrontRightBottomCorner() const
Alias of nodesFrontBottomRightCorner().
auto nodesLeftBackEdge() const
Alias of nodesBackLeftEdge()
auto nodesLeftBackBottomCorner() const
Alias of nodesBackBottomLeftCorner().
auto nodesRightFrontOpenEdge() const
Alias of nodesFrontRightOpenEdge().
auto nodesTopBackEdge() const
Alias of nodesBackTopEdge()
auto nodesBackTopRightCorner() const
Back-Top-Right corner node.
auto nodesRightFrontBottomCorner() const
Alias of nodesFrontBottomRightCorner().
auto nodesRightFrontTopCorner() const
Alias of nodesFrontTopRightCorner().
auto nodesTop() const
Nodes along the top face (y = nely * h).
CRTP base class for regular meshes.
auto ndim() const
Number of dimensions == 2.
D derived_type
Underlying type.
auto nely() const
Number of elements in y-direction == height of the mesh, in units of h,.
auto nnode() const
Number of nodes.
auto nodesPeriodic() const
Periodic node pairs, in two columns: (independent, dependent).
auto dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
auto coor() const
Nodal coordinates [nnode, ndim].
auto dofsPeriodic() const
DOF-numbers for the case that the periodicity if fully eliminated.
auto getElementType() const
The ElementType().
auto nodesOrigin() const
Reference node to use for periodicity, because all corners are tied to it.
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 nne() const
Number of nodes-per-element == 4.
auto nelem() const
Number of elements.
auto conn() const
Connectivity [nelem, nne].
Renumber indices to lowest possible index.
T apply(const T &list) const
Apply renumbering to other set.
const array_type::tensor< size_t, 1 > & index() const
Get the list needed to renumber, e.g.:
Reorder to lowest possible index, in specific order.
T apply(const T &list) const
Apply reordering to other set.
const array_type::tensor< size_t, 1 > & index() const
Get the list needed to reorder, e.g.:
Reorder(const std::initializer_list< T > args)
Stitch mesh objects, automatically searching for overlapping nodes.
T nodeset(const std::vector< T > &set) const
Combine set of node numbers for an original to the final mesh (removes duplicates).
size_t ndim() const
Number of dimensions.
T nodeset(std::initializer_list< T > set) const
Combine set of node numbers for an original to the final mesh (removes duplicates).
std::vector< size_t > m_nel
Number of elements per sub-mesh.
std::vector< array_type::tensor< size_t, 1 > > m_map
See nodemap(size_t)
size_t nne() const
Number of nodes-per-element.
const array_type::tensor< size_t, 2 > & conn() const
Connectivity [nelem, nne].
array_type::tensor< size_t, 2 > dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
size_t nelem() const
Number of elements.
double m_atol
Absolute tolerance to find overlapping nodes.
const array_type::tensor< double, 2 > & coor() const
Nodal coordinates [nnode, ndim].
T elemset(const std::vector< T > &set) const
Combine set of element numbers for an original to the final mesh.
std::vector< array_type::tensor< size_t, 1 > > nodemap() const
Node-map per sub-mesh.
T elemset(std::initializer_list< T > set) const
Combine set of element numbers for an original to the final mesh.
size_t nmesh() const
Number of sub meshes.
std::vector< size_t > m_el_offset
First element of every sub-mesh.
T nodeset(const T &set, size_t mesh_index) const
Convert set of node-numbers for a sub-mesh to the stitched mesh.
double m_rtol
Relative tolerance to find overlapping nodes.
Stitch(double rtol=1e-5, double atol=1e-8)
void push_back(const C &coor, const E &conn)
Add mesh to be stitched.
T elemset(const T &set, size_t mesh_index) const
Convert set of element-numbers for a sub-mesh to the stitched mesh.
array_type::tensor< size_t, 2 > m_conn
Connectivity [nelem, nne].
array_type::tensor< size_t, 1 > elemmap(size_t mesh_index) const
The element numbers in the stitched mesh that are coming from a specific sub-mesh.
size_t nnode() const
Number of nodes.
std::vector< array_type::tensor< size_t, 1 > > elemmap() const
Element-map per sub-mesh.
array_type::tensor< size_t, 1 > nodemap(size_t mesh_index) const
The node numbers in the stitched mesh that are coming from a specific sub-mesh.
array_type::tensor< double, 2 > m_coor
Nodal coordinates [nnode, ndim].
void push_back(const C &coor, const E &conn, const N &nodes_bot, const N &nodes_top)
Add a mesh to the top of the current stack.
Vstack(bool check_overlap=true, double rtol=1e-5, double atol=1e-8)
Class to switch between storage types.
array_type::tensor< double, 2 > AsNode(const T &arg) const
Convert "dofval" or "elemvec" to "nodevec" (overwrite entries that occur more than once).
array_type::tensor< double, 3 > AsElement(const T &arg) const
Convert "dofval" or "nodevec" to "elemvec" (overwrite entries that occur more than once).
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
#define GOOSEFEM_CHECK(expr)
Assertion that cannot be switched off.
array_type::tensor< double, 1 > w()
Integration point weights.
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
array_type::tensor< size_t, 1 > coordination(const E &conn)
Number of elements connected to each node.
array_type::tensor< size_t, 2 > dofs(size_t nnode, size_t ndim)
List with DOF-numbers in sequential order.
array_type::tensor< size_t, 2 > overlapping(const S &coor_a, const T &coor_b, double rtol=1e-5, double atol=1e-8)
Find overlapping nodes.
array_type::tensor< double, 2 > nodal_mass(const C &coor, const E &conn, ElementType type)
Compute the mass of each node in the mesh.
std::vector< std::vector< size_t > > node2dof(const D &dofs, bool sorted=true)
Nodes connected to each DOF.
array_type::tensor< double, 2 > edgesize(const C &coor, const E &conn, ElementType type)
Return size of each element edge.
T renumber(const T &dofs)
Renumber to lowest possible index (see GooseFEM::Mesh::Renumber).
array_type::tensor< double, 2 > centers(const C &coor, const E &conn, ElementType type)
Coordinates of the center of each element.
array_type::tensor< double, 1 > center_of_gravity(const C &coor, const E &conn, ElementType type)
Compute the center of gravity of a mesh.
ElementType
Enumerator for element-types.
@ Quad4
Quadrilateral: 4-noded element in 2-d.
@ Unknown
Unknown element-type.
@ Tri3
Triangle: 3-noded element in 2-d.
@ Hex8
Hexahedron: 8-noded element in 3-d.
std::vector< std::vector< size_t > > elem2node(const E &conn, bool sorted=true)
Elements connected to each node.
std::vector< std::vector< size_t > > nodaltyings(const D &dofs)
List nodal tyings based on DOF-numbers per node.
array_type::tensor< size_t, 1 > elemmap2nodemap(const T &elem_map, const C &coor, const E &conn, ElementType type)
Convert an element-map to a node-map.
ElementType defaultElementType(const S &coor, const T &conn)
Extract the element type based on the connectivity.
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.
bool is_unique(const T &arg)
Returns true is a list is unique (has not duplicate items).