FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
Mesh.h
Go to the documentation of this file.
1
9#ifndef GOOSEFEM_MESH_H
10#define GOOSEFEM_MESH_H
11
12#include "ElementQuad4.h"
13#include "MatrixDiagonal.h"
14#include "Vector.h"
15#include "assertions.h"
16#include "config.h"
17
18namespace GooseFEM {
19
23namespace Mesh {
24
25template <class D>
26inline std::vector<std::vector<size_t>> nodaltyings(const D& dofs);
27
31enum class ElementType {
32 Unknown,
33 Quad4,
34 Hex8,
35 Tri3
36};
37
45template <class S, class T>
46inline ElementType defaultElementType(const S& coor, const T& conn)
47{
48 GOOSEFEM_ASSERT(coor.dimension() == 2);
49 GOOSEFEM_ASSERT(conn.dimension() == 2);
50
51 if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
52 return ElementType::Tri3;
53 }
54 if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
55 return ElementType::Quad4;
56 }
57 if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
58 return ElementType::Hex8;
59 }
60
61 throw std::runtime_error("Element-type not implemented");
62}
63
64namespace detail {
65
66template <class T, class R>
67inline T renum(const T& arg, const R& mapping)
68{
69 T ret = T::from_shape(arg.shape());
70
71 auto jt = ret.begin();
72
73 for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) {
74 *jt = mapping(*it);
75 }
76
77 return ret;
78}
79
80} // namespace detail
81
93inline array_type::tensor<size_t, 2> dofs(size_t nnode, size_t ndim)
94{
95 return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
96}
97
114class Renumber {
115public:
116 Renumber() = default;
117
121 template <class T>
122 Renumber(const T& dofs)
123 {
124 size_t n = xt::amax(dofs)() + 1;
125 size_t i = 0;
126
127 array_type::tensor<size_t, 1> unique = xt::unique(dofs);
128
129 m_renum = xt::empty<size_t>({n});
130
131 for (auto& j : unique) {
132 m_renum(j) = i;
133 ++i;
134 }
135 }
136
143 template <class T>
144 T apply(const T& list) const
145 {
146 return detail::renum(list, m_renum);
147 }
148
157 {
158 return m_renum;
159 }
160
161private:
163};
164
171template <class T>
172inline T renumber(const T& dofs)
173{
174 return Renumber(dofs).apply(dofs);
175}
176
180template <class D>
182public:
186 using derived_type = D;
187
192 auto nelem() const
193 {
194 return derived_cast().m_nelem;
195 }
196
201 auto nnode() const
202 {
203 return derived_cast().m_nnode;
204 }
205
210 auto nne() const
211 {
212 return derived_cast().m_nne;
213 }
214
219 auto ndim() const
220 {
221 return derived_cast().m_ndim;
222 }
223
228 auto nelx() const
229 {
230 return derived_cast().nelx_impl();
231 }
232
237 auto nely() const
238 {
239 return derived_cast().nely_impl();
240 }
241
246 auto h() const
247 {
248 return derived_cast().m_h;
249 }
250
255 auto getElementType() const
256 {
257 return derived_cast().getElementType_impl();
258 }
259
264 auto coor() const
265 {
266 return derived_cast().coor_impl();
267 }
268
273 auto conn() const
274 {
275 return derived_cast().conn_impl();
276 }
277
282 auto dofs() const
283 {
284 return GooseFEM::Mesh::dofs(this->nnode(), this->ndim());
285 }
286
291 auto dofsPeriodic() const
292 {
295 array_type::tensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
296 array_type::tensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
297
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);
300 }
301
302 return GooseFEM::Mesh::renumber(ret);
303 }
304
309 auto nodesPeriodic() const
310 {
311 return derived_cast().nodesPeriodic_impl();
312 }
313
318 auto nodesOrigin() const
319 {
320 return derived_cast().nodesOrigin_impl();
321 }
322
323private:
324 auto derived_cast() -> derived_type&
325 {
326 return *static_cast<derived_type*>(this);
327 }
328
329 auto derived_cast() const -> const derived_type&
330 {
331 return *static_cast<const derived_type*>(this);
332 }
333};
334
338template <class D>
339class RegularBase2d : public RegularBase<D> {
340public:
344 using derived_type = D;
345
350 auto nodesBottomEdge() const
351 {
352 return derived_cast().nodesBottomEdge_impl();
353 }
354
359 auto nodesTopEdge() const
360 {
361 return derived_cast().nodesTopEdge_impl();
362 }
363
368 auto nodesLeftEdge() const
369 {
370 return derived_cast().nodesLeftEdge_impl();
371 }
372
377 auto nodesRightEdge() const
378 {
379 return derived_cast().nodesRightEdge_impl();
380 }
381
388 {
389 return derived_cast().nodesBottomOpenEdge_impl();
390 }
391
397 auto nodesTopOpenEdge() const
398 {
399 return derived_cast().nodesTopOpenEdge_impl();
400 }
401
407 auto nodesLeftOpenEdge() const
408 {
409 return derived_cast().nodesLeftOpenEdge_impl();
410 }
411
418 {
419 return derived_cast().nodesRightOpenEdge_impl();
420 }
421
428 {
429 return derived_cast().nodesBottomLeftCorner_impl();
430 }
431
438 {
439 return derived_cast().nodesBottomRightCorner_impl();
440 }
441
448 {
449 return derived_cast().nodesTopLeftCorner_impl();
450 }
451
458 {
459 return derived_cast().nodesTopRightCorner_impl();
460 }
461
467 {
468 return derived_cast().nodesBottomLeftCorner_impl();
469 }
470
476 {
477 return derived_cast().nodesTopLeftCorner_impl();
478 }
479
485 {
486 return derived_cast().nodesBottomRightCorner_impl();
487 }
488
494 {
495 return derived_cast().nodesTopRightCorner_impl();
496 }
497
498private:
499 auto derived_cast() -> derived_type&
500 {
501 return *static_cast<derived_type*>(this);
502 }
503
504 auto derived_cast() const -> const derived_type&
505 {
506 return *static_cast<const derived_type*>(this);
507 }
508
509 friend class RegularBase<D>;
510
511 array_type::tensor<size_t, 2> nodesPeriodic_impl() const
512 {
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);
519
520 ret(0, 0) = derived_cast().nodesBottomLeftCorner_impl();
521 ret(0, 1) = derived_cast().nodesBottomRightCorner_impl();
522
523 ret(1, 0) = derived_cast().nodesBottomLeftCorner_impl();
524 ret(1, 1) = derived_cast().nodesTopRightCorner_impl();
525
526 ret(2, 0) = derived_cast().nodesBottomLeftCorner_impl();
527 ret(2, 1) = derived_cast().nodesTopLeftCorner_impl();
528
529 size_t i = 3;
530
531 xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
532 xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
533
534 i += bot.size();
535
536 xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
537 xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
538
539 return ret;
540 }
541
542 auto nodesOrigin_impl() const
543 {
544 return derived_cast().nodesBottomLeftCorner_impl();
545 }
546};
547
551template <class D>
552class RegularBase3d : public RegularBase<D> {
553public:
557 using derived_type = D;
558
563 auto nelz() const
564 {
565 return derived_cast().nelz_impl();
566 }
567
572 auto nodesBottom() const
573 {
574 return derived_cast().nodesBottom_impl();
575 }
576
581 auto nodesTop() const
582 {
583 return derived_cast().nodesTop_impl();
584 }
585
590 auto nodesLeft() const
591 {
592 return derived_cast().nodesLeft_impl();
593 }
594
599 auto nodesRight() const
600 {
601 return derived_cast().nodesRight_impl();
602 }
603
608 auto nodesFront() const
609 {
610 return derived_cast().nodesFront_impl();
611 }
612
617 auto nodesBack() const
618 {
619 return derived_cast().nodesBack_impl();
620 }
621
628 {
629 return derived_cast().nodesFrontBottomEdge_impl();
630 }
631
637 auto nodesFrontTopEdge() const
638 {
639 return derived_cast().nodesFrontTopEdge_impl();
640 }
641
648 {
649 return derived_cast().nodesFrontLeftEdge_impl();
650 }
651
658 {
659 return derived_cast().nodesFrontRightEdge_impl();
660 }
661
668 {
669 return derived_cast().nodesBackBottomEdge_impl();
670 }
671
677 auto nodesBackTopEdge() const
678 {
679 return derived_cast().nodesBackTopEdge_impl();
680 }
681
687 auto nodesBackLeftEdge() const
688 {
689 return derived_cast().nodesBackLeftEdge_impl();
690 }
691
698 {
699 return derived_cast().nodesBackRightEdge_impl();
700 }
701
708 {
709 return derived_cast().nodesBottomLeftEdge_impl();
710 }
711
718 {
719 return derived_cast().nodesBottomRightEdge_impl();
720 }
721
727 auto nodesTopLeftEdge() const
728 {
729 return derived_cast().nodesTopLeftEdge_impl();
730 }
731
737 auto nodesTopRightEdge() const
738 {
739 return derived_cast().nodesTopRightEdge_impl();
740 }
741
747 {
748 return derived_cast().nodesFrontBottomEdge_impl();
749 }
750
756 {
757 return derived_cast().nodesBackBottomEdge_impl();
758 }
759
764 auto nodesTopFrontEdge() const
765 {
766 return derived_cast().nodesFrontTopEdge_impl();
767 }
768
773 auto nodesTopBackEdge() const
774 {
775 return derived_cast().nodesBackTopEdge_impl();
776 }
777
783 {
784 return derived_cast().nodesBottomLeftEdge_impl();
785 }
786
792 {
793 return derived_cast().nodesFrontLeftEdge_impl();
794 }
795
800 auto nodesLeftBackEdge() const
801 {
802 return derived_cast().nodesBackLeftEdge_impl();
803 }
804
809 auto nodesLeftTopEdge() const
810 {
811 return derived_cast().nodesTopLeftEdge_impl();
812 }
813
819 {
820 return derived_cast().nodesBottomRightEdge_impl();
821 }
822
827 auto nodesRightTopEdge() const
828 {
829 return derived_cast().nodesTopRightEdge_impl();
830 }
831
837 {
838 return derived_cast().nodesFrontRightEdge_impl();
839 }
840
846 {
847 return derived_cast().nodesBackRightEdge_impl();
848 }
849
856 auto nodesFrontFace() const
857 {
858 return derived_cast().nodesFrontFace_impl();
859 }
860
867 auto nodesBackFace() const
868 {
869 return derived_cast().nodesBackFace_impl();
870 }
871
878 auto nodesLeftFace() const
879 {
880 return derived_cast().nodesLeftFace_impl();
881 }
882
889 auto nodesRightFace() const
890 {
891 return derived_cast().nodesRightFace_impl();
892 }
893
900 auto nodesBottomFace() const
901 {
902 return derived_cast().nodesBottomFace_impl();
903 }
904
911 auto nodesTopFace() const
912 {
913 return derived_cast().nodesTopFace_impl();
914 }
915
921 {
922 return derived_cast().nodesFrontBottomOpenEdge_impl();
923 }
924
930 {
931 return derived_cast().nodesFrontTopOpenEdge_impl();
932 }
933
939 {
940 return derived_cast().nodesFrontLeftOpenEdge_impl();
941 }
942
948 {
949 return derived_cast().nodesFrontRightOpenEdge_impl();
950 }
951
957 {
958 return derived_cast().nodesBackBottomOpenEdge_impl();
959 }
960
966 {
967 return derived_cast().nodesBackTopOpenEdge_impl();
968 }
969
975 {
976 return derived_cast().nodesBackLeftOpenEdge_impl();
977 }
978
984 {
985 return derived_cast().nodesBackRightOpenEdge_impl();
986 }
987
993 {
994 return derived_cast().nodesBottomLeftOpenEdge_impl();
995 }
996
1002 {
1003 return derived_cast().nodesBottomRightOpenEdge_impl();
1004 }
1005
1011 {
1012 return derived_cast().nodesTopLeftOpenEdge_impl();
1013 }
1014
1020 {
1021 return derived_cast().nodesTopRightOpenEdge_impl();
1022 }
1023
1029 {
1030 return derived_cast().nodesFrontBottomOpenEdge_impl();
1031 }
1032
1038 {
1039 return derived_cast().nodesBackBottomOpenEdge_impl();
1040 }
1041
1047 {
1048 return derived_cast().nodesFrontTopOpenEdge_impl();
1049 }
1050
1056 {
1057 return derived_cast().nodesBackTopOpenEdge_impl();
1058 }
1059
1065 {
1066 return derived_cast().nodesBottomLeftOpenEdge_impl();
1067 }
1068
1074 {
1075 return derived_cast().nodesFrontLeftOpenEdge_impl();
1076 }
1077
1083 {
1084 return derived_cast().nodesBackLeftOpenEdge_impl();
1085 }
1086
1092 {
1093 return derived_cast().nodesTopLeftOpenEdge_impl();
1094 }
1095
1101 {
1102 return derived_cast().nodesBottomRightOpenEdge_impl();
1103 }
1104
1110 {
1111 return derived_cast().nodesTopRightOpenEdge_impl();
1112 }
1113
1119 {
1120 return derived_cast().nodesFrontRightOpenEdge_impl();
1121 }
1122
1128 {
1129 return derived_cast().nodesBackRightOpenEdge_impl();
1130 }
1131
1137 {
1138 return derived_cast().nodesFrontBottomLeftCorner_impl();
1139 }
1140
1146 {
1147 return derived_cast().nodesFrontBottomRightCorner_impl();
1148 }
1149
1155 {
1156 return derived_cast().nodesFrontTopLeftCorner_impl();
1157 }
1158
1164 {
1165 return derived_cast().nodesFrontTopRightCorner_impl();
1166 }
1167
1173 {
1174 return derived_cast().nodesBackBottomLeftCorner_impl();
1175 }
1176
1182 {
1183 return derived_cast().nodesBackBottomRightCorner_impl();
1184 }
1185
1191 {
1192 return derived_cast().nodesBackTopLeftCorner_impl();
1193 }
1194
1200 {
1201 return derived_cast().nodesBackTopRightCorner_impl();
1202 }
1203
1209 {
1210 return derived_cast().nodesFrontBottomLeftCorner_impl();
1211 }
1212
1218 {
1219 return derived_cast().nodesFrontBottomLeftCorner_impl();
1220 }
1221
1227 {
1228 return derived_cast().nodesFrontBottomLeftCorner_impl();
1229 }
1230
1236 {
1237 return derived_cast().nodesFrontBottomLeftCorner_impl();
1238 }
1239
1245 {
1246 return derived_cast().nodesFrontBottomLeftCorner_impl();
1247 }
1248
1254 {
1255 return derived_cast().nodesFrontBottomRightCorner_impl();
1256 }
1257
1263 {
1264 return derived_cast().nodesFrontBottomRightCorner_impl();
1265 }
1266
1272 {
1273 return derived_cast().nodesFrontBottomRightCorner_impl();
1274 }
1275
1281 {
1282 return derived_cast().nodesFrontBottomRightCorner_impl();
1283 }
1284
1290 {
1291 return derived_cast().nodesFrontBottomRightCorner_impl();
1292 }
1293
1299 {
1300 return derived_cast().nodesFrontTopLeftCorner_impl();
1301 }
1302
1308 {
1309 return derived_cast().nodesFrontTopLeftCorner_impl();
1310 }
1311
1317 {
1318 return derived_cast().nodesFrontTopLeftCorner_impl();
1319 }
1320
1326 {
1327 return derived_cast().nodesFrontTopLeftCorner_impl();
1328 }
1329
1335 {
1336 return derived_cast().nodesFrontTopLeftCorner_impl();
1337 }
1338
1344 {
1345 return derived_cast().nodesFrontTopRightCorner_impl();
1346 }
1347
1353 {
1354 return derived_cast().nodesFrontTopRightCorner_impl();
1355 }
1356
1362 {
1363 return derived_cast().nodesFrontTopRightCorner_impl();
1364 }
1365
1371 {
1372 return derived_cast().nodesFrontTopRightCorner_impl();
1373 }
1374
1380 {
1381 return derived_cast().nodesFrontTopRightCorner_impl();
1382 }
1383
1389 {
1390 return derived_cast().nodesBackBottomLeftCorner_impl();
1391 }
1392
1398 {
1399 return derived_cast().nodesBackBottomLeftCorner_impl();
1400 }
1401
1407 {
1408 return derived_cast().nodesBackBottomLeftCorner_impl();
1409 }
1410
1416 {
1417 return derived_cast().nodesBackBottomLeftCorner_impl();
1418 }
1419
1425 {
1426 return derived_cast().nodesBackBottomLeftCorner_impl();
1427 }
1428
1434 {
1435 return derived_cast().nodesBackBottomRightCorner_impl();
1436 }
1437
1443 {
1444 return derived_cast().nodesBackBottomRightCorner_impl();
1445 }
1446
1452 {
1453 return derived_cast().nodesBackBottomRightCorner_impl();
1454 }
1455
1461 {
1462 return derived_cast().nodesBackBottomRightCorner_impl();
1463 }
1464
1470 {
1471 return derived_cast().nodesBackBottomRightCorner_impl();
1472 }
1473
1479 {
1480 return derived_cast().nodesBackTopLeftCorner_impl();
1481 }
1482
1488 {
1489 return derived_cast().nodesBackTopLeftCorner_impl();
1490 }
1491
1497 {
1498 return derived_cast().nodesBackTopLeftCorner_impl();
1499 }
1500
1506 {
1507 return derived_cast().nodesBackTopLeftCorner_impl();
1508 }
1509
1515 {
1516 return derived_cast().nodesBackTopLeftCorner_impl();
1517 }
1518
1524 {
1525 return derived_cast().nodesBackTopRightCorner_impl();
1526 }
1527
1533 {
1534 return derived_cast().nodesBackTopRightCorner_impl();
1535 }
1536
1542 {
1543 return derived_cast().nodesBackTopRightCorner_impl();
1544 }
1545
1551 {
1552 return derived_cast().nodesBackTopRightCorner_impl();
1553 }
1554
1560 {
1561 return derived_cast().nodesBackTopRightCorner_impl();
1562 }
1563
1564private:
1565 auto derived_cast() -> derived_type&
1566 {
1567 return *static_cast<derived_type*>(this);
1568 }
1569
1570 auto derived_cast() const -> const derived_type&
1571 {
1572 return *static_cast<const derived_type*>(this);
1573 }
1574
1575 friend class RegularBase<D>;
1576
1577 array_type::tensor<size_t, 2> nodesPeriodic_impl() const
1578 {
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();
1585
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();
1598
1599 size_t tface = fro.size() + lft.size() + bot.size();
1600 size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size();
1601 size_t tnode = 7;
1602 array_type::tensor<size_t, 2> ret =
1603 xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
1604
1605 size_t i = 0;
1606
1607 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1608 ret(i, 1) = derived_cast().nodesFrontBottomRightCorner_impl();
1609 ++i;
1610 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1611 ret(i, 1) = derived_cast().nodesBackBottomRightCorner_impl();
1612 ++i;
1613 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1614 ret(i, 1) = derived_cast().nodesBackBottomLeftCorner_impl();
1615 ++i;
1616 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1617 ret(i, 1) = derived_cast().nodesFrontTopLeftCorner_impl();
1618 ++i;
1619 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1620 ret(i, 1) = derived_cast().nodesFrontTopRightCorner_impl();
1621 ++i;
1622 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1623 ret(i, 1) = derived_cast().nodesBackTopRightCorner_impl();
1624 ++i;
1625 ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
1626 ret(i, 1) = derived_cast().nodesBackTopLeftCorner_impl();
1627 ++i;
1628
1629 for (size_t j = 0; j < froBot.size(); ++j) {
1630 ret(i, 0) = froBot(j);
1631 ret(i, 1) = bckBot(j);
1632 ++i;
1633 }
1634 for (size_t j = 0; j < froBot.size(); ++j) {
1635 ret(i, 0) = froBot(j);
1636 ret(i, 1) = bckTop(j);
1637 ++i;
1638 }
1639 for (size_t j = 0; j < froBot.size(); ++j) {
1640 ret(i, 0) = froBot(j);
1641 ret(i, 1) = froTop(j);
1642 ++i;
1643 }
1644 for (size_t j = 0; j < botLft.size(); ++j) {
1645 ret(i, 0) = botLft(j);
1646 ret(i, 1) = botRgt(j);
1647 ++i;
1648 }
1649 for (size_t j = 0; j < botLft.size(); ++j) {
1650 ret(i, 0) = botLft(j);
1651 ret(i, 1) = topRgt(j);
1652 ++i;
1653 }
1654 for (size_t j = 0; j < botLft.size(); ++j) {
1655 ret(i, 0) = botLft(j);
1656 ret(i, 1) = topLft(j);
1657 ++i;
1658 }
1659 for (size_t j = 0; j < froLft.size(); ++j) {
1660 ret(i, 0) = froLft(j);
1661 ret(i, 1) = froRgt(j);
1662 ++i;
1663 }
1664 for (size_t j = 0; j < froLft.size(); ++j) {
1665 ret(i, 0) = froLft(j);
1666 ret(i, 1) = bckRgt(j);
1667 ++i;
1668 }
1669 for (size_t j = 0; j < froLft.size(); ++j) {
1670 ret(i, 0) = froLft(j);
1671 ret(i, 1) = bckLft(j);
1672 ++i;
1673 }
1674
1675 for (size_t j = 0; j < fro.size(); ++j) {
1676 ret(i, 0) = fro(j);
1677 ret(i, 1) = bck(j);
1678 ++i;
1679 }
1680 for (size_t j = 0; j < lft.size(); ++j) {
1681 ret(i, 0) = lft(j);
1682 ret(i, 1) = rgt(j);
1683 ++i;
1684 }
1685 for (size_t j = 0; j < bot.size(); ++j) {
1686 ret(i, 0) = bot(j);
1687 ret(i, 1) = top(j);
1688 ++i;
1689 }
1690
1691 return ret;
1692 }
1693
1694 auto nodesOrigin_impl() const
1695 {
1696 return derived_cast().nodesFrontBottomLeftCorner_impl();
1697 }
1698};
1699
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)
1717{
1718 GOOSEFEM_ASSERT(coor_a.dimension() == 2);
1719 GOOSEFEM_ASSERT(coor_b.dimension() == 2);
1720 GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
1721
1722 std::vector<size_t> ret_a;
1723 std::vector<size_t> ret_b;
1724
1725 for (size_t i = 0; i < coor_a.shape(0); ++i) {
1726
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)));
1729
1730 for (auto& j : idx) {
1731 ret_a.push_back(i);
1732 ret_b.push_back(j);
1733 }
1734 }
1735
1736 array_type::tensor<size_t, 2> ret = xt::empty<size_t>({size_t(2), ret_a.size()});
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];
1740 }
1741
1742 return ret;
1743}
1744
1749public:
1750 ManualStitch() = default;
1751
1763 template <class CA, class EA, class NA, class CB, class EB, class NB>
1765 const CA& coor_a,
1766 const EA& conn_a,
1767 const NA& overlapping_nodes_a,
1768 const CB& coor_b,
1769 const EB& conn_b,
1770 const NB& overlapping_nodes_b,
1771 bool check_position = true,
1772 double rtol = 1e-5,
1773 double atol = 1e-8)
1774 {
1775 UNUSED(rtol);
1776 UNUSED(atol);
1777
1778 GOOSEFEM_ASSERT(coor_a.dimension() == 2);
1779 GOOSEFEM_ASSERT(conn_a.dimension() == 2);
1780 GOOSEFEM_ASSERT(overlapping_nodes_a.dimension() == 1);
1781 GOOSEFEM_ASSERT(coor_b.dimension() == 2);
1782 GOOSEFEM_ASSERT(conn_b.dimension() == 2);
1783 GOOSEFEM_ASSERT(overlapping_nodes_b.dimension() == 1);
1784 GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape()));
1785 GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
1786 GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1));
1787
1788 if (check_position) {
1789 GOOSEFEM_CHECK(xt::allclose(
1790 xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()),
1791 xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()),
1792 rtol,
1793 atol));
1794 }
1795
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();
1800
1801 size_t nela = conn_a.shape(0);
1802 size_t nelb = conn_b.shape(0);
1803 size_t nne = conn_a.shape(1);
1804
1805 m_nel_a = nela;
1806 m_nel_b = nelb;
1807 m_nnd_a = nnda;
1808
1810 xt::setdiff1d(xt::arange<size_t>(nndb), overlapping_nodes_b);
1811
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;
1815
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);
1819
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());
1824 }
1825
1830 size_t nmesh() const
1831 {
1832 return 2;
1833 }
1834
1839 size_t nelem() const
1840 {
1841 return m_conn.shape(0);
1842 }
1843
1848 size_t nnode() const
1849 {
1850 return m_coor.shape(0);
1851 }
1852
1857 size_t nne() const
1858 {
1859 return m_conn.shape(1);
1860 }
1861
1866 size_t ndim() const
1867 {
1868 return m_coor.shape(1);
1869 }
1870
1876 {
1877 return m_coor;
1878 }
1879
1885 {
1886 return m_conn;
1887 }
1888
1894 {
1895 size_t nnode = this->nnode();
1896 size_t ndim = this->ndim();
1897 return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
1898 }
1899
1904 std::vector<array_type::tensor<size_t, 1>> nodemap() const
1905 {
1906 std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
1907 for (size_t i = 0; i < this->nmesh(); ++i) {
1908 ret[i] = this->nodemap(i);
1909 }
1910 return ret;
1911 }
1912
1917 std::vector<array_type::tensor<size_t, 1>> elemmap() const
1918 {
1919 std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
1920 for (size_t i = 0; i < this->nmesh(); ++i) {
1921 ret[i] = this->elemmap(i);
1922 }
1923 return ret;
1924 }
1925
1931 {
1932 GOOSEFEM_ASSERT(mesh_index <= 1);
1933
1934 if (mesh_index == 0) {
1935 return xt::arange<size_t>(m_nnd_a);
1936 }
1937
1938 return m_map_b;
1939 }
1940
1946 {
1947 GOOSEFEM_ASSERT(mesh_index <= 1);
1948
1949 if (mesh_index == 0) {
1950 return xt::arange<size_t>(m_nel_a);
1951 }
1952
1953 return xt::arange<size_t>(m_nel_b) + m_nel_a;
1954 }
1955
1963 template <class T>
1964 T nodeset(const T& set, size_t mesh_index) const
1965 {
1966 GOOSEFEM_ASSERT(mesh_index <= 1);
1967
1968 if (mesh_index == 0) {
1969 GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a);
1970 return set;
1971 }
1972
1973 GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size());
1974 return detail::renum(set, m_map_b);
1975 }
1976
1984 template <class T>
1985 T elemset(const T& set, size_t mesh_index) const
1986 {
1987 GOOSEFEM_ASSERT(mesh_index <= 1);
1988
1989 if (mesh_index == 0) {
1990 GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a);
1991 return set;
1992 }
1993
1994 GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b);
1995 return set + m_nel_a;
1996 }
1997
1998private:
2002 size_t m_nnd_a;
2003 size_t m_nel_a;
2004 size_t m_nel_b;
2005};
2006
2010class Stitch {
2011public:
2016 Stitch(double rtol = 1e-5, double atol = 1e-8)
2017 {
2018 m_rtol = rtol;
2019 m_atol = atol;
2020 }
2021
2028 template <class C, class E>
2029 void push_back(const C& coor, const E& conn)
2030 {
2031 GOOSEFEM_ASSERT(coor.dimension() == 2);
2032 GOOSEFEM_ASSERT(conn.dimension() == 2);
2033
2034 if (m_map.size() == 0) {
2035 m_coor = coor;
2036 m_conn = conn;
2037 m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
2038 m_nel.push_back(conn.shape(0));
2039 m_el_offset.push_back(0);
2040 return;
2041 }
2042
2043 auto overlap = overlapping(m_coor, coor, m_rtol, m_atol);
2044 size_t index = m_map.size();
2045
2046 ManualStitch stitch(
2047 m_coor,
2048 m_conn,
2049 xt::eval(xt::view(overlap, 0, xt::all())),
2050 coor,
2051 conn,
2052 xt::eval(xt::view(overlap, 1, xt::all())),
2053 false);
2054
2055 m_coor = stitch.coor();
2056 m_conn = stitch.conn();
2057 m_map.push_back(stitch.nodemap(1));
2058 m_nel.push_back(conn.shape(0));
2059 m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
2060 }
2061
2066 size_t nmesh() const
2067 {
2068 return m_map.size();
2069 }
2070
2075 size_t nelem() const
2076 {
2077 return m_conn.shape(0);
2078 }
2079
2084 size_t nnode() const
2085 {
2086 return m_coor.shape(0);
2087 }
2088
2093 size_t nne() const
2094 {
2095 return m_conn.shape(1);
2096 }
2097
2102 size_t ndim() const
2103 {
2104 return m_coor.shape(1);
2105 }
2106
2112 {
2113 return m_coor;
2114 }
2115
2121 {
2122 return m_conn;
2123 }
2124
2130 {
2131 size_t nnode = this->nnode();
2132 size_t ndim = this->ndim();
2133 return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
2134 }
2135
2140 std::vector<array_type::tensor<size_t, 1>> nodemap() const
2141 {
2142 std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
2143 for (size_t i = 0; i < this->nmesh(); ++i) {
2144 ret[i] = this->nodemap(i);
2145 }
2146 return ret;
2147 }
2148
2153 std::vector<array_type::tensor<size_t, 1>> elemmap() const
2154 {
2155 std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
2156 for (size_t i = 0; i < this->nmesh(); ++i) {
2157 ret[i] = this->elemmap(i);
2158 }
2159 return ret;
2160 }
2161
2169 {
2170 GOOSEFEM_ASSERT(mesh_index < m_map.size());
2171 return m_map[mesh_index];
2172 }
2173
2181 {
2182 GOOSEFEM_ASSERT(mesh_index < m_map.size());
2183 return xt::arange<size_t>(m_nel[mesh_index]) + m_el_offset[mesh_index];
2184 }
2185
2193 template <class T>
2194 T nodeset(const T& set, size_t mesh_index) const
2195 {
2196 GOOSEFEM_ASSERT(mesh_index < m_map.size());
2197 GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size());
2198 return detail::renum(set, m_map[mesh_index]);
2199 }
2200
2208 template <class T>
2209 T elemset(const T& set, size_t mesh_index) const
2210 {
2211 GOOSEFEM_ASSERT(mesh_index < m_map.size());
2212 GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]);
2213 return set + m_el_offset[mesh_index];
2214 }
2215
2222 template <class T>
2223 T nodeset(const std::vector<T>& set) const
2224 {
2225 GOOSEFEM_ASSERT(set.size() == m_map.size());
2226
2227 size_t n = 0;
2228
2229 for (size_t i = 0; i < set.size(); ++i) {
2230 n += set[i].size();
2231 }
2232
2233 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({n});
2234
2235 n = 0;
2236
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);
2239 n += set[i].size();
2240 }
2241
2242 return xt::unique(ret);
2243 }
2244
2248 template <class T>
2249 T nodeset(std::initializer_list<T> set) const
2250 {
2251 return this->nodeset(std::vector<T>(set));
2252 }
2253
2260 template <class T>
2261 T elemset(const std::vector<T>& set) const
2262 {
2263 GOOSEFEM_ASSERT(set.size() == m_map.size());
2264
2265 size_t n = 0;
2266
2267 for (size_t i = 0; i < set.size(); ++i) {
2268 n += set[i].size();
2269 }
2270
2271 array_type::tensor<size_t, 1> ret = xt::empty<size_t>({n});
2272
2273 n = 0;
2274
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);
2277 n += set[i].size();
2278 }
2279
2280 return ret;
2281 }
2282
2286 template <class T>
2287 T elemset(std::initializer_list<T> set) const
2288 {
2289 return this->elemset(std::vector<T>(set));
2290 }
2291
2292protected:
2295 std::vector<array_type::tensor<size_t, 1>> m_map;
2296 std::vector<size_t> m_nel;
2297 std::vector<size_t> m_el_offset;
2298 double m_rtol;
2299 double m_atol;
2300};
2301
2305class Vstack : public Stitch {
2306public:
2312 Vstack(bool check_overlap = true, double rtol = 1e-5, double atol = 1e-8)
2313 {
2314 m_check_overlap = check_overlap;
2315 m_rtol = rtol;
2316 m_atol = atol;
2317 }
2318
2328 template <class C, class E, class N>
2329 void push_back(const C& coor, const E& conn, const N& nodes_bot, const N& nodes_top)
2330 {
2331 if (m_map.size() == 0) {
2332 m_coor = coor;
2333 m_conn = conn;
2334 m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
2335 m_nel.push_back(conn.shape(0));
2336 m_el_offset.push_back(0);
2337 m_nodes_bot.push_back(nodes_bot);
2338 m_nodes_top.push_back(nodes_top);
2339 return;
2340 }
2341
2342 GOOSEFEM_ASSERT(nodes_bot.size() == m_nodes_top.back().size());
2343
2344 size_t index = m_map.size();
2345
2346 double shift = xt::amax(xt::view(m_coor, xt::all(), 1))();
2347 auto x = coor;
2348 xt::view(x, xt::all(), 1) += shift;
2349
2350 ManualStitch stitch(
2351 m_coor,
2352 m_conn,
2353 m_nodes_top.back(),
2354 x,
2355 conn,
2356 nodes_bot,
2357 m_check_overlap,
2358 m_rtol,
2359 m_atol);
2360
2361 m_nodes_bot.push_back(stitch.nodeset(nodes_bot, 1));
2362 m_nodes_top.push_back(stitch.nodeset(nodes_top, 1));
2363
2364 m_coor = stitch.coor();
2365 m_conn = stitch.conn();
2366 m_map.push_back(stitch.nodemap(1));
2367 m_nel.push_back(conn.shape(0));
2368 m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
2369 }
2370
2371private:
2372 std::vector<array_type::tensor<size_t, 1>>
2373 m_nodes_bot;
2374 std::vector<array_type::tensor<size_t, 1>>
2375 m_nodes_top;
2376 bool m_check_overlap;
2377};
2378
2387class Reorder {
2388public:
2389 Reorder() = default;
2390
2394 template <class T>
2395 Reorder(const std::initializer_list<T> args)
2396 {
2397 size_t n = 0;
2398 size_t i = 0;
2399
2400 for (auto& arg : args) {
2401 if (arg.size() == 0) {
2402 continue;
2403 }
2404 n = std::max(n, xt::amax(arg)() + 1);
2405 }
2406
2407#ifdef GOOSEFEM_ENABLE_ASSERT
2408 for (auto& arg : args) {
2410 }
2411#endif
2412
2413 m_renum = xt::empty<size_t>({n});
2414
2415 for (auto& arg : args) {
2416 for (auto& j : arg) {
2417 m_renum(j) = i;
2418 ++i;
2419 }
2420 }
2421 }
2422
2429 template <class T>
2430 T apply(const T& list) const
2431 {
2432 T ret = T::from_shape(list.shape());
2433
2434 auto jt = ret.begin();
2435
2436 for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
2437 *jt = m_renum(*it);
2438 }
2439
2440 return ret;
2441 }
2442
2451 {
2452 return m_renum;
2453 }
2454
2455private:
2457};
2458
2465template <class E>
2467{
2468 GOOSEFEM_ASSERT(conn.dimension() == 2);
2469
2470 size_t nnode = xt::amax(conn)() + 1;
2471
2472 array_type::tensor<size_t, 1> N = xt::zeros<size_t>({nnode});
2473
2474 for (auto it = conn.begin(); it != conn.end(); ++it) {
2475 N(*it) += 1;
2476 }
2477
2478 return N;
2479}
2480
2488template <class E>
2489inline std::vector<std::vector<size_t>> elem2node(const E& conn, bool sorted = true)
2490{
2491 auto N = coordination(conn);
2492 auto nnode = N.size();
2493
2494 std::vector<std::vector<size_t>> ret(nnode);
2495 for (size_t i = 0; i < nnode; ++i) {
2496 ret[i].reserve(N(i));
2497 }
2498
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);
2502 }
2503 }
2504
2505 if (sorted) {
2506 for (auto& row : ret) {
2507 std::sort(row.begin(), row.end());
2508 }
2509 }
2510
2511 return ret;
2512}
2513
2519template <class E, class D>
2520inline std::vector<std::vector<size_t>> elem2node(const E& conn, const D& dofs, bool sorted = true)
2521{
2522 size_t nnode = dofs.shape(0);
2523 auto ret = elem2node(conn, sorted);
2524 auto nties = nodaltyings(dofs);
2525
2526 for (size_t m = 0; m < nnode; ++m) {
2527 if (nties[m].size() <= 1) {
2528 continue;
2529 }
2530 if (m > nties[m][0]) {
2531 continue;
2532 }
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();
2536 }
2537 ret[m].reserve(n);
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());
2540 }
2541 if (sorted) {
2542 std::sort(ret[m].begin(), ret[m].end());
2543 }
2544 for (size_t j = 1; j < nties[m].size(); ++j) {
2545 ret[nties[m][j]] = ret[m];
2546 }
2547 }
2548
2549 return ret;
2550}
2551
2559template <class D>
2560inline std::vector<std::vector<size_t>> node2dof(const D& dofs, bool sorted = true)
2561{
2562 return elem2node(dofs, sorted);
2563}
2564
2573template <class C, class E>
2574inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn, ElementType type)
2575{
2576 GOOSEFEM_ASSERT(coor.dimension() == 2);
2577 GOOSEFEM_ASSERT(conn.dimension() == 2);
2578 GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
2579
2580 if (type == ElementType::Quad4) {
2581 GOOSEFEM_ASSERT(coor.shape(1) == 2ul);
2582 GOOSEFEM_ASSERT(conn.shape(1) == 4ul);
2583 array_type::tensor<size_t, 1> n0 = xt::view(conn, xt::all(), 0);
2584 array_type::tensor<size_t, 1> n1 = xt::view(conn, xt::all(), 1);
2585 array_type::tensor<size_t, 1> n2 = xt::view(conn, xt::all(), 2);
2586 array_type::tensor<size_t, 1> n3 = xt::view(conn, xt::all(), 3);
2587 array_type::tensor<double, 1> x0 = xt::view(coor, xt::keep(n0), 0);
2588 array_type::tensor<double, 1> x1 = xt::view(coor, xt::keep(n1), 0);
2589 array_type::tensor<double, 1> x2 = xt::view(coor, xt::keep(n2), 0);
2590 array_type::tensor<double, 1> x3 = xt::view(coor, xt::keep(n3), 0);
2591 array_type::tensor<double, 1> y0 = xt::view(coor, xt::keep(n0), 1);
2592 array_type::tensor<double, 1> y1 = xt::view(coor, xt::keep(n1), 1);
2593 array_type::tensor<double, 1> y2 = xt::view(coor, xt::keep(n2), 1);
2594 array_type::tensor<double, 1> y3 = xt::view(coor, xt::keep(n3), 1);
2595 array_type::tensor<double, 2> ret = xt::empty<double>(conn.shape());
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));
2600 return ret;
2601 }
2602
2603 throw std::runtime_error("Element-type not implemented");
2604}
2605
2614template <class C, class E>
2615inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn)
2616{
2617 return edgesize(coor, conn, defaultElementType(coor, conn));
2618}
2619
2628template <class C, class E>
2629inline array_type::tensor<double, 2> centers(const C& coor, const E& conn, ElementType type)
2630{
2631 GOOSEFEM_ASSERT(coor.dimension() == 2);
2632 GOOSEFEM_ASSERT(conn.dimension() == 2);
2633 GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
2634 array_type::tensor<double, 2> ret = xt::zeros<double>({conn.shape(0), coor.shape(1)});
2635
2636 if (type == ElementType::Quad4) {
2637 GOOSEFEM_ASSERT(coor.shape(1) == 2);
2638 GOOSEFEM_ASSERT(conn.shape(1) == 4);
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());
2642 }
2643 ret /= 4.0;
2644 return ret;
2645 }
2646
2647 throw std::runtime_error("Element-type not implemented");
2648}
2649
2658template <class C, class E>
2659inline array_type::tensor<double, 2> centers(const C& coor, const E& conn)
2660{
2661 return centers(coor, conn, defaultElementType(coor, conn));
2662}
2663
2673template <class T, class C, class E>
2675elemmap2nodemap(const T& elem_map, const C& coor, const E& conn, ElementType type)
2676{
2677 GOOSEFEM_ASSERT(elem_map.dimension() == 1);
2678 GOOSEFEM_ASSERT(coor.dimension() == 2);
2679 GOOSEFEM_ASSERT(conn.dimension() == 2);
2680 GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
2681 GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0));
2682 size_t N = coor.shape(0);
2683
2684 array_type::tensor<size_t, 1> ret = N * xt::ones<size_t>({N});
2685
2686 if (type == ElementType::Quad4) {
2687 GOOSEFEM_ASSERT(coor.shape(1) == 2);
2688 GOOSEFEM_ASSERT(conn.shape(1) == 4);
2689
2690 for (size_t i = 0; i < 4; ++i) {
2691 array_type::tensor<size_t, 1> t = N * xt::ones<size_t>({N});
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);
2696 }
2697
2698 return ret;
2699 }
2700
2701 throw std::runtime_error("Element-type not implemented");
2702}
2703
2713template <class T, class C, class E>
2715elemmap2nodemap(const T& elem_map, const C& coor, const E& conn)
2716{
2717 return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn));
2718}
2719
2734template <class C, class E>
2735inline array_type::tensor<double, 2> nodal_mass(const C& coor, const E& conn, ElementType type)
2736{
2737 auto dof = dofs(coor.shape(0), coor.shape(1));
2738 GooseFEM::MatrixDiagonal M(conn, dof);
2739 GooseFEM::Vector vector(conn, dof);
2740 array_type::tensor<double, 2> rho = xt::ones<double>(conn.shape());
2741
2742 if (type == ElementType::Quad4) {
2744 vector.AsElement(coor),
2747 M.assemble(quad.Int_N_scalar_NT_dV(rho));
2748 }
2749 else {
2750 throw std::runtime_error("Element-type not implemented");
2751 }
2752
2753 return vector.AsNode(M.data());
2754}
2755
2769template <class C, class E>
2770inline array_type::tensor<double, 2> nodal_mass(const C& coor, const E& conn)
2771{
2772 return nodal_mass(coor, conn, defaultElementType(coor, conn));
2773}
2774
2775namespace detail {
2776
2777// todo: remove after upstream fix
2778template <class T>
2779array_type::tensor<double, 1> average_axis_0(const T& data, const T& weights)
2780{
2781 GOOSEFEM_ASSERT(data.dimension() == 2);
2782 GOOSEFEM_ASSERT(xt::has_shape(data, weights.shape()));
2783
2784 array_type::tensor<double, 1> ret = xt::zeros<double>({weights.shape(1)});
2785
2786 for (size_t j = 0; j < weights.shape(1); ++j) {
2787 double norm = 0.0;
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);
2791 }
2792 ret(j) /= norm;
2793 }
2794 return ret;
2795}
2796
2797} // namespace detail
2798
2809template <class C, class E>
2810inline array_type::tensor<double, 1>
2811center_of_gravity(const C& coor, const E& conn, ElementType type)
2812{
2813 // todo: remove after upstream fix
2814 return detail::average_axis_0(coor, nodal_mass(coor, conn, type));
2815 // return xt::average(coor, nodal_mass(coor, conn, type), 0);
2816}
2817
2827template <class C, class E>
2828inline array_type::tensor<double, 1> center_of_gravity(const C& coor, const E& conn)
2829{
2830 // todo: remove after upstream fix
2831 return detail::average_axis_0(coor, nodal_mass(coor, conn, defaultElementType(coor, conn)));
2832 // return xt::average(coor, nodal_mass(coor, conn, defaultElementType(coor, conn)), 0);
2833}
2834
2841template <class D>
2842inline std::vector<std::vector<size_t>> nodaltyings(const D& dofs)
2843{
2844 size_t nnode = dofs.shape(0);
2845 size_t ndim = dofs.shape(1);
2846 auto nodemap = node2dof(dofs);
2847 std::vector<std::vector<size_t>> ret(nnode);
2848
2849 for (size_t m = 0; m < nnode; ++m) {
2850 auto r = nodemap[dofs(m, 0)];
2851 std::sort(r.begin(), r.end());
2852 ret[m] = r;
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());
2857 GOOSEFEM_ASSERT(std::equal(r.begin(), r.end(), t.begin()));
2858 }
2859#endif
2860 }
2861
2862 return ret;
2863}
2864
2865} // namespace Mesh
2866} // namespace GooseFEM
2867
2868#endif
Quadrature for 4-noded quadrilateral element in 2d (GooseFEM::Mesh::ElementType::Quad4),...
Diagonal matrix.
Methods to switch between storage types based on a mesh.
Interpolation and quadrature.
Definition: ElementQuad4.h:237
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.
Definition: Element.h:768
void assemble(const T &elemmat)
Assemble from "elemmat".
Definition: Matrix.h:327
const array_type::tensor< double, 1 > & data() const
Underlying matrix.
Stitch two mesh objects, specifying overlapping nodes by hand.
Definition: Mesh.h:1748
T elemset(const T &set, size_t mesh_index) const
Convert set of element numbers for an original mesh to the stitched mesh.
Definition: Mesh.h:1985
array_type::tensor< size_t, 1 > nodemap(size_t mesh_index) const
Definition: Mesh.h:1930
std::vector< array_type::tensor< size_t, 1 > > elemmap() const
Element-map per sub-mesh.
Definition: Mesh.h:1917
const array_type::tensor< size_t, 2 > & conn() const
Connectivity [nelem, nne].
Definition: Mesh.h:1884
array_type::tensor< size_t, 2 > dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
Definition: Mesh.h:1893
T nodeset(const T &set, size_t mesh_index) const
Convert set of node numbers for an original mesh to the stitched mesh.
Definition: Mesh.h:1964
size_t nne() const
Number of nodes-per-element.
Definition: Mesh.h:1857
size_t nmesh() const
Number of sub meshes == 2.
Definition: Mesh.h:1830
array_type::tensor< size_t, 1 > elemmap(size_t mesh_index) const
Definition: Mesh.h:1945
std::vector< array_type::tensor< size_t, 1 > > nodemap() const
Node-map per sub-mesh.
Definition: Mesh.h:1904
const array_type::tensor< double, 2 > & coor() const
Nodal coordinates [nnode, ndim].
Definition: Mesh.h:1875
size_t nnode() const
Number of nodes.
Definition: Mesh.h:1848
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)
Definition: Mesh.h:1764
size_t nelem() const
Number of elements.
Definition: Mesh.h:1839
size_t ndim() const
Number of dimensions.
Definition: Mesh.h:1866
CRTP base class for regular meshes in 2d.
Definition: Mesh.h:339
auto nodesTopRightCorner() const
The top-right corner node (at x = nelx * h, y = nely * h).
Definition: Mesh.h:457
auto nodesTopLeftCorner() const
The top-left corner node (at x = 0, y = nely * h).
Definition: Mesh.h:447
auto nodesLeftEdge() const
Nodes along the left edge (x = 0), in order of increasing y.
Definition: Mesh.h:368
auto nodesBottomEdge() const
Nodes along the bottom edge (y = 0), in order of increasing x.
Definition: Mesh.h:350
auto nodesRightOpenEdge() const
Nodes along the right edge (x = nelx * h), without the corners (at y = 0 and y = nely * h).
Definition: Mesh.h:417
auto nodesRightBottomCorner() const
Alias of nodesBottomRightCorner().
Definition: Mesh.h:484
auto nodesRightEdge() const
Nodes along the right edge (x = nelx * h), in order of increasing y.
Definition: Mesh.h:377
auto nodesBottomLeftCorner() const
The bottom-left corner node (at x = 0, y = 0).
Definition: Mesh.h:427
auto nodesTopEdge() const
Nodes along the top edge (y = nely * h), in order of increasing x.
Definition: Mesh.h:359
auto nodesBottomRightCorner() const
The bottom-right corner node (at x = nelx * h, y = 0).
Definition: Mesh.h:437
auto nodesRightTopCorner() const
Alias of nodesTopRightCorner().
Definition: Mesh.h:493
auto nodesLeftBottomCorner() const
Alias of nodesBottomLeftCorner().
Definition: Mesh.h:466
auto nodesLeftOpenEdge() const
Nodes along the left edge (x = 0), without the corners (at y = 0 and y = nely * h).
Definition: Mesh.h:407
auto nodesBottomOpenEdge() const
Nodes along the bottom edge (y = 0), without the corners (at x = 0 and x = nelx * h).
Definition: Mesh.h:387
auto nodesLeftTopCorner() const
Alias of nodesTopLeftCorner().
Definition: Mesh.h:475
D derived_type
Underlying type.
Definition: Mesh.h:344
auto nodesTopOpenEdge() const
Nodes along the top edge (y = nely * h), without the corners (at x = 0 and x = nelx * h).
Definition: Mesh.h:397
CRTP base class for regular meshes in 3d.
Definition: Mesh.h:552
auto nodesRightBackOpenEdge() const
Alias of nodesBackRightOpenEdge().
Definition: Mesh.h:1127
auto nodesFrontRightTopCorner() const
Alias of nodesFrontTopRightCorner().
Definition: Mesh.h:1343
auto nodesRightBottomFrontCorner() const
Alias of nodesFrontBottomRightCorner().
Definition: Mesh.h:1289
auto nodesFrontRightOpenEdge() const
Same as nodesFrontRightEdge() but without corners.
Definition: Mesh.h:947
auto nodesFrontTopLeftCorner() const
Front-Top-Left corner node.
Definition: Mesh.h:1154
auto nodesFrontBottomRightCorner() const
Front-Bottom-Right corner node.
Definition: Mesh.h:1145
auto nodesFrontLeftTopCorner() const
Alias of nodesFrontTopLeftCorner().
Definition: Mesh.h:1298
auto nodesTopBackRightCorner() const
Alias of nodesBackTopRightCorner().
Definition: Mesh.h:1532
auto nodesFrontTopEdge() const
Nodes along the edge at the intersection of the front and top faces (z = 0 and y = nely * h).
Definition: Mesh.h:637
auto nodesFrontLeftOpenEdge() const
Same as nodesFrontLeftEdge() but without corners.
Definition: Mesh.h:938
auto nodesTopLeftEdge() const
Nodes along the edge at the intersection of the top and left faces (y = 0 and x = nelx * h).
Definition: Mesh.h:727
auto nodesBackLeftEdge() const
Nodes along the edge at the intersection of the back and left faces (z = nelz * h and x = nelx * h).
Definition: Mesh.h:687
auto nodesLeftBottomEdge() const
Alias of nodesBottomLeftEdge()
Definition: Mesh.h:782
auto nodesFrontLeftBottomCorner() const
Alias of nodesFrontBottomLeftCorner().
Definition: Mesh.h:1208
auto nodesRightTopEdge() const
Alias of nodesTopRightEdge()
Definition: Mesh.h:827
auto nodesBackRightBottomCorner() const
Alias of nodesBackBottomRightCorner().
Definition: Mesh.h:1433
auto nodesBottomRightOpenEdge() const
Same as nodesBottomRightEdge() but without corners.
Definition: Mesh.h:1001
auto nodesFrontBottomEdge() const
Nodes along the edge at the intersection of the front and bottom faces (z = 0 and y = 0).
Definition: Mesh.h:627
auto nodesTopRightOpenEdge() const
Same as nodesTopRightEdge() but without corners.
Definition: Mesh.h:1019
auto nodesFrontBottomOpenEdge() const
Same as nodesFrontBottomEdge() but without corners.
Definition: Mesh.h:920
auto nodesBack() const
Nodes along the back face (z = nelz * h).
Definition: Mesh.h:617
auto nodesBottomRightFrontCorner() const
Alias of nodesFrontBottomRightCorner().
Definition: Mesh.h:1271
auto nodesTopRightBackCorner() const
Alias of nodesBackTopRightCorner().
Definition: Mesh.h:1541
auto nelz() const
Number of elements in y-direction == height of the mesh, in units of h,.
Definition: Mesh.h:563
auto nodesRightBottomOpenEdge() const
Alias of nodesBottomRightOpenEdge().
Definition: Mesh.h:1100
auto nodesLeftBottomFrontCorner() const
Alias of nodesFrontBottomLeftCorner().
Definition: Mesh.h:1244
auto nodesLeftTopEdge() const
Alias of nodesTopLeftEdge()
Definition: Mesh.h:809
auto nodesBackRightTopCorner() const
Alias of nodesBackTopRightCorner().
Definition: Mesh.h:1523
auto nodesTopLeftFrontCorner() const
Alias of nodesFrontTopLeftCorner().
Definition: Mesh.h:1316
auto nodesLeftFrontEdge() const
Alias of nodesFrontLeftEdge()
Definition: Mesh.h:791
auto nodesFrontTopOpenEdge() const
Same as nodesFrontTopEdge() but without corners.
Definition: Mesh.h:929
auto nodesTopLeftOpenEdge() const
Same as nodesTopLeftEdge() but without corners.
Definition: Mesh.h:1010
auto nodesTopLeftBackCorner() const
Alias of nodesBackTopLeftCorner().
Definition: Mesh.h:1496
auto nodesLeftFrontBottomCorner() const
Alias of nodesFrontBottomLeftCorner().
Definition: Mesh.h:1235
auto nodesTopFrontLeftCorner() const
Alias of nodesFrontTopLeftCorner().
Definition: Mesh.h:1307
auto nodesLeftFrontOpenEdge() const
Alias of nodesFrontLeftOpenEdge().
Definition: Mesh.h:1073
auto nodesBottomBackLeftCorner() const
Alias of nodesBackBottomLeftCorner().
Definition: Mesh.h:1397
auto nodesTopRightFrontCorner() const
Alias of nodesFrontTopRightCorner().
Definition: Mesh.h:1361
auto nodesBackLeftBottomCorner() const
Alias of nodesBackBottomLeftCorner().
Definition: Mesh.h:1388
auto nodesBottomLeftEdge() const
Nodes along the edge at the intersection of the bottom and left faces (y = 0 and x = 0).
Definition: Mesh.h:707
auto nodesTopBackOpenEdge() const
Alias of nodesBackTopOpenEdge().
Definition: Mesh.h:1055
auto nodesBottomBackOpenEdge() const
Alias of nodesBackBottomOpenEdge().
Definition: Mesh.h:1037
auto nodesBottomLeftOpenEdge() const
Same as nodesBottomLeftEdge() but without corners.
Definition: Mesh.h:992
auto nodesBackTopOpenEdge() const
Same as nodesBackTopEdge() but without corners.
Definition: Mesh.h:965
auto nodesLeftFace() const
Nodes along the left face excluding edges.
Definition: Mesh.h:878
auto nodesRightBottomEdge() const
Alias of nodesBottomRightEdge()
Definition: Mesh.h:818
auto nodesBackRightEdge() const
Nodes along the edge at the intersection of the back and right faces (? = nelz * h and ?...
Definition: Mesh.h:697
auto nodesRightTopBackCorner() const
Alias of nodesBackTopRightCorner().
Definition: Mesh.h:1559
auto nodesBottomLeftBackCorner() const
Alias of nodesBackBottomLeftCorner().
Definition: Mesh.h:1406
auto nodesTopFace() const
Nodes along the top face excluding edges.
Definition: Mesh.h:911
auto nodesLeftBackTopCorner() const
Alias of nodesBackTopLeftCorner().
Definition: Mesh.h:1505
auto nodesLeft() const
Nodes along the left face (x = 0).
Definition: Mesh.h:590
auto nodesTopBackLeftCorner() const
Alias of nodesBackTopLeftCorner().
Definition: Mesh.h:1487
auto nodesBackTopLeftCorner() const
Back-Top-Left corner node.
Definition: Mesh.h:1190
auto nodesTopFrontOpenEdge() const
Alias of nodesFrontTopOpenEdge().
Definition: Mesh.h:1046
auto nodesRightTopFrontCorner() const
Alias of nodesFrontTopRightCorner().
Definition: Mesh.h:1379
auto nodesTopFrontEdge() const
Alias of nodesFrontTopEdge()
Definition: Mesh.h:764
auto nodesRightFace() const
Nodes along the right face excluding edges.
Definition: Mesh.h:889
auto nodesLeftBottomBackCorner() const
Alias of nodesBackBottomLeftCorner().
Definition: Mesh.h:1424
auto nodesBottomRightEdge() const
Nodes along the edge at the intersection of the bottom and right faces (y = 0 and x = nelx * h).
Definition: Mesh.h:717
auto nodesFrontRightEdge() const
Nodes along the edge at the intersection of the front and right faces (z = 0 and x = nelx * h).
Definition: Mesh.h:657
auto nodesTopRightEdge() const
Nodes along the edge at the intersection of the top and right faces (y = nely * h and x = nelx * h).
Definition: Mesh.h:737
auto nodesBackTopEdge() const
Nodes along the edge at the intersection of the back and top faces (z = nelz * h and x = 0).
Definition: Mesh.h:677
auto nodesLeftBottomOpenEdge() const
Alias of nodesBottomLeftOpenEdge().
Definition: Mesh.h:1064
auto nodesRightFrontEdge() const
Alias of nodesFrontRightEdge()
Definition: Mesh.h:836
auto nodesBackBottomLeftCorner() const
Back-Bottom-Left corner node.
Definition: Mesh.h:1172
auto nodesFrontFace() const
Nodes along the front face excluding edges.
Definition: Mesh.h:856
auto nodesBackLeftTopCorner() const
Alias of nodesBackTopLeftCorner().
Definition: Mesh.h:1478
auto nodesFrontTopRightCorner() const
Front-Top-Right corner node.
Definition: Mesh.h:1163
auto nodesRightBackBottomCorner() const
Alias of nodesBackBottomRightCorner().
Definition: Mesh.h:1460
auto nodesLeftFrontTopCorner() const
Alias of nodesFrontTopLeftCorner().
Definition: Mesh.h:1325
auto nodesRightBackTopCorner() const
Alias of nodesBackTopRightCorner().
Definition: Mesh.h:1550
auto nodesLeftTopFrontCorner() const
Alias of nodesFrontTopLeftCorner().
Definition: Mesh.h:1334
auto nodesRightTopOpenEdge() const
Alias of nodesTopRightOpenEdge().
Definition: Mesh.h:1109
auto nodesBackRightOpenEdge() const
Same as nodesBackRightEdge() but without corners.
Definition: Mesh.h:983
auto nodesBackBottomRightCorner() const
Back-Bottom-Right corner node.
Definition: Mesh.h:1181
auto nodesBottomFrontRightCorner() const
Alias of nodesFrontBottomRightCorner().
Definition: Mesh.h:1262
auto nodesBottomRightBackCorner() const
Alias of nodesBackBottomRightCorner().
Definition: Mesh.h:1451
auto nodesLeftTopOpenEdge() const
Alias of nodesTopLeftOpenEdge().
Definition: Mesh.h:1091
auto nodesBottomFrontLeftCorner() const
Alias of nodesFrontBottomLeftCorner().
Definition: Mesh.h:1217
D derived_type
Underlying type.
Definition: Mesh.h:557
auto nodesBottomFrontEdge() const
Alias of nodesFrontBottomEdge()
Definition: Mesh.h:746
auto nodesBottomLeftFrontCorner() const
Alias of nodesFrontBottomLeftCorner().
Definition: Mesh.h:1226
auto nodesFront() const
Nodes along the front face (z = 0).
Definition: Mesh.h:608
auto nodesBackBottomEdge() const
Nodes along the edge at the intersection of the back and bottom faces (z = nelz * h and y = nely * h)...
Definition: Mesh.h:667
auto nodesBottomFrontOpenEdge() const
Alias of nodesFrontBottomOpenEdge().
Definition: Mesh.h:1028
auto nodesBackBottomOpenEdge() const
Same as nodesBackBottomEdge() but without corners.
Definition: Mesh.h:956
auto nodesBackFace() const
Nodes along the back face excluding edges.
Definition: Mesh.h:867
auto nodesBottomBackEdge() const
Alias of nodesBackBottomEdge()
Definition: Mesh.h:755
auto nodesBackLeftOpenEdge() const
Same as nodesBackLeftEdge() but without corners.
Definition: Mesh.h:974
auto nodesBottomBackRightCorner() const
Alias of nodesBackBottomRightCorner().
Definition: Mesh.h:1442
auto nodesRightBackEdge() const
Alias of nodesBackRightEdge()
Definition: Mesh.h:845
auto nodesBottomFace() const
Nodes along the bottom face excluding edges.
Definition: Mesh.h:900
auto nodesLeftBackOpenEdge() const
Alias of nodesBackLeftOpenEdge().
Definition: Mesh.h:1082
auto nodesFrontLeftEdge() const
Nodes along the edge at the intersection of the front and left faces (z = 0 and x = 0).
Definition: Mesh.h:647
auto nodesTopFrontRightCorner() const
Alias of nodesFrontTopRightCorner().
Definition: Mesh.h:1352
auto nodesFrontBottomLeftCorner() const
Front-Bottom-Left corner node.
Definition: Mesh.h:1136
auto nodesRight() const
Nodes along the right face (x = nelx * h).
Definition: Mesh.h:599
auto nodesLeftTopBackCorner() const
Alias of nodesBackTopLeftCorner().
Definition: Mesh.h:1514
auto nodesRightBottomBackCorner() const
Alias of nodesBackBottomRightCorner().
Definition: Mesh.h:1469
auto nodesBottom() const
Nodes along the bottom face (y = 0).
Definition: Mesh.h:572
auto nodesFrontRightBottomCorner() const
Alias of nodesFrontBottomRightCorner().
Definition: Mesh.h:1253
auto nodesLeftBackEdge() const
Alias of nodesBackLeftEdge()
Definition: Mesh.h:800
auto nodesLeftBackBottomCorner() const
Alias of nodesBackBottomLeftCorner().
Definition: Mesh.h:1415
auto nodesRightFrontOpenEdge() const
Alias of nodesFrontRightOpenEdge().
Definition: Mesh.h:1118
auto nodesTopBackEdge() const
Alias of nodesBackTopEdge()
Definition: Mesh.h:773
auto nodesBackTopRightCorner() const
Back-Top-Right corner node.
Definition: Mesh.h:1199
auto nodesRightFrontBottomCorner() const
Alias of nodesFrontBottomRightCorner().
Definition: Mesh.h:1280
auto nodesRightFrontTopCorner() const
Alias of nodesFrontTopRightCorner().
Definition: Mesh.h:1370
auto nodesTop() const
Nodes along the top face (y = nely * h).
Definition: Mesh.h:581
CRTP base class for regular meshes.
Definition: Mesh.h:181
auto ndim() const
Number of dimensions == 2.
Definition: Mesh.h:219
D derived_type
Underlying type.
Definition: Mesh.h:186
auto nely() const
Number of elements in y-direction == height of the mesh, in units of h,.
Definition: Mesh.h:237
auto nnode() const
Number of nodes.
Definition: Mesh.h:201
auto nodesPeriodic() const
Periodic node pairs, in two columns: (independent, dependent).
Definition: Mesh.h:309
auto dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
Definition: Mesh.h:282
auto coor() const
Nodal coordinates [nnode, ndim].
Definition: Mesh.h:264
auto dofsPeriodic() const
DOF-numbers for the case that the periodicity if fully eliminated.
Definition: Mesh.h:291
auto getElementType() const
The ElementType().
Definition: Mesh.h:255
auto nodesOrigin() const
Reference node to use for periodicity, because all corners are tied to it.
Definition: Mesh.h:318
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 nne() const
Number of nodes-per-element == 4.
Definition: Mesh.h:210
auto nelem() const
Number of elements.
Definition: Mesh.h:192
auto conn() const
Connectivity [nelem, nne].
Definition: Mesh.h:273
Renumber indices to lowest possible index.
Definition: Mesh.h:114
Renumber(const T &dofs)
Definition: Mesh.h:122
T apply(const T &list) const
Apply renumbering to other set.
Definition: Mesh.h:144
const array_type::tensor< size_t, 1 > & index() const
Get the list needed to renumber, e.g.:
Definition: Mesh.h:156
Reorder to lowest possible index, in specific order.
Definition: Mesh.h:2387
T apply(const T &list) const
Apply reordering to other set.
Definition: Mesh.h:2430
const array_type::tensor< size_t, 1 > & index() const
Get the list needed to reorder, e.g.:
Definition: Mesh.h:2450
Reorder(const std::initializer_list< T > args)
Definition: Mesh.h:2395
Stitch mesh objects, automatically searching for overlapping nodes.
Definition: Mesh.h:2010
T nodeset(const std::vector< T > &set) const
Combine set of node numbers for an original to the final mesh (removes duplicates).
Definition: Mesh.h:2223
size_t ndim() const
Number of dimensions.
Definition: Mesh.h:2102
T nodeset(std::initializer_list< T > set) const
Combine set of node numbers for an original to the final mesh (removes duplicates).
Definition: Mesh.h:2249
std::vector< size_t > m_nel
Number of elements per sub-mesh.
Definition: Mesh.h:2296
std::vector< array_type::tensor< size_t, 1 > > m_map
See nodemap(size_t)
Definition: Mesh.h:2295
size_t nne() const
Number of nodes-per-element.
Definition: Mesh.h:2093
const array_type::tensor< size_t, 2 > & conn() const
Connectivity [nelem, nne].
Definition: Mesh.h:2120
array_type::tensor< size_t, 2 > dofs() const
DOF numbers for each node (numbered sequentially) [nnode, ndim].
Definition: Mesh.h:2129
size_t nelem() const
Number of elements.
Definition: Mesh.h:2075
double m_atol
Absolute tolerance to find overlapping nodes.
Definition: Mesh.h:2299
const array_type::tensor< double, 2 > & coor() const
Nodal coordinates [nnode, ndim].
Definition: Mesh.h:2111
T elemset(const std::vector< T > &set) const
Combine set of element numbers for an original to the final mesh.
Definition: Mesh.h:2261
std::vector< array_type::tensor< size_t, 1 > > nodemap() const
Node-map per sub-mesh.
Definition: Mesh.h:2140
T elemset(std::initializer_list< T > set) const
Combine set of element numbers for an original to the final mesh.
Definition: Mesh.h:2287
size_t nmesh() const
Number of sub meshes.
Definition: Mesh.h:2066
std::vector< size_t > m_el_offset
First element of every sub-mesh.
Definition: Mesh.h:2297
T nodeset(const T &set, size_t mesh_index) const
Convert set of node-numbers for a sub-mesh to the stitched mesh.
Definition: Mesh.h:2194
double m_rtol
Relative tolerance to find overlapping nodes.
Definition: Mesh.h:2298
Stitch(double rtol=1e-5, double atol=1e-8)
Definition: Mesh.h:2016
void push_back(const C &coor, const E &conn)
Add mesh to be stitched.
Definition: Mesh.h:2029
T elemset(const T &set, size_t mesh_index) const
Convert set of element-numbers for a sub-mesh to the stitched mesh.
Definition: Mesh.h:2209
array_type::tensor< size_t, 2 > m_conn
Connectivity [nelem, nne].
Definition: Mesh.h:2294
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.
Definition: Mesh.h:2180
size_t nnode() const
Number of nodes.
Definition: Mesh.h:2084
std::vector< array_type::tensor< size_t, 1 > > elemmap() const
Element-map per sub-mesh.
Definition: Mesh.h:2153
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.
Definition: Mesh.h:2168
array_type::tensor< double, 2 > m_coor
Nodal coordinates [nnode, ndim].
Definition: Mesh.h:2293
Vertically stack meshes.
Definition: Mesh.h:2305
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.
Definition: Mesh.h:2329
Vstack(bool check_overlap=true, double rtol=1e-5, double atol=1e-8)
Definition: Mesh.h:2312
Class to switch between storage types.
Definition: Vector.h:23
array_type::tensor< double, 2 > AsNode(const T &arg) const
Convert "dofval" or "elemvec" to "nodevec" (overwrite entries that occur more than once).
Definition: Vector.h:168
array_type::tensor< double, 3 > AsElement(const T &arg) const
Convert "dofval" or "nodevec" to "elemvec" (overwrite entries that occur more than once).
Definition: Vector.h:194
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
Definition: config.h:96
#define GOOSEFEM_CHECK(expr)
Assertion that cannot be switched off.
Definition: config.h:106
array_type::tensor< double, 1 > w()
Integration point weights.
Definition: ElementQuad4.h:149
array_type::tensor< double, 2 > xi()
Integration point coordinates (local coordinates).
Definition: ElementQuad4.h:122
array_type::tensor< size_t, 1 > coordination(const E &conn)
Number of elements connected to each node.
Definition: Mesh.h:2466
array_type::tensor< size_t, 2 > dofs(size_t nnode, size_t ndim)
List with DOF-numbers in sequential order.
Definition: Mesh.h:93
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.
Definition: Mesh.h:1716
array_type::tensor< double, 2 > nodal_mass(const C &coor, const E &conn, ElementType type)
Compute the mass of each node in the mesh.
Definition: Mesh.h:2735
std::vector< std::vector< size_t > > node2dof(const D &dofs, bool sorted=true)
Nodes connected to each DOF.
Definition: Mesh.h:2560
array_type::tensor< double, 2 > edgesize(const C &coor, const E &conn, ElementType type)
Return size of each element edge.
Definition: Mesh.h:2574
T renumber(const T &dofs)
Renumber to lowest possible index (see GooseFEM::Mesh::Renumber).
Definition: Mesh.h:172
array_type::tensor< double, 2 > centers(const C &coor, const E &conn, ElementType type)
Coordinates of the center of each element.
Definition: Mesh.h:2629
array_type::tensor< double, 1 > center_of_gravity(const C &coor, const E &conn, ElementType type)
Compute the center of gravity of a mesh.
Definition: Mesh.h:2811
ElementType
Enumerator for element-types.
Definition: Mesh.h:31
@ 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.
Definition: Mesh.h:2489
std::vector< std::vector< size_t > > nodaltyings(const D &dofs)
List nodal tyings based on DOF-numbers per node.
Definition: Mesh.h:2842
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.
Definition: Mesh.h:2675
ElementType defaultElementType(const S &coor, const T &conn)
Extract the element type based on the connectivity.
Definition: Mesh.h:46
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:176
Toolbox to perform finite element computations.
Definition: Allocate.h:14
bool is_unique(const T &arg)
Returns true is a list is unique (has not duplicate items).
Definition: assertions.h:20