18#include <boost/math/constants/constants.hpp>
21#ifdef ENABLE_PYTHON_BINDING
22 #include <boost/python.hpp>
23 #include <boost/python/def.hpp>
24 #include <boost/python/numpy.hpp>
25namespace bp = boost::python;
26namespace np = boost::python::numpy;
28 #pragma message "With ENABLE_PYTHON_BINDING"
32 #pragma message "Without ENABLE_PYTHON_BINDING"
58 const EntityType type) {
62 auto dim = CN::Dimension(type);
64 std::vector<int> sendcounts(pcomm->size());
65 std::vector<int> displs(pcomm->size());
66 std::vector<int> sendbuf(r.size());
67 if (pcomm->rank() == 0) {
68 for (
auto p = 1; p != pcomm->size(); p++) {
70 ->getPartEntities(m_field.
get_moab(), p)
73 CHKERR m_field.
get_moab().get_adjacencies(part_ents, dim,
true, faces,
74 moab::Interface::UNION);
75 faces = intersect(faces, r);
76 sendcounts[p] = faces.size();
77 displs[p] = sendbuf.size();
78 for (
auto f : faces) {
80 sendbuf.push_back(
id);
86 MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
88 std::vector<int> recvbuf(recv_data);
89 MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
90 recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
92 if (pcomm->rank() > 0) {
94 for (
auto &f : recvbuf) {
104 const std::string block_name,
int dim) {
110 std::regex((boost::format(
"%s(.*)") % block_name).str())
114 for (
auto bc : bcs) {
126 const std::string block_name,
int dim) {
127 std::map<std::string, Range> r;
132 std::regex((boost::format(
"%s(.*)") % block_name).str())
136 for (
auto bc : bcs) {
141 r[bc->getName()] = faces;
148 const unsigned int cubit_bc_type) {
151 CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
155static auto save_range(moab::Interface &moab,
const std::string name,
156 const Range r, std::vector<Tag> tags = {}) {
159 CHKERR moab.add_entities(*out_meshset, r);
161 CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1,
162 tags.data(), tags.size());
164 MOFEM_LOG(
"SELF", Sev::warning) <<
"Empty range for " << name;
171 ParallelComm *pcomm =
174 PSTATUS_SHARED | PSTATUS_MULTISHARED,
175 PSTATUS_NOT, -1, &boundary_ents),
177 return boundary_ents;
182 ParallelComm *pcomm =
184 CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
193 CHK_MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents),
"find_skin");
199 ParallelComm *pcomm =
202 Range crack_skin_without_bdy;
203 if (pcomm->rank() == 0) {
205 CHKERR moab.get_adjacencies(crack_faces, 1,
true, crack_edges,
206 moab::Interface::UNION);
207 auto crack_skin =
get_skin(m_field, crack_faces);
211 "get_entities_by_dimension");
212 auto body_skin =
get_skin(m_field, body_ents);
213 Range body_skin_edges;
214 CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1,
true, body_skin_edges,
215 moab::Interface::UNION),
217 crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
219 for (
auto &
m : front_edges_map) {
220 auto add_front = subtract(
m.second, crack_edges);
221 auto i = intersect(
m.second, crack_edges);
223 crack_skin_without_bdy.merge(add_front);
227 CHKERR moab.get_adjacencies(i_skin, 1,
true, adj_i_skin,
228 moab::Interface::UNION);
229 adj_i_skin = subtract(intersect(adj_i_skin,
m.second), crack_edges);
230 crack_skin_without_bdy.merge(adj_i_skin);
234 return send_type(m_field, crack_skin_without_bdy, MBEDGE);
240 ParallelComm *pcomm =
243 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface";
245 if (!pcomm->rank()) {
247 auto impl = [&](
auto &saids) {
252 auto get_adj = [&](
auto e,
auto dim) {
255 e, dim,
true, adj, moab::Interface::UNION),
260 auto get_conn = [&](
auto e) {
267 constexpr bool debug =
false;
271 auto body_skin =
get_skin(m_field, body_ents);
272 auto body_skin_edges = get_adj(body_skin, 1);
275 subtract(
get_skin(m_field, crack_faces), body_skin_edges);
276 auto crack_skin_conn = get_conn(crack_skin);
277 auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
278 auto crack_edges = get_adj(crack_faces, 1);
279 crack_edges = subtract(crack_edges, crack_skin);
280 auto all_tets = get_adj(crack_edges, 3);
281 crack_edges = subtract(crack_edges, crack_skin_conn_edges);
282 auto crack_conn = get_conn(crack_edges);
283 all_tets.merge(get_adj(crack_conn, 3));
292 if (crack_faces.size()) {
293 auto grow = [&](
auto r) {
294 auto crack_faces_conn = get_conn(crack_faces);
297 while (size_r != r.size() && r.size() > 0) {
299 CHKERR moab.get_connectivity(r,
v,
true);
300 v = subtract(
v, crack_faces_conn);
303 moab::Interface::UNION);
304 r = intersect(r, all_tets);
313 Range all_tets_ord = all_tets;
314 while (all_tets.size()) {
315 Range faces = get_adj(unite(saids.first, saids.second), 2);
316 faces = subtract(crack_faces, faces);
319 auto fit = faces.begin();
320 for (; fit != faces.end(); ++fit) {
321 tets = intersect(get_adj(
Range(*fit, *fit), 3), all_tets);
322 if (tets.size() == 2) {
329 saids.first.insert(tets[0]);
330 saids.first = grow(saids.first);
331 all_tets = subtract(all_tets, saids.first);
332 if (tets.size() == 2) {
333 saids.second.insert(tets[1]);
334 saids.second = grow(saids.second);
335 all_tets = subtract(all_tets, saids.second);
343 saids.first = subtract(all_tets_ord, saids.second);
344 saids.second = subtract(all_tets_ord, saids.first);
350 std::pair<Range, Range> saids;
351 if (crack_faces.size())
356 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface <- done";
358 return std::pair<Range, Range>();
366 boost::shared_ptr<Range> front_nodes,
367 boost::shared_ptr<Range> front_edges,
368 boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi =
nullptr)
372 int order_col,
int order_data) {
375 constexpr bool debug =
false;
377 constexpr int numNodes = 4;
378 constexpr int numEdges = 6;
379 constexpr int refinementLevels = 6;
381 auto &m_field = fe_raw_ptr->
mField;
382 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
385 auto set_base_quadrature = [&]() {
387 int rule = 2 * order_data + 1;
398 auto &gauss_pts = fe_ptr->gaussPts;
399 gauss_pts.resize(4, nb_gauss_pts,
false);
400 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
401 &gauss_pts(0, 0), 1);
402 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
403 &gauss_pts(1, 0), 1);
404 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
405 &gauss_pts(2, 0), 1);
407 &gauss_pts(3, 0), 1);
408 auto &data = fe_ptr->dataOnElement[
H1];
409 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
412 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
413 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
422 CHKERR set_base_quadrature();
426 auto get_singular_nodes = [&]() {
429 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
431 std::bitset<numNodes> singular_nodes;
432 for (
auto nn = 0; nn != numNodes; ++nn) {
434 singular_nodes.set(nn);
436 singular_nodes.reset(nn);
439 return singular_nodes;
442 auto get_singular_edges = [&]() {
443 std::bitset<numEdges> singular_edges;
444 for (
int ee = 0; ee != numEdges; ee++) {
446 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
448 singular_edges.set(ee);
450 singular_edges.reset(ee);
453 return singular_edges;
456 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
458 fe_ptr->gaussPts.swap(ref_gauss_pts);
459 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
460 auto &data = fe_ptr->dataOnElement[
H1];
461 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
463 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
465 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
470 auto singular_nodes = get_singular_nodes();
471 if (singular_nodes.count()) {
472 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
474 CHKERR set_gauss_pts(it_map_ref_coords->second);
478 auto refine_quadrature = [&]() {
481 const int max_level = refinementLevels;
485 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
487 for (
int nn = 0; nn != 4; nn++)
488 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
489 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
493 Range tets(tet, tet);
496 tets, 1,
true, edges, moab::Interface::UNION);
501 Range nodes_at_front;
502 for (
int nn = 0; nn != numNodes; nn++) {
503 if (singular_nodes[nn]) {
505 CHKERR moab_ref.side_element(tet, 0, nn, ent);
506 nodes_at_front.insert(ent);
510 auto singular_edges = get_singular_edges();
513 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
514 for (
int ee = 0; ee != numEdges; ee++) {
515 if (singular_edges[ee]) {
517 CHKERR moab_ref.side_element(tet, 1, ee, ent);
518 CHKERR moab_ref.add_entities(meshset, &ent, 1);
524 for (
int ll = 0; ll != max_level; ll++) {
527 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
531 CHKERR moab_ref.get_adjacencies(
532 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
533 ref_edges = intersect(ref_edges, edges);
535 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
536 ref_edges = intersect(ref_edges, ents);
539 ->getEntitiesByTypeAndRefLevel(
541 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
545 ->updateMeshsetByEntitiesChildren(meshset,
547 meshset, MBEDGE,
true);
553 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
563 for (Range::iterator tit = tets.begin(); tit != tets.end();
567 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
568 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
571 auto &data = fe_ptr->dataOnElement[
H1];
572 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
573 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
575 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
577 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
578 double *tet_coords = &ref_coords(tt, 0);
581 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
582 for (
int dd = 0; dd != 3; dd++) {
583 ref_gauss_pts(dd, gg) =
584 shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
585 shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
586 shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
587 shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
589 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
593 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
600 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
601 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
606 CHKERR refine_quadrature();
616 using ForcesAndSourcesCore::dataOnElement;
619 using ForcesAndSourcesCore::ForcesAndSourcesCore;
625 boost::shared_ptr<CGGUserPolynomialBase::CachePhi>
cachePhi;
636 boost::shared_ptr<Range> front_edges)
640 boost::shared_ptr<Range> front_edges,
645 int order_col,
int order_data) {
648 constexpr bool debug =
false;
650 constexpr int numNodes = 3;
651 constexpr int numEdges = 3;
652 constexpr int refinementLevels = 6;
654 auto &m_field = fe_raw_ptr->
mField;
655 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
658 auto set_base_quadrature = [&]() {
660 int rule =
funRule(order_data);
670 fe_ptr->gaussPts.resize(3, nb_gauss_pts,
false);
671 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
672 &fe_ptr->gaussPts(0, 0), 1);
673 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
674 &fe_ptr->gaussPts(1, 0), 1);
676 &fe_ptr->gaussPts(2, 0), 1);
677 auto &data = fe_ptr->dataOnElement[
H1];
678 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
681 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
682 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
684 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
687 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
696 CHKERR set_base_quadrature();
700 auto get_singular_nodes = [&]() {
703 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
705 std::bitset<numNodes> singular_nodes;
706 for (
auto nn = 0; nn != numNodes; ++nn) {
708 singular_nodes.set(nn);
710 singular_nodes.reset(nn);
713 return singular_nodes;
716 auto get_singular_edges = [&]() {
717 std::bitset<numEdges> singular_edges;
718 for (
int ee = 0; ee != numEdges; ee++) {
720 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
722 singular_edges.set(ee);
724 singular_edges.reset(ee);
727 return singular_edges;
730 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
732 fe_ptr->gaussPts.swap(ref_gauss_pts);
733 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
734 auto &data = fe_ptr->dataOnElement[
H1];
735 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
737 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
739 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
743 auto singular_nodes = get_singular_nodes();
744 if (singular_nodes.count()) {
745 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
747 CHKERR set_gauss_pts(it_map_ref_coords->second);
751 auto refine_quadrature = [&]() {
754 const int max_level = refinementLevels;
757 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
759 for (
int nn = 0; nn != numNodes; nn++)
760 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
762 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
766 Range tris(tri, tri);
769 tris, 1,
true, edges, moab::Interface::UNION);
774 Range nodes_at_front;
775 for (
int nn = 0; nn != numNodes; nn++) {
776 if (singular_nodes[nn]) {
778 CHKERR moab_ref.side_element(tri, 0, nn, ent);
779 nodes_at_front.insert(ent);
783 auto singular_edges = get_singular_edges();
786 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
787 for (
int ee = 0; ee != numEdges; ee++) {
788 if (singular_edges[ee]) {
790 CHKERR moab_ref.side_element(tri, 1, ee, ent);
791 CHKERR moab_ref.add_entities(meshset, &ent, 1);
797 for (
int ll = 0; ll != max_level; ll++) {
800 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
804 CHKERR moab_ref.get_adjacencies(
805 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
806 ref_edges = intersect(ref_edges, edges);
808 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
809 ref_edges = intersect(ref_edges, ents);
812 ->getEntitiesByTypeAndRefLevel(
814 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
818 ->updateMeshsetByEntitiesChildren(meshset,
820 meshset, MBEDGE,
true);
826 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
837 for (Range::iterator tit = tris.begin(); tit != tris.end();
841 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
842 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
845 auto &data = fe_ptr->dataOnElement[
H1];
846 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
847 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
849 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
851 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
852 double *tri_coords = &ref_coords(tt, 0);
855 auto det = t_normal.
l2();
856 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
857 for (
int dd = 0; dd != 2; dd++) {
858 ref_gauss_pts(dd, gg) =
859 shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
860 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
861 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
863 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
867 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
873 CHKERR refine_quadrature();
883 using ForcesAndSourcesCore::dataOnElement;
886 using ForcesAndSourcesCore::ForcesAndSourcesCore;
935 const char *list_rots[] = {
"small",
"moderate",
"large",
"no_h1"};
936 const char *list_symm[] = {
"symm",
"not_symm"};
937 const char *list_release[] = {
"griffith_force",
"griffith_skeleton"};
938 const char *list_stretches[] = {
"linear",
"log",
"log_quadratic"};
943 PetscInt choice_stretch = StretchSelector::LOG;
944 char analytical_expr_file_name[255] =
"analytical_expr.py";
946 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
948 CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
950 CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
952 CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
954 CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
956 CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
958 CHKERR PetscOptionsScalar(
"-viscosity_alpha_omega",
"rot viscosity",
"",
960 CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
964 CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
965 LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
967 CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
968 list_rots, NO_H1_CONFIGURATION + 1,
969 list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
970 CHKERR PetscOptionsEList(
"-symm",
"symmetric variant",
"", list_symm, 2,
971 list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
975 CHKERR PetscOptionsEList(
"-stretches",
"stretches",
"", list_stretches,
976 StretchSelector::STRETCH_SELECTOR_LAST,
977 list_stretches[choice_stretch], &choice_stretch,
980 CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
982 CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
984 CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
988 CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
992 CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
999 CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
1001 CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
1003 CHKERR PetscOptionsEList(
"-energy_release_variant",
"energy release variant",
1004 "", list_release, 2, list_release[choice_release],
1005 &choice_release, PETSC_NULLPTR);
1006 CHKERR PetscOptionsInt(
"-nb_J_integral_levels",
"Number of J integarl levels",
1010 "-nb_J_integral_contours",
"Number of J integral contours",
"",
1014 char tag_name[255] =
"";
1015 CHKERR PetscOptionsString(
"-internal_stress_tag_name",
1016 "internal stress tag name",
"",
"", tag_name, 255,
1020 "-internal_stress_interp_order",
"internal stress interpolation order",
1022 CHKERR PetscOptionsBool(
"-internal_stress_voigt",
"Voigt index notation",
"",
1027 analytical_expr_file_name, 255, PETSC_NULLPTR);
1033 "Unsupported internal stress interpolation order %d",
1046 static_cast<EnergyReleaseSelector
>(choice_release);
1049 case StretchSelector::LINEAR:
1057 case StretchSelector::LOG:
1059 std::numeric_limits<float>::epsilon()) {
1075 case StretchSelector::LOG_QUADRATIC:
1081 "No logarithmic quadratic stretch for this case");
1097 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaU: -viscosity_alpha_u " <<
alphaU;
1098 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaW: -viscosity_alpha_w " <<
alphaW;
1100 <<
"alphaOmega: -viscosity_alpha_omega " <<
alphaOmega;
1105 MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
1113 MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
1115 <<
"Stretch: -stretches " << list_stretches[choice_stretch];
1141 <<
"Energy release variant: -energy_release_variant "
1144 <<
"Number of J integral contours: -nb_J_integral_contours "
1147#ifdef ENABLE_PYTHON_BINDING
1148 auto file_exists = [](std::string myfile) {
1149 std::ifstream file(myfile.c_str());
1156 if (file_exists(analytical_expr_file_name)) {
1157 MOFEM_LOG(
"EP", Sev::inform) << analytical_expr_file_name <<
" file found";
1161 analytical_expr_file_name);
1165 << analytical_expr_file_name <<
" file NOT found";
1178 auto get_tets = [&]() {
1184 auto get_tets_skin = [&]() {
1185 Range tets_skin_part;
1187 CHKERR skin.find_skin(0, get_tets(),
false, tets_skin_part);
1188 ParallelComm *pcomm =
1191 CHKERR pcomm->filter_pstatus(tets_skin_part,
1192 PSTATUS_SHARED | PSTATUS_MULTISHARED,
1193 PSTATUS_NOT, -1, &tets_skin);
1197 auto subtract_boundary_conditions = [&](
auto &&tets_skin) {
1203 tets_skin = subtract(tets_skin,
v.faces);
1207 tets_skin = subtract(tets_skin,
v.faces);
1212 tets_skin = subtract(tets_skin,
v.faces);
1218 auto add_blockset = [&](
auto block_name,
auto &&tets_skin) {
1221 tets_skin.merge(crack_faces);
1225 auto subtract_blockset = [&](
auto block_name,
auto &&tets_skin) {
1226 auto contact_range =
1228 tets_skin = subtract(tets_skin, contact_range);
1232 auto get_stress_trace_faces = [&](
auto &&tets_skin) {
1235 faces, moab::Interface::UNION);
1236 Range trace_faces = subtract(faces, tets_skin);
1240 auto tets = get_tets();
1244 auto trace_faces = get_stress_trace_faces(
1246 subtract_blockset(
"CONTACT",
1247 subtract_boundary_conditions(get_tets_skin()))
1254 boost::make_shared<Range>(subtract(trace_faces, *
contactFaces));
1269 auto add_broken_hdiv_field = [
this, meshset](
const std::string
field_name,
1275 auto get_side_map_hdiv = [&]() {
1278 std::pair<EntityType,
1293 get_side_map_hdiv(), MB_TAG_DENSE,
MF_ZERO);
1299 auto add_l2_field = [
this, meshset](
const std::string
field_name,
1300 const int order,
const int dim) {
1309 auto add_h1_field = [
this, meshset](
const std::string
field_name,
1310 const int order,
const int dim) {
1322 auto add_l2_field_by_range = [
this](
const std::string
field_name,
1323 const int order,
const int dim,
1324 const int field_dim,
Range &&r) {
1334 auto add_bubble_field = [
this, meshset](
const std::string
field_name,
1335 const int order,
const int dim) {
1341 auto field_order_table =
1342 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1343 auto get_cgg_bubble_order_zero = [](
int p) {
return 0; };
1344 auto get_cgg_bubble_order_tet = [](
int p) {
1347 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
1348 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
1349 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
1350 field_order_table[MBTET] = get_cgg_bubble_order_tet;
1357 auto add_user_l2_field = [
this, meshset](
const std::string
field_name,
1358 const int order,
const int dim) {
1364 auto field_order_table =
1365 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1366 auto zero_dofs = [](
int p) {
return 0; };
1368 field_order_table[MBVERTEX] = zero_dofs;
1369 field_order_table[MBEDGE] = zero_dofs;
1370 field_order_table[MBTRI] = zero_dofs;
1371 field_order_table[MBTET] = dof_l2_tet;
1389 auto get_hybridised_disp = [&]() {
1391 auto skin = subtract_boundary_conditions(get_tets_skin());
1393 faces.merge(intersect(bc.faces, skin));
1399 get_hybridised_disp());
1425 auto project_ho_geometry = [&](
auto field) {
1431 auto get_adj_front_edges = [&](
auto &front_edges) {
1432 Range front_crack_nodes;
1433 Range crack_front_edges_with_both_nodes_not_at_front;
1438 moab.get_connectivity(front_edges, front_crack_nodes,
true),
1439 "get_connectivity failed");
1440 Range crack_front_edges;
1442 false, crack_front_edges,
1443 moab::Interface::UNION),
1444 "get_adjacencies failed");
1445 Range crack_front_edges_nodes;
1447 crack_front_edges_nodes,
true),
1448 "get_connectivity failed");
1450 crack_front_edges_nodes =
1451 subtract(crack_front_edges_nodes, front_crack_nodes);
1452 Range crack_front_edges_with_both_nodes_not_at_front;
1454 moab.get_adjacencies(crack_front_edges_nodes, 1,
false,
1455 crack_front_edges_with_both_nodes_not_at_front,
1456 moab::Interface::UNION),
1457 "get_adjacencies failed");
1459 crack_front_edges_with_both_nodes_not_at_front =
1460 intersect(crack_front_edges,
1461 crack_front_edges_with_both_nodes_not_at_front);
1465 crack_front_edges_with_both_nodes_not_at_front =
send_type(
1466 mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
1468 return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
1469 boost::make_shared<Range>(
1470 crack_front_edges_with_both_nodes_not_at_front));
1477 auto [front_vertices, front_adj_edges] = get_adj_front_edges(*
frontEdges);
1482 <<
"Number of crack faces: " <<
crackFaces->size();
1484 <<
"Number of front edges: " <<
frontEdges->size();
1488 <<
"Number of front adjacent edges: " <<
frontAdjEdges->size();
1497 (boost::format(
"crack_faces_%d.vtk") % rank).str(),
1500 (boost::format(
"front_edges_%d.vtk") % rank).str(),
1511 auto set_singular_dofs = [&](
auto &front_adj_edges,
auto &front_vertices) {
1519 MOFEM_LOG(
"EP", Sev::inform) <<
"Singularity eps " << beta;
1524 [&](boost::shared_ptr<FieldEntity> field_entity_ptr) ->
MoFEMErrorCode {
1529 auto nb_dofs = field_entity_ptr->getEntFieldData().size();
1535 if (field_entity_ptr->getNbOfCoeffs() != 3)
1537 "Expected 3 coefficients per edge");
1538 if (nb_dofs % 3 != 0)
1540 "Expected multiple of 3 coefficients per edge");
1543 auto get_conn = [&]() {
1546 CHKERR moab.get_connectivity(field_entity_ptr->getEnt(), conn,
1548 return std::make_pair(conn, num_nodes);
1551 auto get_dir = [&](
auto &&conn_p) {
1552 auto [conn, num_nodes] = conn_p;
1554 CHKERR moab.get_coords(conn, num_nodes, coords);
1556 coords[4] - coords[1],
1557 coords[5] - coords[2]};
1561 auto get_singularity_dof = [&](
auto &&conn_p,
auto &&t_edge_dir) {
1562 auto [conn, num_nodes] = conn_p;
1564 if (front_vertices.find(conn[0]) != front_vertices.end()) {
1565 t_singularity_dof(
i) = t_edge_dir(
i) * (-
eps);
1566 }
else if (front_vertices.find(conn[1]) != front_vertices.end()) {
1567 t_singularity_dof(
i) = t_edge_dir(
i) *
eps;
1569 return t_singularity_dof;
1572 auto t_singularity_dof =
1573 get_singularity_dof(get_conn(), get_dir(get_conn()));
1575 auto field_data = field_entity_ptr->getEntFieldData();
1577 &field_data[0], &field_data[1], &field_data[2]};
1579 t_dof(
i) = t_singularity_dof(
i);
1581 for (
auto n = 1;
n < field_data.size() / 3; ++
n) {
1606 auto add_field_to_fe = [
this](
const std::string fe,
1645 auto set_fe_adjacency = [&](
auto fe_name) {
1648 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
1656 auto add_field_to_fe = [
this](
const std::string fe,
1668 Range natural_bc_elements;
1671 natural_bc_elements.merge(
v.faces);
1676 natural_bc_elements.merge(
v.faces);
1681 natural_bc_elements.merge(
v.faces);
1686 natural_bc_elements.merge(
v.faces);
1691 natural_bc_elements.merge(
v.faces);
1696 natural_bc_elements.merge(
v.faces);
1701 natural_bc_elements.merge(
v.faces);
1704 natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
1715 auto get_skin = [&](
auto &body_ents) {
1718 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
1723 Range boundary_ents;
1724 ParallelComm *pcomm =
1726 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
1727 PSTATUS_NOT, -1, &boundary_ents);
1728 return boundary_ents;
1813 auto remove_dofs_on_broken_skin = [&](
const std::string prb_name) {
1815 for (
int d : {0, 1, 2}) {
1816 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1818 ->getSideDofsOnBrokenSpaceEntities(
1829 CHKERR remove_dofs_on_broken_skin(
"ESHELBY_PLASTICITY");
1865 auto set_zero_block = [&]() {
1897 auto set_section = [&]() {
1899 PetscSection section;
1904 CHKERR PetscSectionDestroy(§ion);
1927BcDisp::BcDisp(std::string name, std::vector<double> attr,
Range faces)
1928 : blockName(name), faces(faces) {
1929 vals.resize(3,
false);
1930 flags.resize(3,
false);
1931 for (
int ii = 0; ii != 3; ++ii) {
1932 vals[ii] = attr[ii];
1933 flags[ii] =
static_cast<int>(attr[ii + 3]);
1936 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
1938 <<
"Add BCDisp vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1940 <<
"Add BCDisp flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1941 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " <<
faces.size();
1945 : blockName(name), faces(faces) {
1946 vals.resize(attr.size(),
false);
1947 for (
int ii = 0; ii != attr.size(); ++ii) {
1948 vals[ii] = attr[ii];
1954 : blockName(name), faces(faces) {
1955 vals.resize(3,
false);
1956 flags.resize(3,
false);
1957 for (
int ii = 0; ii != 3; ++ii) {
1958 vals[ii] = attr[ii];
1959 flags[ii] =
static_cast<int>(attr[ii + 3]);
1962 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
1964 <<
"Add BCForce vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1966 <<
"Add BCForce flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1967 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " <<
faces.size();
1971 std::vector<double> attr,
1973 : blockName(name), faces(faces) {
1976 if (attr.size() < 1) {
1978 "Wrong size of normal displacement BC");
1983 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc " << name;
1984 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc val " <<
val;
1986 <<
"Add NormalDisplacementBc nb. of faces " <<
faces.size();
1990 : blockName(name), faces(faces) {
1993 if (attr.size() < 1) {
1995 "Wrong size of normal displacement BC");
2000 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc " << name;
2001 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc val " <<
val;
2003 <<
"Add PressureBc nb. of faces " <<
faces.size();
2009 : blockName(name), ents(ents) {
2012 if (attr.size() < 2) {
2014 "Wrong size of external strain attribute");
2020 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain " << name;
2021 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain val " <<
val;
2022 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain bulk modulus K "
2025 <<
"Add ExternalStrain nb. of tets " <<
ents.size();
2031 std::vector<double> attr,
2033 : blockName(name), faces(faces) {
2036 if (attr.size() < 3) {
2038 "Wrong size of analytical displacement BC");
2041 flags.resize(3,
false);
2042 for (
int ii = 0; ii != 3; ++ii) {
2043 flags[ii] = attr[ii];
2046 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalDisplacementBc " << name;
2048 <<
"Add AnalyticalDisplacementBc flags " <<
flags[0] <<
" " <<
flags[1]
2051 <<
"Add AnalyticalDisplacementBc nb. of faces " <<
faces.size();
2055 std::vector<double> attr,
2057 : blockName(name), faces(faces) {
2060 if (attr.size() < 3) {
2062 "Wrong size of analytical traction BC");
2065 flags.resize(3,
false);
2066 for (
int ii = 0; ii != 3; ++ii) {
2067 flags[ii] = attr[ii];
2070 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalTractionBc " << name;
2072 <<
"Add AnalyticalTractionBc flags " <<
flags[0] <<
" " <<
flags[1]
2075 <<
"Add AnalyticalTractionBc nb. of faces " <<
faces.size();
2081 boost::shared_ptr<TractionFreeBc> &bc_ptr,
2082 const std::string contact_set_name) {
2087 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
2088 Range tets_skin_part;
2089 Skinner skin(&mField.get_moab());
2090 CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
2091 ParallelComm *pcomm =
2094 CHKERR pcomm->filter_pstatus(tets_skin_part,
2095 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2096 PSTATUS_NOT, -1, &tets_skin);
2099 for (
int dd = 0; dd != 3; ++dd)
2100 (*bc_ptr)[dd] = tets_skin;
2103 if (bcSpatialDispVecPtr)
2104 for (
auto &
v : *bcSpatialDispVecPtr) {
2106 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2108 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2110 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2114 if (bcSpatialRotationVecPtr)
2115 for (
auto &
v : *bcSpatialRotationVecPtr) {
2116 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2117 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2118 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2121 if (bcSpatialNormalDisplacementVecPtr)
2122 for (
auto &
v : *bcSpatialNormalDisplacementVecPtr) {
2123 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2124 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2125 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2128 if (bcSpatialAnalyticalDisplacementVecPtr)
2129 for (
auto &
v : *bcSpatialAnalyticalDisplacementVecPtr) {
2131 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2133 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2135 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2138 if (bcSpatialTractionVecPtr)
2139 for (
auto &
v : *bcSpatialTractionVecPtr) {
2140 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2141 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2142 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2145 if (bcSpatialAnalyticalTractionVecPtr)
2146 for (
auto &
v : *bcSpatialAnalyticalTractionVecPtr) {
2147 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2148 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2149 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2152 if (bcSpatialPressureVecPtr)
2153 for (
auto &
v : *bcSpatialPressureVecPtr) {
2154 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2155 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2156 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2161 std::regex((boost::format(
"%s(.*)") % contact_set_name).str()))) {
2163 CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
2165 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
2166 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
2167 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
2184 return 2 * p_data + 1;
2190 return 2 * (p_data + 1);
2197 const int tag,
const bool do_rhs,
const bool do_lhs,
const bool calc_rates,
2199 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
2203 boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0,
MatrixDouble());
2204 fe->getUserPolynomialBase() =
2205 boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
2206 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2207 fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
2210 fe->getRuleHook = [](int, int, int) {
return -1; };
2212 SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges, bubble_cache);
2218 dataAtPts->physicsPtr = physicalEquations;
2223 piolaStress, dataAtPts->getApproxPAtPts()));
2225 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
2227 piolaStress, dataAtPts->getDivPAtPts()));
2230 fe->getOpPtrVector().push_back(
2231 physicalEquations->returnOpCalculateStretchFromStress(
2232 dataAtPts, physicalEquations));
2234 fe->getOpPtrVector().push_back(
2236 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
2240 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
2242 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
2244 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
2248 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
2250 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
2255 spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
2258 fe->getOpPtrVector().push_back(
2260 stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
2261 fe->getOpPtrVector().push_back(
2263 stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
2267 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
2269 rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
2272 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2274 spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
2281 piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2284 bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2285 var_vec, MBMAXTYPE));
2287 rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
2289 piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
2291 spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
2294 fe->getOpPtrVector().push_back(
2295 physicalEquations->returnOpCalculateVarStretchFromStress(
2296 dataAtPts, physicalEquations));
2298 fe->getOpPtrVector().push_back(
2300 stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
2306 dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
2311 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
2312 tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
2319 const int tag,
const bool add_elastic,
const bool add_material,
2320 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
2321 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
2325 auto get_body_range = [
this](
auto name,
int dim) {
2326 std::map<int, Range> map;
2329 mField.getInterface<
MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2331 (boost::format(
"%s(.*)") % name).str()
2337 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2340 map[m_ptr->getMeshsetId()] = ents;
2346 constexpr bool stab_tau_dot_variant =
false;
2348 auto local_tau_sacale = boost::make_shared<double>(1.0);
2351 using BdyEleOp = BoundaryEle::UserDataOperator;
2353 set_scale_op->doWorkRhsHook =
2355 local_tau_sacale](
DataOperator *raw_op_ptr,
int side, EntityType type,
2358 auto op_ptr =
static_cast<BdyEleOp *
>(raw_op_ptr);
2359 auto h = std::sqrt(op_ptr->getFTensor1Normal().l2());
2360 *local_tau_sacale = (alphaTau /
h);
2364 auto rule_contact = [](int, int,
int o) {
return -1; };
2367 auto set_rule_contact = [refine](
2370 int order_col,
int order_data
2374 auto rule = 2 * order_data;
2380 fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2387 fe_rhs->getOpPtrVector().push_back(
2389 fe_rhs->getOpPtrVector().push_back(
2394 if (!internalStressTagName.empty()) {
2395 switch (internalStressInterpOrder) {
2397 fe_rhs->getOpPtrVector().push_back(
2401 fe_rhs->getOpPtrVector().push_back(
2406 "Unsupported internal stress interpolation order %d",
2407 internalStressInterpOrder);
2409 if (internalStressVoigt) {
2410 fe_rhs->getOpPtrVector().push_back(
2414 fe_rhs->getOpPtrVector().push_back(
2419 if (
auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
2420 stretchTensor, dataAtPts, externalStrainVecPtr, timeScaleMap)) {
2421 fe_rhs->getOpPtrVector().push_back(op);
2422 }
else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
2424 "OpSpatialPhysicalExternalStrain not implemented for this "
2428 fe_rhs->getOpPtrVector().push_back(
2429 physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
2432 fe_rhs->getOpPtrVector().push_back(
2434 fe_rhs->getOpPtrVector().push_back(
2436 fe_rhs->getOpPtrVector().push_back(
2439 auto set_hybridisation = [&](
auto &pip) {
2446 using SideEleOp = EleOnSide::UserDataOperator;
2447 using BdyEleOp = BoundaryEle::UserDataOperator;
2452 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2454 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2457 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2458 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2460 CHKERR EshelbianPlasticity::
2461 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2462 op_loop_skeleton_side->getOpPtrVector(), {L2},
2463 materialH1Positions, frontAdjEdges);
2467 auto broken_data_ptr =
2468 boost::make_shared<std::vector<BrokenBaseSideData>>();
2471 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2472 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2473 boost::make_shared<CGGUserPolynomialBase>();
2475 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2476 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2477 materialH1Positions, frontAdjEdges);
2478 op_loop_domain_side->getOpPtrVector().push_back(
2480 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
2481 op_loop_domain_side->getOpPtrVector().push_back(
2484 op_loop_domain_side->getOpPtrVector().push_back(
2488 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2490 GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
2492 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2493 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dHybrid(
2494 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
2495 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2496 op_loop_skeleton_side->getOpPtrVector().push_back(
2499 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dBroken(
2500 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
2503 pip.push_back(op_loop_skeleton_side);
2508 auto set_tau_stabilsation = [&](
auto &pip) {
2515 using SideEleOp = EleOnSide::UserDataOperator;
2516 using BdyEleOp = BoundaryEle::UserDataOperator;
2521 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2523 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2526 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2527 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2530 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2531 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2532 boost::make_shared<CGGUserPolynomialBase>();
2534 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2535 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2536 materialH1Positions, frontAdjEdges);
2539 auto broken_disp_data_ptr =
2540 boost::make_shared<std::vector<BrokenBaseSideData>>();
2541 op_loop_domain_side->getOpPtrVector().push_back(
2543 broken_disp_data_ptr));
2544 auto disp_mat_ptr = boost::make_shared<MatrixDouble>();
2545 if (stab_tau_dot_variant) {
2546 op_loop_domain_side->getOpPtrVector().push_back(
2550 op_loop_domain_side->getOpPtrVector().push_back(
2556 LinearForm<GAUSS>::OpBaseTimesVector<1,
SPACE_DIM, 1>;
2559 op_loop_domain_side->getOpPtrVector().push_back(
2560 new OpBaseTimesVector(spatialL2Disp, disp_mat_ptr,
2561 [local_tau_sacale](
double,
double,
double) {
2562 return *local_tau_sacale;
2566 op_loop_domain_side->getOpPtrVector().push_back(
2569 op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
2570 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2573 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2574 if (stab_tau_dot_variant) {
2575 op_loop_skeleton_side->getOpPtrVector().push_back(
2579 op_loop_skeleton_side->getOpPtrVector().push_back(
2585 op_loop_skeleton_side->getOpPtrVector().push_back(
2587 hybridSpatialDisp, hybrid_ptr,
2588 [local_tau_sacale, broken_disp_data_ptr](
double,
double,
double) {
2589 return broken_disp_data_ptr->size() * (*local_tau_sacale);
2592 op_loop_skeleton_side->getOpPtrVector().push_back(
2594 hybridSpatialDisp, broken_disp_data_ptr,
2595 [local_tau_sacale](
double,
double,
double) {
2596 return -(*local_tau_sacale);
2599 op_loop_skeleton_side->getOpPtrVector().push_back(
2601 broken_disp_data_ptr, hybrid_ptr,
2602 [local_tau_sacale](
double,
double,
double) {
2603 return -(*local_tau_sacale);
2607 pip.push_back(op_loop_skeleton_side);
2612 auto set_contact = [&](
auto &pip) {
2619 using SideEleOp = EleOnSide::UserDataOperator;
2620 using BdyEleOp = BoundaryEle::UserDataOperator;
2625 mField, contactElement,
SPACE_DIM - 1, Sev::noisy);
2627 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2628 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2629 CHKERR EshelbianPlasticity::
2630 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2631 op_loop_skeleton_side->getOpPtrVector(), {L2},
2632 materialH1Positions, frontAdjEdges);
2636 auto broken_data_ptr =
2637 boost::make_shared<std::vector<BrokenBaseSideData>>();
2640 auto contact_common_data_ptr =
2641 boost::make_shared<ContactOps::CommonData>();
2643 auto add_ops_domain_side = [&](
auto &pip) {
2647 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2648 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2649 boost::make_shared<CGGUserPolynomialBase>();
2651 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2652 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2653 materialH1Positions, frontAdjEdges);
2654 op_loop_domain_side->getOpPtrVector().push_back(
2657 op_loop_domain_side->getOpPtrVector().push_back(
2659 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2660 pip.push_back(op_loop_domain_side);
2664 auto add_ops_contact_rhs = [&](
auto &pip) {
2667 auto contact_sfd_map_range_ptr =
2668 boost::make_shared<std::map<int, Range>>(
2669 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2672 contactDisp, contact_common_data_ptr->contactDispPtr()));
2673 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2676 pip.push_back(
new OpTreeSearch(
2677 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2681 contactDisp, contact_common_data_ptr, contactTreeRhs,
2682 contact_sfd_map_range_ptr));
2685 broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
2691 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2692 CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
2695 pip.push_back(op_loop_skeleton_side);
2700 CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
2701 CHKERR set_contact(fe_rhs->getOpPtrVector());
2702 if (alphaTau > 0.0) {
2703 CHKERR set_tau_stabilsation(fe_rhs->getOpPtrVector());
2707 using BodyNaturalBC =
2709 Assembly<PETSC>::LinearForm<
GAUSS>;
2711 BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
2713 auto body_time_scale =
2714 boost::make_shared<DynamicRelaxationTimeScale>(
"body_force.txt");
2715 CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
2716 fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
2717 "BODY_FORCE", Sev::inform);
2721 fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2729 fe_lhs->getOpPtrVector().push_back(
2732 bubbleField, piolaStress, dataAtPts));
2733 fe_lhs->getOpPtrVector().push_back(
2737 fe_lhs->getOpPtrVector().push_back(
2738 physicalEquations->returnOpSpatialPhysical_du_du(
2739 stretchTensor, stretchTensor, dataAtPts, alphaU));
2741 stretchTensor, piolaStress, dataAtPts,
true));
2743 stretchTensor, bubbleField, dataAtPts,
true));
2745 stretchTensor, rotAxis, dataAtPts,
2746 symmetrySelector ==
SYMMETRIC ?
true :
false));
2750 spatialL2Disp, piolaStress, dataAtPts,
true));
2752 spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
2755 piolaStress, rotAxis, dataAtPts,
2756 symmetrySelector ==
SYMMETRIC ?
true :
false));
2758 bubbleField, rotAxis, dataAtPts,
2759 symmetrySelector ==
SYMMETRIC ?
true :
false));
2764 rotAxis, stretchTensor, dataAtPts,
false));
2767 rotAxis, piolaStress, dataAtPts,
false));
2769 rotAxis, bubbleField, dataAtPts,
false));
2772 rotAxis, rotAxis, dataAtPts, alphaOmega));
2774 auto set_hybridisation = [&](
auto &pip) {
2781 using SideEleOp = EleOnSide::UserDataOperator;
2782 using BdyEleOp = BoundaryEle::UserDataOperator;
2787 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2788 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2791 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2792 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2793 CHKERR EshelbianPlasticity::
2794 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2795 op_loop_skeleton_side->getOpPtrVector(), {L2},
2796 materialH1Positions, frontAdjEdges);
2800 auto broken_data_ptr =
2801 boost::make_shared<std::vector<BrokenBaseSideData>>();
2804 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2805 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2806 boost::make_shared<CGGUserPolynomialBase>();
2808 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2809 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2810 materialH1Positions, frontAdjEdges);
2811 op_loop_domain_side->getOpPtrVector().push_back(
2814 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2816 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2817 op_loop_skeleton_side->getOpPtrVector().push_back(
2818 new OpC(hybridSpatialDisp, broken_data_ptr,
2819 boost::make_shared<double>(1.0),
true,
false));
2822 pip.push_back(op_loop_skeleton_side);
2827 auto set_tau_stabilsation = [&](
auto &pip) {
2834 using SideEleOp = EleOnSide::UserDataOperator;
2835 using BdyEleOp = BoundaryEle::UserDataOperator;
2840 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2841 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2844 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2845 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2846 CHKERR EshelbianPlasticity::
2847 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2848 op_loop_skeleton_side->getOpPtrVector(), {L2},
2849 materialH1Positions, frontAdjEdges);
2853 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2854 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2855 boost::make_shared<CGGUserPolynomialBase>();
2857 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2858 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2859 materialH1Positions, frontAdjEdges);
2861 auto broken_disp_data_ptr =
2862 boost::make_shared<std::vector<BrokenBaseSideData>>();
2863 op_loop_domain_side->getOpPtrVector().push_back(
2865 broken_disp_data_ptr));
2868 BiLinearForm<GAUSS>::OpMass<1,
SPACE_DIM>;
2870 op_loop_domain_side->getOpPtrVector().push_back(
2871 new OpMassVectorDomain(spatialL2Disp, spatialL2Disp,
2872 [local_tau_sacale](
double,
double,
double) {
2873 return *local_tau_sacale;
2875 auto time_shift = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a; };
2876 if (stab_tau_dot_variant) {
2877 static_cast<OpMassVectorDomain &
>(
2878 op_loop_domain_side->getOpPtrVector().back())
2879 .feScalingFun = time_shift;
2881 op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
2882 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2886 hybridSpatialDisp, hybridSpatialDisp,
2887 [local_tau_sacale, broken_disp_data_ptr](
double,
double,
double) {
2888 return broken_disp_data_ptr->size() * (*local_tau_sacale);
2890 if (stab_tau_dot_variant) {
2892 op_loop_skeleton_side->getOpPtrVector().back())
2893 .feScalingFun = time_shift;
2896 op_loop_skeleton_side->getOpPtrVector().push_back(
2898 hybridSpatialDisp, broken_disp_data_ptr,
2899 [local_tau_sacale](
double,
double,
double) {
2900 return -(*local_tau_sacale);
2903 if (stab_tau_dot_variant) {
2905 op_loop_skeleton_side->getOpPtrVector().back())
2906 .feScalingFun = time_shift;
2910 op_loop_skeleton_side->getOpPtrVector().push_back(
2912 hybridSpatialDisp, broken_disp_data_ptr,
2913 [local_tau_sacale](
double,
double,
double) {
2914 return -(*local_tau_sacale);
2917 if (stab_tau_dot_variant) {
2919 op_loop_skeleton_side->getOpPtrVector().back())
2920 .feScalingFun = time_shift;
2924 pip.push_back(op_loop_skeleton_side);
2929 auto set_contact = [&](
auto &pip) {
2936 using SideEleOp = EleOnSide::UserDataOperator;
2937 using BdyEleOp = BoundaryEle::UserDataOperator;
2942 mField, contactElement,
SPACE_DIM - 1, Sev::noisy);
2944 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2945 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2946 CHKERR EshelbianPlasticity::
2947 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2948 op_loop_skeleton_side->getOpPtrVector(), {L2},
2949 materialH1Positions, frontAdjEdges);
2953 auto broken_data_ptr =
2954 boost::make_shared<std::vector<BrokenBaseSideData>>();
2957 auto contact_common_data_ptr =
2958 boost::make_shared<ContactOps::CommonData>();
2960 auto add_ops_domain_side = [&](
auto &pip) {
2964 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2965 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2966 boost::make_shared<CGGUserPolynomialBase>();
2968 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2969 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2970 materialH1Positions, frontAdjEdges);
2971 op_loop_domain_side->getOpPtrVector().push_back(
2974 op_loop_domain_side->getOpPtrVector().push_back(
2976 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2977 pip.push_back(op_loop_domain_side);
2981 auto add_ops_contact_lhs = [&](
auto &pip) {
2984 contactDisp, contact_common_data_ptr->contactDispPtr()));
2985 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2988 pip.push_back(
new OpTreeSearch(
2989 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2994 auto contact_sfd_map_range_ptr =
2995 boost::make_shared<std::map<int, Range>>(
2996 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2999 contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
3000 contact_sfd_map_range_ptr));
3003 contactDisp, broken_data_ptr, contact_common_data_ptr,
3004 contactTreeRhs, contact_sfd_map_range_ptr));
3007 broken_data_ptr, contactDisp, contact_common_data_ptr,
3014 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
3015 CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
3018 pip.push_back(op_loop_skeleton_side);
3023 CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
3024 CHKERR set_contact(fe_lhs->getOpPtrVector());
3025 if (alphaTau > 0.0) {
3026 CHKERR set_tau_stabilsation(fe_lhs->getOpPtrVector());
3037 const bool add_elastic,
const bool add_material,
3038 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
3039 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
3042 fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
3043 fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
3048 fe_rhs->getRuleHook = [](int, int, int) {
return -1; };
3049 fe_lhs->getRuleHook = [](int, int, int) {
return -1; };
3050 fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
3051 fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
3054 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3055 fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
3057 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3058 fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
3062 auto get_broken_op_side = [
this](
auto &pip) {
3065 using SideEleOp = EleOnSide::UserDataOperator;
3067 auto broken_data_ptr =
3068 boost::make_shared<std::vector<BrokenBaseSideData>>();
3071 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
3072 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3073 boost::make_shared<CGGUserPolynomialBase>();
3075 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3076 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3077 materialH1Positions, frontAdjEdges);
3078 op_loop_domain_side->getOpPtrVector().push_back(
3080 boost::shared_ptr<double> piola_scale_ptr(
new double);
3081 op_loop_domain_side->getOpPtrVector().push_back(
3082 physicalEquations->returnOpSetScale(piola_scale_ptr,
3083 physicalEquations));
3084 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
3085 op_loop_domain_side->getOpPtrVector().push_back(
3088 op_loop_domain_side->getOpPtrVector().push_back(
3090 pip.push_back(op_loop_domain_side);
3091 return std::make_tuple(broken_data_ptr, piola_scale_ptr);
3094 auto set_rhs = [&]() {
3097 auto [broken_data_ptr, piola_scale_ptr] =
3098 get_broken_op_side(fe_rhs->getOpPtrVector());
3100 fe_rhs->getOpPtrVector().push_back(
3101 new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr, timeScaleMap));
3103 broken_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
3106 broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
3108 fe_rhs->getOpPtrVector().push_back(
3110 piola_scale_ptr, timeScaleMap));
3111 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
3113 fe_rhs->getOpPtrVector().push_back(
3115 hybridSpatialDisp, hybrid_grad_ptr));
3117 hybridSpatialDisp, bcSpatialPressureVecPtr, piola_scale_ptr,
3118 hybrid_grad_ptr, timeScaleMap));
3120 hybridSpatialDisp, bcSpatialAnalyticalTractionVecPtr, piola_scale_ptr,
3123 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
3124 fe_rhs->getOpPtrVector().push_back(
3128 hybridSpatialDisp, hybrid_ptr, broken_data_ptr,
3129 bcSpatialNormalDisplacementVecPtr, timeScaleMap));
3131 auto get_normal_disp_bc_faces = [&]() {
3134 return boost::make_shared<Range>(faces);
3139 using BdyEleOp = BoundaryEle::UserDataOperator;
3141 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
3142 fe_rhs->getOpPtrVector().push_back(
new OpC_dBroken(
3143 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
3144 get_normal_disp_bc_faces()));
3149 auto set_lhs = [&]() {
3152 auto [broken_data_ptr, piola_scale_ptr] =
3153 get_broken_op_side(fe_lhs->getOpPtrVector());
3156 hybridSpatialDisp, bcSpatialNormalDisplacementVecPtr, timeScaleMap));
3158 hybridSpatialDisp, broken_data_ptr, bcSpatialNormalDisplacementVecPtr,
3161 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
3163 fe_lhs->getOpPtrVector().push_back(
3165 hybridSpatialDisp, hybrid_grad_ptr));
3167 hybridSpatialDisp, bcSpatialPressureVecPtr, hybrid_grad_ptr,
3170 auto get_normal_disp_bc_faces = [&]() {
3173 return boost::make_shared<Range>(faces);
3178 using BdyEleOp = BoundaryEle::UserDataOperator;
3180 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
3181 fe_lhs->getOpPtrVector().push_back(
new OpC(
3182 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
3183 true,
true, get_normal_disp_bc_faces()));
3198 boost::shared_ptr<ContactTree> &fe_contact_tree
3204 auto get_body_range = [
this](
auto name,
int dim,
auto sev) {
3205 std::map<int, Range> map;
3208 mField.getInterface<
MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
3210 (boost::format(
"%s(.*)") % name).str()
3216 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
3219 map[m_ptr->getMeshsetId()] = ents;
3220 MOFEM_LOG(
"EPSYNC", sev) <<
"Meshset: " << m_ptr->getMeshsetId() <<
" "
3221 << ents.size() <<
" entities";
3228 auto get_map_skin = [
this](
auto &&map) {
3229 ParallelComm *pcomm =
3232 Skinner skin(&mField.get_moab());
3233 for (
auto &
m : map) {
3235 CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
3237 PSTATUS_SHARED | PSTATUS_MULTISHARED,
3238 PSTATUS_NOT, -1,
nullptr),
3240 m.second.swap(skin_faces);
3250 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3252 auto calcs_side_traction = [&](
auto &pip) {
3256 using SideEleOp = EleOnSide::UserDataOperator;
3258 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
3259 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3260 boost::make_shared<CGGUserPolynomialBase>();
3261 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3262 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3263 materialH1Positions, frontAdjEdges);
3264 op_loop_domain_side->getOpPtrVector().push_back(
3266 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3267 boost::make_shared<double>(1.0)));
3268 pip.push_back(op_loop_domain_side);
3272 auto add_contact_three = [&]() {
3274 auto tree_moab_ptr = boost::make_shared<moab::Core>();
3275 fe_contact_tree = boost::make_shared<ContactTree>(
3276 mField, tree_moab_ptr, spaceOrder,
3277 get_body_range(
"CONTACT",
SPACE_DIM - 1, Sev::inform));
3278 fe_contact_tree->getOpPtrVector().push_back(
3280 contactDisp, contact_common_data_ptr->contactDispPtr()));
3281 CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
3282 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3283 fe_contact_tree->getOpPtrVector().push_back(
3285 fe_contact_tree->getOpPtrVector().push_back(
3286 new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
3290 CHKERR add_contact_three();
3300 CHKERR setContactElementRhsOps(contactTreeRhs);
3302 CHKERR setVolumeElementOps(tag,
true,
false, elasticFeRhs, elasticFeLhs);
3303 CHKERR setFaceElementOps(
true,
false, elasticBcRhs, elasticBcLhs);
3306 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
3308 auto get_op_contact_bc = [&]() {
3311 mField, contactElement,
SPACE_DIM - 1, Sev::noisy, adj_cache);
3312 return op_loop_side;
3320 boost::shared_ptr<FEMethod> null;
3322 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3355 bool set_ts_monitor) {
3357#ifdef ENABLE_PYTHON_BINDING
3358 auto setup_sdf = [&]() {
3359 boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
3361 auto file_exists = [](std::string myfile) {
3362 std::ifstream file(myfile.c_str());
3369 char sdf_file_name[255] =
"sdf.py";
3371 sdf_file_name, 255, PETSC_NULLPTR);
3373 if (file_exists(sdf_file_name)) {
3374 MOFEM_LOG(
"EP", Sev::inform) << sdf_file_name <<
" file found";
3375 sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
3376 CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
3377 ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
3378 MOFEM_LOG(
"EP", Sev::inform) <<
"SdfPython initialized";
3380 MOFEM_LOG(
"EP", Sev::warning) << sdf_file_name <<
" file NOT found";
3382 return sdf_python_ptr;
3384 auto sdf_python_ptr = setup_sdf();
3387 auto setup_ts_monitor = [&]() {
3388 boost::shared_ptr<TsCtx>
ts_ctx;
3390 if (set_ts_monitor) {
3392 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
3396 MOFEM_LOG(
"EP", Sev::inform) <<
"TS monitor setup";
3397 return std::make_tuple(
ts_ctx);
3400 auto setup_snes_monitor = [&]() {
3403 CHKERR TSGetSNES(ts, &snes);
3405 CHKERR SNESMonitorSet(snes,
3408 (
void *)(snes_ctx.get()), PETSC_NULLPTR);
3409 MOFEM_LOG(
"EP", Sev::inform) <<
"SNES monitor setup";
3413 auto setup_snes_conergence_test = [&]() {
3416 auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
3417 PetscReal snorm, PetscReal fnorm,
3418 SNESConvergedReason *reason,
void *cctx) {
3421 CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
3425 CHKERR SNESGetSolutionUpdate(snes, &x_update);
3426 CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
3494 auto setup_section = [&]() {
3495 PetscSection section_raw;
3498 CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
3499 for (
int ff = 0; ff != num_fields; ff++) {
3507 auto set_vector_on_mesh = [&]() {
3511 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
3512 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
3513 MOFEM_LOG(
"EP", Sev::inform) <<
"Vector set on mesh";
3517 auto setup_schur_block_solver = [&]() {
3518 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver";
3519 CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
3520 CHKERR TSSetFromOptions(ts);
3523 boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
3524 if constexpr (
A == AssemblyType::BLOCK_MAT) {
3529 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver done";
3536#ifdef ENABLE_PYTHON_BINDING
3537 return std::make_tuple(setup_sdf(), setup_ts_monitor(),
3538 setup_snes_monitor(), setup_snes_conergence_test(),
3539 setup_section(), set_vector_on_mesh(),
3540 setup_schur_block_solver());
3542 return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
3543 setup_snes_conergence_test(), setup_section(),
3544 set_vector_on_mesh(), setup_schur_block_solver());
3552 PetscBool debug_model = PETSC_FALSE;
3556 if (debug_model == PETSC_TRUE) {
3558 auto post_proc = [&](TS ts, PetscReal
t, Vec u, Vec u_t, Vec u_tt, Vec
F,
3563 CHKERR TSGetSNES(ts, &snes);
3565 CHKERR SNESGetIterationNumber(snes, &it);
3566 std::string file_name =
"snes_iteration_" + std::to_string(it) +
".h5m";
3567 CHKERR postProcessResults(1, file_name,
F, u_t);
3568 std::string file_skel_name =
3569 "snes_iteration_skel_" + std::to_string(it) +
".h5m";
3571 auto get_material_force_tag = [&]() {
3572 auto &moab = mField.get_moab();
3579 CHKERR calculateFaceMaterialForce(1, ts);
3580 CHKERR postProcessSkeletonResults(1, file_skel_name,
F,
3581 {get_material_force_tag()});
3585 ts_ctx_ptr->tsDebugHook = post_proc;
3590 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3592 CHKERR VecDuplicate(x, &xx);
3593 CHKERR VecZeroEntries(xx);
3594 CHKERR TS2SetSolution(ts, x, xx);
3597 CHKERR TSSetSolution(ts, x);
3600 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3601 {elasticFeLhs.get(), elasticFeRhs.get()});
3606 CHKERR TSSolve(ts, PETSC_NULLPTR);
3608 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3609 {elasticFeLhs.get(), elasticFeRhs.get()});
3612 CHKERR TSGetSNES(ts, &snes);
3613 int lin_solver_iterations;
3614 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
3616 <<
"Number of linear solver iterations " << lin_solver_iterations;
3618 PetscBool test_cook_flg = PETSC_FALSE;
3621 if (test_cook_flg) {
3622 constexpr int expected_lin_solver_iterations = 11;
3623 if (lin_solver_iterations > expected_lin_solver_iterations)
3626 "Expected number of iterations is different than expected %d > %d",
3627 lin_solver_iterations, expected_lin_solver_iterations);
3630 PetscBool test_sslv116_flag = PETSC_FALSE;
3632 &test_sslv116_flag, PETSC_NULLPTR);
3634 if (test_sslv116_flag) {
3635 double max_val = 0.0;
3636 double min_val = 0.0;
3637 auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3639 auto ent_type = ent_ptr->getEntType();
3640 if (ent_type == MBVERTEX) {
3641 max_val = std::max(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], max_val);
3642 min_val = std::min(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], min_val);
3647 field_min_max, spatialH1Disp);
3649 double global_max_val = 0.0;
3650 double global_min_val = 0.0;
3651 MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
3653 MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
3656 <<
"Max " << spatialH1Disp <<
" value: " << global_max_val;
3658 <<
"Min " << spatialH1Disp <<
" value: " << global_min_val;
3660 double ref_max_val = 0.00767;
3661 double ref_min_val = -0.00329;
3662 if (std::abs(global_max_val - ref_max_val) > 1e-5) {
3664 "Incorrect max value of the displacement field: %f != %f",
3665 global_max_val, ref_max_val);
3667 if (std::abs(global_min_val - ref_min_val) > 4e-5) {
3669 "Incorrect min value of the displacement field: %f != %f",
3670 global_min_val, ref_min_val);
3681 double start_time) {
3686 double final_time = 1;
3687 double delta_time = 0.1;
3689 PetscBool ts_h1_update = PETSC_FALSE;
3691 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Dynamic Relaxation Options",
3694 CHKERR PetscOptionsScalar(
"-dynamic_final_time",
3695 "dynamic relaxation final time",
"", final_time,
3696 &final_time, PETSC_NULLPTR);
3697 CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
3698 "dynamic relaxation final time",
"", delta_time,
3699 &delta_time, PETSC_NULLPTR);
3700 CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
3701 max_it, &max_it, PETSC_NULLPTR);
3702 CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
3703 ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
3707 auto setup_ts_monitor = [&]() {
3708 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*
this);
3711 auto monitor_ptr = setup_ts_monitor();
3713 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3714 {elasticFeLhs.get(), elasticFeRhs.get()});
3718 double ts_delta_time;
3719 CHKERR TSGetTimeStep(ts, &ts_delta_time);
3729 dynamicTime = start_time;
3730 dynamicStep = start_step;
3731 monitor_ptr->ts_u = PETSC_NULLPTR;
3732 monitor_ptr->ts_t = dynamicTime;
3733 monitor_ptr->ts_step = dynamicStep;
3735 MOFEM_LOG(
"EP", Sev::inform) <<
"IT " << dynamicStep <<
" Time "
3736 << dynamicTime <<
" delta time " << delta_time;
3737 dynamicTime += delta_time;
3740 for (; dynamicTime < final_time; dynamicTime += delta_time) {
3741 MOFEM_LOG(
"EP", Sev::inform) <<
"IT " << dynamicStep <<
" Time "
3742 << dynamicTime <<
" delta time " << delta_time;
3744 CHKERR TSSetStepNumber(ts, 0);
3746 CHKERR TSSetTimeStep(ts, ts_delta_time);
3747 if (!ts_h1_update) {
3750 CHKERR TSSolve(ts, PETSC_NULLPTR);
3751 if (!ts_h1_update) {
3757 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
3758 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
3760 monitor_ptr->ts_u = x;
3761 monitor_ptr->ts_t = dynamicTime;
3762 monitor_ptr->ts_step = dynamicStep;
3766 if (dynamicStep > max_it)
3771 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3772 {elasticFeLhs.get(), elasticFeRhs.get()});
3780 auto set_block = [&](
auto name,
int dim) {
3781 std::map<int, Range> map;
3782 auto set_tag_impl = [&](
auto name) {
3787 std::regex((boost::format(
"%s(.*)") % name).str())
3790 for (
auto bc : bcs) {
3792 CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
3794 map[bc->getMeshsetId()] = r;
3799 CHKERR set_tag_impl(name);
3801 return std::make_pair(name, map);
3804 auto set_skin = [&](
auto &&map) {
3805 for (
auto &
m : map.second) {
3812 auto set_tag = [&](
auto &&map) {
3814 auto name = map.first;
3815 int def_val[] = {-1};
3817 mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER,
th,
3818 MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
3820 for (
auto &
m : map.second) {
3828 listTagsToTransfer.push_back(set_tag(set_skin(set_block(
"BODY", 3))));
3829 listTagsToTransfer.push_back(set_tag(set_skin(set_block(
"MAT_ELASTIC", 3))));
3830 listTagsToTransfer.push_back(
3831 set_tag(set_skin(set_block(
"MAT_NEOHOOKEAN", 3))));
3832 listTagsToTransfer.push_back(set_tag(set_block(
"CONTACT", 2)));
3839 Vec f_residual, Vec var_vector,
3840 std::vector<Tag> tags_to_transfer) {
3846 rval = mField.get_moab().tag_get_handle(
"CRACK",
th);
3847 if (
rval == MB_SUCCESS) {
3848 CHKERR mField.get_moab().tag_delete(
th);
3850 int def_val[] = {0};
3851 CHKERR mField.get_moab().tag_get_handle(
3852 "CRACK", 1, MB_TYPE_INTEGER,
th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
3854 CHKERR mField.get_moab().tag_clear_data(
th, *crackFaces, mark);
3855 tags_to_transfer.push_back(
th);
3857 auto get_tag = [&](
auto name,
auto dim) {
3858 auto &mob = mField.get_moab();
3860 double def_val[] = {0., 0., 0.};
3861 CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3862 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3867 tags_to_transfer.push_back(get_tag(
"MaterialForce", 3));
3872 for (
auto t : listTagsToTransfer) {
3873 tags_to_transfer.push_back(
t);
3881 struct exclude_sdf {
3882 exclude_sdf(
Range &&r) : map(r) {}
3883 bool operator()(
FEMethod *fe_method_ptr) {
3885 if (map.find(ent) != map.end()) {
3895 contactTreeRhs->exeTestHook =
3900 auto get_post_proc = [&](
auto &post_proc_mesh,
auto sense) {
3902 auto post_proc_ptr =
3903 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3904 mField, post_proc_mesh);
3905 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3906 post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
3909 auto domain_ops = [&](
auto &fe,
int sense) {
3911 fe.getUserPolynomialBase() =
3913 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3914 fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
3916 auto piola_scale_ptr = boost::make_shared<double>(1.0);
3917 fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
3918 piola_scale_ptr, physicalEquations));
3920 piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
3922 bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
3925 fe.getOpPtrVector().push_back(
3926 physicalEquations->returnOpCalculateStretchFromStress(
3927 dataAtPts, physicalEquations));
3929 fe.getOpPtrVector().push_back(
3931 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
3936 piolaStress, dataAtPts->getVarPiolaPts(),
3937 boost::make_shared<double>(1), vec));
3939 bubbleField, dataAtPts->getVarPiolaPts(),
3940 boost::make_shared<double>(1), vec, MBMAXTYPE));
3942 rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
3944 fe.getOpPtrVector().push_back(
3945 physicalEquations->returnOpCalculateVarStretchFromStress(
3946 dataAtPts, physicalEquations));
3948 fe.getOpPtrVector().push_back(
3950 stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
3955 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
3957 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
3960 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
3962 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
3964 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
3966 fe.getOpPtrVector().push_back(
3970 fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
3971 tag,
true,
false, dataAtPts, physicalEquations));
3973 physicalEquations->returnOpCalculateEnergy(dataAtPts,
nullptr)) {
3974 fe.getOpPtrVector().push_back(op);
3983 struct OpSidePPMap :
public OpPPMap {
3984 OpSidePPMap(moab::Interface &post_proc_mesh,
3985 std::vector<EntityHandle> &map_gauss_pts,
3986 DataMapVec data_map_scalar, DataMapMat data_map_vec,
3987 DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
3989 :
OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
3990 data_map_vec, data_map_mat, data_symm_map_mat),
3997 if (tagSense != OpPPMap::getSkeletonSense())
4009 vec_fields[
"SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
4010 vec_fields[
"SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
4011 vec_fields[
"Omega"] = dataAtPts->getRotAxisAtPts();
4012 vec_fields[
"ContactDisplacement"] = dataAtPts->getContactL2AtPts();
4013 vec_fields[
"AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
4014 vec_fields[
"X"] = dataAtPts->getLargeXH1AtPts();
4016 vec_fields[
"EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
4020 vec_fields[
"VarOmega"] = dataAtPts->getVarRotAxisPts();
4021 vec_fields[
"VarSpatialDisplacementL2"] =
4022 boost::make_shared<MatrixDouble>();
4024 spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], vec, MBTET));
4028 vec_fields[
"ResSpatialDisplacementL2"] =
4029 boost::make_shared<MatrixDouble>();
4031 spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], vec, MBTET));
4032 vec_fields[
"ResOmega"] = boost::make_shared<MatrixDouble>();
4034 rotAxis, vec_fields[
"ResOmega"], vec, MBTET));
4038 mat_fields[
"PiolaStress"] = dataAtPts->getApproxPAtPts();
4040 mat_fields[
"VarPiolaStress"] = dataAtPts->getVarPiolaPts();
4044 mat_fields[
"ResPiolaStress"] = boost::make_shared<MatrixDouble>();
4046 piolaStress, mat_fields[
"ResPiolaStress"],
4047 boost::make_shared<double>(1), vec));
4049 bubbleField, mat_fields[
"ResPiolaStress"],
4050 boost::make_shared<double>(1), vec, MBMAXTYPE));
4052 if (!internalStressTagName.empty()) {
4053 mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
4054 switch (internalStressInterpOrder) {
4056 fe.getOpPtrVector().push_back(
4060 fe.getOpPtrVector().push_back(
4065 "Unsupported internal stress interpolation order %d",
4066 internalStressInterpOrder);
4071 mat_fields_symm[
"LogSpatialStretch"] =
4072 dataAtPts->getLogStretchTensorAtPts();
4073 mat_fields_symm[
"SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
4075 mat_fields_symm[
"VarLogSpatialStretch"] =
4076 dataAtPts->getVarLogStreachPts();
4081 mat_fields_symm[
"ResLogSpatialStretch"] =
4082 boost::make_shared<MatrixDouble>();
4083 fe.getOpPtrVector().push_back(
4085 stretchTensor, mat_fields_symm[
"ResLogSpatialStretch"], vec,
4090 fe.getOpPtrVector().push_back(
4094 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4111 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4117 post_proc_ptr->getOpPtrVector().push_back(
4119 dataAtPts->getContactL2AtPts()));
4121 auto X_h1_ptr = boost::make_shared<MatrixDouble>();
4123 post_proc_ptr->getOpPtrVector().push_back(
4125 dataAtPts->getLargeXH1AtPts()));
4130 domain_ops(*(op_loop_side->getSideFEPtr()), sense);
4131 post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
4133 return post_proc_ptr;
4137 auto calcs_side_traction_and_displacements = [&](
auto &post_proc_ptr,
4140 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
4144 using SideEleOp = EleOnSide::UserDataOperator;
4146 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
4147 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4149 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4150 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4151 materialH1Positions, frontAdjEdges);
4152 op_loop_domain_side->getOpPtrVector().push_back(
4154 piolaStress, contact_common_data_ptr->contactTractionPtr(),
4155 boost::make_shared<double>(1.0)));
4156 pip.push_back(op_loop_domain_side);
4158 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
4161 contactDisp, contact_common_data_ptr->contactDispPtr()));
4162 pip.push_back(
new OpTreeSearch(
4163 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
4165 &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
4172 pip.push_back(op_this);
4173 auto contact_residual = boost::make_shared<MatrixDouble>();
4174 op_this->getOpPtrVector().push_back(
4176 contactDisp, contact_residual,
4179 op_this->getOpPtrVector().push_back(
4183 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4187 {{
"res_contact", contact_residual}},
4201 auto post_proc_mesh = boost::make_shared<moab::Core>();
4202 auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
4203 auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
4204 CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
4205 post_proc_ptr->getOpPtrVector());
4211 CHKERR mField.get_moab().get_adjacencies(own_tets,
SPACE_DIM - 1,
true,
4212 own_faces, moab::Interface::UNION);
4214 auto get_post_negative = [&](
auto &&ents) {
4215 auto crack_faces_pos = ents;
4216 auto crack_faces_neg = crack_faces_pos;
4217 auto skin =
get_skin(mField, own_tets);
4218 auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
4219 for (
auto f : crack_on_proc_skin) {
4221 CHKERR mField.get_moab().get_adjacencies(&f, 1,
SPACE_DIM,
false, tet);
4222 tet = intersect(tet, own_tets);
4223 int side_number, sense, offset;
4224 CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
4227 crack_faces_neg.erase(f);
4229 crack_faces_pos.erase(f);
4232 return std::make_pair(crack_faces_pos, crack_faces_neg);
4235 auto get_crack_faces = [&](
auto crack_faces) {
4236 auto get_adj = [&](
auto e,
auto dim) {
4238 CHKERR mField.get_moab().get_adjacencies(e, dim,
true, adj,
4239 moab::Interface::UNION);
4242 auto tets = get_adj(crack_faces, 3);
4243 auto faces = subtract(get_adj(tets, 2), crack_faces);
4244 tets = subtract(tets, get_adj(faces, 3));
4245 return subtract(crack_faces, get_adj(tets, 2));
4248 auto [crack_faces_pos, crack_faces_neg] =
4249 get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
4251 auto only_crack_faces = [&](
FEMethod *fe_method_ptr) {
4252 auto ent = fe_method_ptr->getFEEntityHandle();
4253 if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
4259 auto only_negative_crack_faces = [&](
FEMethod *fe_method_ptr) {
4260 auto ent = fe_method_ptr->getFEEntityHandle();
4261 if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
4267 auto get_adj_front = [&]() {
4268 auto skeleton_faces = *skeletonFaces;
4270 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2,
true, adj_front,
4271 moab::Interface::UNION);
4273 adj_front = intersect(adj_front, skeleton_faces);
4274 adj_front = subtract(adj_front, *crackFaces);
4275 adj_front = intersect(own_faces, adj_front);
4279 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4280 post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
4282 auto post_proc_begin =
4286 post_proc_ptr->exeTestHook = only_crack_faces;
4287 post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
4289 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4291 post_proc_negative_sense_ptr, 0,
4292 mField.get_comm_size());
4294 constexpr bool debug =
false;
4297 auto [adj_front_pos, adj_front_neg] =
4300 auto only_front_faces_pos = [adj_front_pos](
FEMethod *fe_method_ptr) {
4301 auto ent = fe_method_ptr->getFEEntityHandle();
4302 if (adj_front_pos.find(ent) == adj_front_pos.end()) {
4308 auto only_front_faces_neg = [adj_front_neg](
FEMethod *fe_method_ptr) {
4309 auto ent = fe_method_ptr->getFEEntityHandle();
4310 if (adj_front_neg.find(ent) == adj_front_neg.end()) {
4316 post_proc_ptr->exeTestHook = only_front_faces_pos;
4318 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4319 post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
4321 post_proc_negative_sense_ptr, 0,
4322 mField.get_comm_size());
4327 CHKERR post_proc_end.writeFile(file.c_str());
4334 std::vector<Tag> tags_to_transfer) {
4339 auto post_proc_mesh = boost::make_shared<moab::Core>();
4340 auto post_proc_ptr =
4341 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
4343 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
4347 auto hybrid_disp = boost::make_shared<MatrixDouble>();
4348 post_proc_ptr->getOpPtrVector().push_back(
4350 post_proc_ptr->getOpPtrVector().push_back(
4354 auto op_loop_domain_side =
4357 post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4360 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4361 boost::make_shared<CGGUserPolynomialBase>();
4362 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4363 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4365 op_loop_domain_side->getOpPtrVector().push_back(
4368 op_loop_domain_side->getOpPtrVector().push_back(
4371 op_loop_domain_side->getOpPtrVector().push_back(
4374 op_loop_domain_side->getOpPtrVector().push_back(
4379 op_loop_domain_side->getOpPtrVector().push_back(
4383 op_loop_domain_side->getOpPtrVector().push_back(
4391 vec_fields[
"HybridDisplacement"] = hybrid_disp;
4393 vec_fields[
"spatialL2Disp"] =
dataAtPts->getSmallWL2AtPts();
4394 vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
4396 mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
4397 mat_fields[
"HybridDisplacementGradient"] =
4400 mat_fields_symm[
"LogSpatialStretch"] =
dataAtPts->getLogStretchTensorAtPts();
4402 post_proc_ptr->getOpPtrVector().push_back(
4406 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4421 auto hybrid_res = boost::make_shared<MatrixDouble>();
4422 post_proc_ptr->getOpPtrVector().push_back(
4427 post_proc_ptr->getOpPtrVector().push_back(
4431 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4435 {{
"res_hybrid", hybrid_res}},
4446 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4448 auto post_proc_begin =
4455 CHKERR post_proc_end.writeFile(file.c_str());
4463 constexpr bool debug =
false;
4465 auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
4466 std::vector<Tag> tags;
4467 tags.reserve(names.size());
4468 auto create_and_clean = [&]() {
4470 for (
auto n : names) {
4471 tags.push_back(
Tag());
4472 auto &tag = tags.back();
4474 auto rval = moab.tag_get_handle(
n.first.c_str(), tag);
4475 if (
rval == MB_SUCCESS) {
4476 moab.tag_delete(tag);
4478 double def_val[] = {0., 0., 0.};
4479 CHKERR moab.tag_get_handle(
n.first.c_str(),
n.second, MB_TYPE_DOUBLE,
4480 tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4488 enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
4490 auto tags = get_tags_vec({{
"MaterialForce", 3},
4492 {
"GriffithForce", 1},
4493 {
"FacePressure", 1}});
4495 auto calculate_material_forces = [&]() {
4501 auto get_face_material_force_fe = [&]() {
4503 auto fe_ptr = boost::make_shared<FaceEle>(
mField);
4504 fe_ptr->getRuleHook = [](int, int, int) {
return -1; };
4505 fe_ptr->setRuleHook =
4509 EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
4513 fe_ptr->getOpPtrVector().push_back(
4516 auto op_loop_domain_side =
4519 fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4523 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4524 boost::make_shared<CGGUserPolynomialBase>();
4526 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4527 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4529 op_loop_domain_side->getOpPtrVector().push_back(
4532 op_loop_domain_side->getOpPtrVector().push_back(
4536 op_loop_domain_side->getOpPtrVector().push_back(
4544 op_loop_domain_side->getOpPtrVector().push_back(
4553 op_loop_domain_side->getOpPtrVector().push_back(
4558 op_loop_domain_side->getOpPtrVector().push_back(
4564 auto integrate_face_material_force_fe = [&](
auto &&face_energy_fe) {
4572 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
4575 CHKERR VecGetLocalSize(
v.second, &size);
4577 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
4578 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( "
4579 << low <<
" " << high <<
" ) ";
4583 CHKERR print_loc_size(face_exchange,
"material face_exchange",
4587 mField.
get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
4594 "front_skin_faces_material_force_" +
4603 CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
4608 auto calculate_front_material_force = [&](
auto nb_J_integral_contours) {
4612 auto get_conn = [&](
auto e) {
4615 "get connectivity");
4619 auto get_conn_range = [&](
auto e) {
4622 "get connectivity");
4626 auto get_adj = [&](
auto e,
auto dim) {
4633 auto get_adj_range = [&](
auto e,
auto dim) {
4636 moab::Interface::UNION),
4641 auto get_material_force = [&](
auto r,
auto th) {
4646 return material_forces;
4651 auto crack_edges = get_adj_range(*
crackFaces, 1);
4652 auto front_nodes = get_conn_range(*
frontEdges);
4659 Range all_skin_faces;
4660 Range all_front_faces;
4663 auto calculate_edge_direction = [&](
auto e) {
4668 "get connectivity");
4669 std::array<double, 6> coords;
4674 &coords[0], &coords[1], &coords[2]};
4676 &coords[3], &coords[4], &coords[5]};
4679 t_dir(
i) = t_p1(
i) - t_p0(
i);
4684 auto calculate_force_through_node = [&]() {
4696 auto body_skin_conn = get_conn_range(body_skin);
4699 for (
auto n : front_nodes) {
4700 auto adj_tets = get_adj(
n, 3);
4701 for (
int ll = 0; ll < nb_J_integral_contours; ++ll) {
4702 auto conn = get_conn_range(adj_tets);
4703 adj_tets = get_adj_range(conn, 3);
4706 auto material_forces = get_material_force(skin_faces, tags[0]);
4710 all_skin_faces.merge(skin_faces);
4714 auto calculate_node_material_force = [&]() {
4716 getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
4718 for (
auto face : skin_faces) {
4721 t_face_force_tmp(
I) = t_face_T(
I);
4724 auto face_tets = intersect(get_adj(face, 3), adj_tets);
4726 if (face_tets.empty()) {
4730 if (face_tets.size() != 1) {
4732 "face_tets.size() != 1");
4735 int side_number, sense, offset;
4739 "moab side number");
4740 t_face_force_tmp(
I) *= sense;
4741 t_node_force(
I) += t_face_force_tmp(
I);
4744 return t_node_force;
4747 auto calculate_crack_area_growth_direction =
4748 [&](
auto n,
auto &t_node_force) {
4751 auto boundary_node = intersect(
Range(
n,
n), body_skin_conn);
4752 if (boundary_node.size()) {
4753 auto faces = intersect(get_adj(
n, 2), body_skin);
4754 for (
auto f : faces) {
4757 f, &t_normal_face(0));
4758 t_project(
I) += t_normal_face(
I);
4760 t_project.normalize();
4767 if (boundary_node.size()) {
4768 t_Q(
I,
J) -= t_project(
I) * t_project(
J);
4771 auto adj_faces = intersect(get_adj(
n, 2), *
crackFaces);
4772 if (adj_faces.empty()) {
4774 intersect(get_adj(
n, 1), unite(*
frontEdges, body_edges));
4776 for (
auto e : adj_edges) {
4777 auto t_dir = calculate_edge_direction(e);
4783 t_node_force_tmp(
I) = t_node_force(
I);
4785 t_area_dir(
I) = -t_node_force_tmp(
I);
4786 t_area_dir(
I) *=
l / 2;
4791 auto front_edges = get_adj(
n, 1);
4793 for (
auto f : adj_faces) {
4798 std::array<double, 9> coords;
4804 coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
4805 auto n_it = std::find(conn, conn + num_nodes,
n);
4806 auto n_index = std::distance(conn, n_it);
4809 t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
4810 t_d_normal(0, n_index, 2),
4812 t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
4813 t_d_normal(1, n_index, 2),
4815 t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
4816 t_d_normal(2, n_index, 2)};
4819 t_projected_hessian(
I,
J) =
4820 t_Q(
I,
K) * (t_face_hessian(
K,
L) * t_Q(
L,
J));
4823 t_face_normal(
I) * t_projected_hessian(
I,
K) / 2.;
4829 auto t_node_force = calculate_node_material_force();
4833 &
n, 1, &t_node_force(0)),
4836 auto get_area_dir = [&]() {
4838 auto adj_edges = intersect(get_adj_range(adj_tets, 1),
4840 auto seed_n = get_conn_range(adj_edges);
4842 skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
4843 seed_n = subtract(seed_n, skin_adj_edges);
4845 for (
auto sn : seed_n) {
4846 auto t_area_dir_sn =
4847 calculate_crack_area_growth_direction(sn, t_node_force);
4848 t_area_dir(
I) += t_area_dir_sn(
I);
4850 for (
auto sn : skin_adj_edges) {
4851 auto t_area_dir_sn =
4852 calculate_crack_area_growth_direction(sn, t_node_force);
4853 t_area_dir(
I) += t_area_dir_sn(
I) / 2;
4858 auto t_area_dir = get_area_dir();
4864 auto griffith = -t_node_force(
I) * t_area_dir(
I) /
4865 (t_area_dir(
K) * t_area_dir(
K));
4873 auto ave_node_force = [&](
auto th) {
4878 auto conn = get_conn(e);
4879 auto data = get_material_force(conn,
th);
4880 auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
4882 for (
auto n : conn) {
4884 t_edge(
I) += t_node(
I);
4887 t_edge(
I) /= conn.size();
4890 calculate_edge_direction(e);
4908 auto ave_node_griffith_energy = [&](
auto th) {
4913 tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
4916 &e, 1, &t_edge_area_dir(0));
4917 double griffith_energy = -t_edge_force(
I) * t_edge_area_dir(
I) /
4918 (t_edge_area_dir(
K) * t_edge_area_dir(
K));
4924 CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
4925 CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
4926 CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
4931 CHKERR calculate_force_through_node();
4935 auto adj_faces = get_adj(e, 2);
4936 auto crack_face = intersect(get_adj(e, 2), *
crackFaces);
4940 all_front_faces.merge(adj_faces);
4946 &e, 1, &t_edge_force(0));
4948 calculate_edge_direction(e);
4957 for (
auto f : adj_faces) {
4961 int side_number, sense, offset;
4964 auto dot = -sense * t_cross(
I) * t_normal(
I);
4966 tags[ExhangeTags::GRIFFITHFORCE], &
f, 1, &dot),
4974 CHKERR TSGetStepNumber(ts, &ts_step);
4976 "front_edges_material_force_" +
4977 std::to_string(ts_step) +
".vtk",
4980 "front_skin_faces_material_force_" +
4981 std::to_string(ts_step) +
".vtk",
4984 "front_faces_material_force_" +
4985 std::to_string(ts_step) +
".vtk",
4994 mField.
get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
4996 mField.
get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
5003 auto print_results = [&](
auto nb_J_integral_conturs) {
5006 auto get_conn_range = [&](
auto e) {
5009 "get connectivity");
5013 auto get_tag_data = [&](
auto &ents,
auto tag,
auto dim) {
5014 std::vector<double> data(ents.size() * dim);
5021 auto at_nodes = [&]() {
5024 auto material_force =
5025 get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
5026 auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
5027 auto griffith_force =
5028 get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
5029 std::vector<double> coords(conn.size() * 3);
5033 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Material force at nodes";
5034 for (
size_t i = 0;
i < conn.size(); ++
i) {
5036 <<
"Node " << conn[
i] <<
" coords " << coords[
i * 3 + 0] <<
" "
5037 << coords[
i * 3 + 1] <<
" " << coords[
i * 3 + 2]
5038 <<
" material force " << material_force[
i * 3 + 0] <<
" "
5039 << material_force[
i * 3 + 1] <<
" " << material_force[
i * 3 + 2]
5040 <<
" area growth " << area_growth[
i * 3 + 0] <<
" "
5041 << area_growth[
i * 3 + 1] <<
" " << area_growth[
i * 3 + 2]
5042 <<
" griffith force " << std::setprecision(12)
5043 << griffith_force[
i] <<
" contour " << nb_J_integral_conturs;
5053 CHKERR calculate_material_forces();
5055 PetscBool all_contours = PETSC_FALSE;
5057 "-calculate_J_integral_all_levels", &all_contours,
5060 "-calculate_J_integral_all_contours", &all_contours,
5064 if (all_contours == PETSC_TRUE) {
5066 CHKERR calculate_front_material_force(
l);
5080 bool set_orientation) {
5083 constexpr bool debug =
false;
5084 constexpr auto sev = Sev::verbose;
5089 Range body_skin_edges;
5091 moab::Interface::UNION);
5092 Range boundary_skin_verts;
5094 boundary_skin_verts,
true);
5097 Range geometry_edges_verts;
5099 geometry_edges_verts,
true);
5100 Range crack_faces_verts;
5103 Range crack_faces_edges;
5105 *
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
5106 Range crack_faces_tets;
5108 *
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
5114 moab::Interface::UNION);
5115 Range front_verts_edges;
5117 front_verts, 1,
true, front_verts_edges, moab::Interface::UNION);
5119 auto get_tags_vec = [&](
auto tag_name,
int dim) {
5120 std::vector<Tag> tags(1);
5125 auto create_and_clean = [&]() {
5128 auto rval = moab.tag_get_handle(tag_name, tags[0]);
5129 if (
rval == MB_SUCCESS) {
5130 moab.tag_delete(tags[0]);
5132 double def_val[] = {0., 0., 0.};
5133 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
5134 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
5143 auto get_adj_front = [&](
bool subtract_crack) {
5146 adj_front, moab::Interface::UNION);
5148 adj_front = subtract(adj_front, *
crackFaces);
5154 auto th_front_position = get_tags_vec(
"FrontPosition", 3);
5155 auto th_max_face_energy = get_tags_vec(
"MaxFaceEnergy", 1);
5159 auto get_crack_adj_tets = [&](
auto r) {
5160 Range crack_faces_conn;
5162 Range crack_faces_conn_tets;
5164 true, crack_faces_conn_tets,
5165 moab::Interface::UNION);
5166 return crack_faces_conn_tets;
5169 auto get_layers_for_sides = [&](
auto &side) {
5170 std::vector<Range> layers;
5174 auto get_adj = [&](
auto &r,
int dim) {
5177 moab::Interface::UNION);
5181 auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
5186 Range front_faces = get_adj(front_nodes, 2);
5187 front_faces = subtract(front_faces, *
crackFaces);
5188 auto front_tets = get_tets(front_nodes);
5189 auto front_side = intersect(side, front_tets);
5190 layers.push_back(front_side);
5193 adj_faces = intersect(adj_faces, front_faces);
5194 auto adj_faces_tets = get_tets(adj_faces);
5195 adj_faces_tets = intersect(adj_faces_tets, front_tets);
5196 layers.push_back(unite(layers.back(), adj_faces_tets));
5197 if (layers.back().size() == layers[layers.size() - 2].size()) {
5208 auto layers_top = get_layers_for_sides(sides_pair.first);
5209 auto layers_bottom = get_layers_for_sides(sides_pair.second);
5221 MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
5223 for (
auto &r : layers_top) {
5224 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
5227 "layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
5232 for (
auto &r : layers_bottom) {
5233 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
5236 "layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
5242 auto get_cross = [&](
auto t_dir,
auto f) {
5254 auto get_sense = [&](
auto f,
auto e) {
5255 int side, sense, offset;
5258 return std::make_tuple(side, sense, offset);
5261 auto calculate_edge_direction = [&](
auto e,
auto normalize =
true) {
5265 std::array<double, 6> coords;
5268 &coords[0], &coords[1], &coords[2]};
5270 &coords[3], &coords[4], &coords[5]};
5273 t_dir(
i) = t_p1(
i) - t_p0(
i);
5279 auto evaluate_face_energy_and_set_orientation = [&](
auto front_edges,
5286 Tag th_material_force;
5288 case GRIFFITH_FORCE:
5289 case GRIFFITH_SKELETON:
5296 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5297 "Unknown energy release selector");
5305 auto find_maximal_face_energy = [&](
auto front_edges,
auto front_faces,
5306 auto &edge_face_max_energy_map) {
5315 for (
auto e : front_edges) {
5317 double griffith_force;
5323 faces = subtract(intersect(faces, front_faces), body_skin);
5324 std::vector<double> face_energy(faces.size());
5326 face_energy.data());
5327 auto max_energy_it =
5328 std::max_element(face_energy.begin(), face_energy.end());
5330 max_energy_it != face_energy.end() ? *max_energy_it : 0;
5332 edge_face_max_energy_map[e] =
5333 std::make_tuple(faces[max_energy_it - face_energy.begin()],
5334 griffith_force,
static_cast<double>(0));
5336 <<
"Edge " << e <<
" griffith force " << griffith_force
5337 <<
" max face energy " << max_energy <<
" factor "
5338 << max_energy / griffith_force;
5340 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
5362 auto calculate_face_orientation = [&](
auto &edge_face_max_energy_map) {
5365 auto up_down_face = [&](
5367 auto &face_angle_map_up,
5368 auto &face_angle_map_down
5373 for (
auto &
m : edge_face_max_energy_map) {
5375 auto [max_face, energy, opt_angle] =
m.second;
5379 faces = intersect(faces, front_faces);
5383 moab::Interface::UNION);
5384 if (adj_tets.size()) {
5389 moab::Interface::UNION);
5390 if (adj_tets.size()) {
5392 Range adj_tets_faces;
5395 adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
5396 moab::Interface::UNION);
5397 adj_tets_faces = intersect(adj_tets_faces, faces);
5402 get_cross(calculate_edge_direction(e,
true), max_face);
5403 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
5404 t_cross_max(
i) *= sense_max;
5406 for (
auto t : adj_tets) {
5407 Range adj_tets_faces;
5409 &
t, 1,
SPACE_DIM - 1,
false, adj_tets_faces);
5410 adj_tets_faces = intersect(adj_tets_faces, faces);
5412 subtract(adj_tets_faces,
Range(max_face, max_face));
5414 if (adj_tets_faces.size() == 1) {
5418 auto t_cross = get_cross(calculate_edge_direction(e,
true),
5420 auto [side, sense, offset] =
5421 get_sense(adj_tets_faces[0], e);
5422 t_cross(
i) *= sense;
5423 double dot = t_cross(
i) * t_cross_max(
i);
5424 auto angle = std::acos(dot);
5428 th_face_energy, adj_tets_faces, &face_energy);
5430 auto [side_face, sense_face, offset_face] =
5431 get_sense(
t, max_face);
5433 if (sense_face > 0) {
5434 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
5438 face_angle_map_down[e] = std::make_tuple(
5439 face_energy, -angle, adj_tets_faces[0]);
5450 auto calc_optimal_angle = [&](
5452 auto &face_angle_map_up,
5453 auto &face_angle_map_down
5458 for (
auto &
m : edge_face_max_energy_map) {
5460 auto &[max_face, e0,
a0] =
m.second;
5462 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
5464 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
5465 face_angle_map_down.find(e) == face_angle_map_down.end()) {
5470 case GRIFFITH_FORCE:
5471 case GRIFFITH_SKELETON: {
5473 Tag th_material_force;
5478 th_material_force, &e, 1, &t_material_force(0));
5479 auto material_force_magnitude = t_material_force.
l2();
5480 if (material_force_magnitude <
5481 std::numeric_limits<double>::epsilon()) {
5486 auto t_edge_dir = calculate_edge_direction(e,
true);
5487 auto t_cross_max = get_cross(t_edge_dir, max_face);
5488 auto [side, sense, offset] = get_sense(max_face, e);
5489 t_cross_max(sense) *= sense;
5496 t_cross_max.normalize();
5499 t_material_force(
J) * t_cross_max(
K);
5500 a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
5503 <<
"Optimal angle " <<
a0 <<
" energy " << e0;
5509 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5510 "Unknown energy release selector");
5520 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5522 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5523 face_angle_map_down;
5524 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
5525 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
5529 auto th_angle = get_tags_vec(
"Angle", 1);
5531 for (
auto &
m : face_angle_map_up) {
5532 auto [e,
a, face] =
m.second;
5537 for (
auto &
m : face_angle_map_down) {
5538 auto [e,
a, face] =
m.second;
5543 Range max_energy_faces;
5544 for (
auto &
m : edge_face_max_energy_map) {
5545 auto [face, e, angle] =
m.second;
5546 max_energy_faces.insert(face);
5562 auto get_conn = [&](
auto e) {
5569 auto get_adj = [&](
auto e,
auto dim) {
5572 e, dim,
false, adj, moab::Interface::UNION),
5577 auto get_coords = [&](
auto v) {
5585 auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
5588 auto t_edge_dir = calculate_edge_direction(e,
true);
5589 auto [side, sense, offset] = get_sense(
f, e);
5590 t_edge_dir(
i) *= sense;
5591 t_edge_dir.normalize();
5592 t_edge_dir(
i) *= angle;
5597 t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
5598 return std::make_tuple(t_normal, t_rotated_normal);
5601 auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
5602 auto &t_move,
auto gamma) {
5603 auto index = adj_vertex_tets_verts.index(
v);
5605 for (
auto ii : {0, 1, 2}) {
5606 coords[3 * index + ii] += gamma * t_move(ii);
5613 auto tets_quality = [&](
auto quality,
auto &adj_vertex_tets_verts,
5614 auto &adj_vertex_tets,
auto &coords) {
5615 for (
auto t : adj_vertex_tets) {
5619 std::array<double, 12> tet_coords;
5620 for (
auto n = 0;
n != 4; ++
n) {
5621 auto index = adj_vertex_tets_verts.index(conn[
n]);
5625 for (
auto ii = 0; ii != 3; ++ii) {
5626 tet_coords[3 *
n + ii] = coords[3 * index + ii];
5630 if (!std::isnormal(q))
5632 quality = std::min(quality, q);
5638 auto calculate_free_face_node_displacement =
5639 [&](
auto &edge_face_max_energy_map) {
5641 auto get_vertex_edges = [&](
auto vertex) {
5648 vertex_edges = subtract(vertex_edges, front_verts_edges);
5650 if (boundary_skin_verts.size() &&
5651 boundary_skin_verts.find(vertex[0]) !=
5652 boundary_skin_verts.end()) {
5653 MOFEM_LOG(
"EP", sev) <<
"Boundary vertex";
5654 vertex_edges = intersect(vertex_edges, body_skin_edges);
5656 if (geometry_edges_verts.size() &&
5657 geometry_edges_verts.find(vertex[0]) !=
5658 geometry_edges_verts.end()) {
5659 MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
5660 vertex_edges = intersect(vertex_edges, geometry_edges);
5662 if (crack_faces_verts.size() &&
5663 crack_faces_verts.find(vertex[0]) !=
5664 crack_faces_verts.end()) {
5665 MOFEM_LOG(
"EP", sev) <<
"Crack face vertex";
5666 vertex_edges = intersect(vertex_edges, crack_faces_edges);
5673 return vertex_edges;
5678 using Bundle = std::vector<
5684 std::map<EntityHandle, Bundle> edge_bundle_map;
5686 for (
auto &
m : edge_face_max_energy_map) {
5688 auto edge =
m.first;
5689 auto &[max_face, energy, opt_angle] =
m.second;
5692 auto [t_normal, t_rotated_normal] =
5693 get_rotated_normal(edge, max_face, opt_angle);
5695 auto front_vertex = get_conn(
Range(
m.first,
m.first));
5696 auto adj_tets = get_adj(
Range(max_face, max_face), 3);
5697 auto adj_tets_faces = get_adj(adj_tets, 2);
5698 auto adj_front_faces = subtract(
5699 intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
5701 if (adj_front_faces.size() > 3)
5703 "adj_front_faces.size()>3");
5707 &t_material_force(0));
5708 std::vector<double> griffith_energy(adj_front_faces.size());
5710 th_face_energy, adj_front_faces, griffith_energy.data());
5713 auto set_edge_bundle = [&](
auto min_gamma) {
5714 for (
auto rotated_f : adj_front_faces) {
5716 double rotated_face_energy =
5717 griffith_energy[adj_front_faces.index(rotated_f)];
5719 auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
5721 if (vertex.size() != 1) {
5723 "Wrong number of vertex to move");
5725 auto front_vertex_edges_vertex = get_conn(
5726 intersect(get_adj(front_vertex, 1), crack_faces_edges));
5728 vertex, front_vertex_edges_vertex);
5729 if (vertex.empty()) {
5733 auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
5736 subtract(body_skin_edges, crack_faces_edges));
5739 for (;
c < 10; ++
c) {
5741 subtract(get_adj(faces, 1), seen_front_edges);
5742 if (front_edges.size() == 0) {
5745 auto front_connected_edges =
5746 intersect(front_edges, whole_front);
5747 if (front_connected_edges.size()) {
5748 seen_front_edges.merge(front_connected_edges);
5751 faces.merge(get_adj(front_edges, 2));
5758 double rotated_face_cardinality = face_cardinality(
5764 rotated_face_cardinality = std::max(rotated_face_cardinality,
5767 auto t_vertex_coords = get_coords(vertex);
5768 auto vertex_edges = get_vertex_edges(vertex);
5777 for (
auto e_used_to_move_detection : vertex_edges) {
5778 auto edge_conn = get_conn(
Range(e_used_to_move_detection,
5779 e_used_to_move_detection));
5780 edge_conn = subtract(edge_conn, vertex);
5790 t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
5794 (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5796 (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5799 constexpr double eps =
5800 std::numeric_limits<double>::epsilon();
5801 if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
5804 t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
5806 auto check_rotated_face_directoon = [&]() {
5808 t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
5811 (t_material_force(
i) / t_material_force.
l2()) *
5813 return -dot > 0 ? true :
false;
5816 if (check_rotated_face_directoon()) {
5819 <<
"Crack edge " << edge <<
" moved face "
5821 <<
" edge: " << e_used_to_move_detection
5822 <<
" face direction/energy " << rotated_face_energy
5823 <<
" face cardinality " << rotated_face_cardinality
5824 <<
" gamma: " << gamma;
5826 auto &bundle = edge_bundle_map[edge];
5827 bundle.emplace_back(rotated_f, e_used_to_move_detection,
5828 vertex[0], t_move, 1,
5829 rotated_face_cardinality, gamma);
5836 set_edge_bundle(std::numeric_limits<double>::epsilon());
5837 if (edge_bundle_map[edge].empty()) {
5838 set_edge_bundle(-1.);
5842 return edge_bundle_map;
5845 auto get_sort_by_energy = [&](
auto &edge_face_max_energy_map) {
5846 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
5849 for (
auto &
m : edge_face_max_energy_map) {
5851 auto &[max_face, energy, opt_angle] =
m.second;
5852 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
5855 return sort_by_energy;
5858 auto set_tag = [&](
auto &&adj_edges_map,
auto &&sort_by_energy) {
5861 Tag th_face_pressure;
5863 mField.
get_moab().tag_get_handle(
"FacePressure", th_face_pressure),
5865 auto get_face_pressure = [&](
auto face) {
5874 <<
"Number of edges to check " << sort_by_energy.size();
5876 enum face_energy { POSITIVE, NEGATIVE };
5877 constexpr bool skip_negative =
true;
5879 for (
auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
5883 for (
auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
5886 auto energy = it->first;
5887 auto [max_edge, max_face, opt_angle] = it->second;
5889 auto face_pressure = get_face_pressure(max_face);
5890 if (skip_negative) {
5891 if (fe == face_energy::POSITIVE) {
5892 if (face_pressure < 0) {
5894 <<
"Skip negative face " << max_face <<
" with energy "
5895 << energy <<
" and pressure " << face_pressure;
5902 <<
"Check face " << max_face <<
" edge " << max_edge
5903 <<
" energy " << energy <<
" optimal angle " << opt_angle
5904 <<
" face pressure " << face_pressure;
5906 auto jt = adj_edges_map.find(max_edge);
5907 if (jt == adj_edges_map.end()) {
5909 <<
"Edge " << max_edge <<
" not found in adj_edges_map";
5912 auto &bundle = jt->second;
5914 auto find_max_in_bundle_impl = [&](
auto edge,
auto &bundle,
5921 double max_quality = -2;
5922 double max_quality_evaluated = -2;
5923 double min_cardinality = std::numeric_limits<double>::max();
5927 for (
auto &b : bundle) {
5928 auto &[face, move_edge, vertex, t_move, quality, cardinality,
5931 auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
5932 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
5933 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
5935 adj_vertex_tets_verts, coords.data()),
5938 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
5939 quality = tets_quality(quality, adj_vertex_tets_verts,
5940 adj_vertex_tets, coords);
5942 auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
5946 return ((edge_gamma < 0) ? (q / 2) : q) / pow(
c, 2);
5950 if (eval_quality(quality, cardinality, edge_gamma) >=
5951 max_quality_evaluated) {
5952 max_quality = quality;
5953 min_cardinality = cardinality;
5954 vertex_max = vertex;
5956 move_edge_max = move_edge;
5957 t_move_last(
i) = t_move(
i);
5958 max_quality_evaluated =
5959 eval_quality(max_quality, min_cardinality, edge_gamma);
5963 return std::make_tuple(vertex_max, face_max, t_move_last,
5964 max_quality, min_cardinality);
5967 auto find_max_in_bundle = [&](
auto edge,
auto &bundle) {
5968 auto b_org_bundle = bundle;
5969 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
5970 auto &[vertex_max, face_max, t_move_last, max_quality,
5972 if (max_quality < 0) {
5973 for (
double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
5974 bundle = b_org_bundle;
5975 r = find_max_in_bundle_impl(edge, bundle, gamma);
5976 auto &[vertex_max, face_max, t_move_last, max_quality,
5979 <<
"Back tracking: gamma " << gamma <<
" edge " << edge
5980 <<
" quality " << max_quality <<
" cardinality "
5982 if (max_quality > 0.01) {
5984 t_move_last(
I) *= gamma;
5995 auto set_tag_to_vertex_and_face = [&](
auto &&r,
auto &quality) {
5997 auto &[
v,
f, t_move, q, cardinality] = r;
5999 if ((q > 0 && std::isnormal(q)) && energy > 0) {
6002 <<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" "
6003 << max_edge <<
" move " << t_move <<
" energy " << energy
6004 <<
" quality " << q <<
" cardinality " << cardinality;
6015 double quality = -2;
6016 CHKERR set_tag_to_vertex_and_face(
6018 find_max_in_bundle(max_edge, bundle),
6024 if (quality > 0 && std::isnormal(quality) && energy > 0) {
6026 <<
"Crack face set with quality: " << quality;
6039 MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
6040 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
6041 edge_face_max_energy_map;
6042 CHKERR find_maximal_face_energy(front_edges, front_faces,
6043 edge_face_max_energy_map);
6044 CHKERR calculate_face_orientation(edge_face_max_energy_map);
6046 MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
6049 calculate_free_face_node_displacement(edge_face_max_energy_map),
6050 get_sort_by_energy(edge_face_max_energy_map)
6058 CHKERR evaluate_face_energy_and_set_orientation(
6059 *
frontEdges, get_adj_front(
false), sides_pair, th_front_position);
6077 auto get_max_moved_faces = [&]() {
6078 Range max_moved_faces;
6079 auto adj_front = get_adj_front(
false);
6080 std::vector<double> face_energy(adj_front.size());
6082 face_energy.data());
6083 for (
int i = 0;
i != adj_front.size(); ++
i) {
6084 if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
6085 max_moved_faces.insert(adj_front[
i]);
6089 return boost::make_shared<Range>(max_moved_faces);
6100 "max_moved_faces_" +
6115 Tag th_front_position;
6117 mField.
get_moab().tag_get_handle(
"FrontPosition", th_front_position);
6122 std::vector<double> coords(3 * verts.size());
6124 std::vector<double> pos(3 * verts.size());
6126 for (
int i = 0;
i != 3 * verts.size(); ++
i) {
6127 coords[
i] += pos[
i];
6130 double zero[] = {0., 0., 0.};
6135 constexpr bool debug =
false;
6140 "set_coords_faces_" +
6151 constexpr bool potential_crack_debug =
false;
6152 if constexpr (potential_crack_debug) {
6155 Range crack_front_verts;
6160 Range crack_front_faces;
6162 true, crack_front_faces,
6163 moab::Interface::UNION);
6164 crack_front_faces = intersect(crack_front_faces, add_ents);
6171 auto get_crack_faces = [&]() {
6179 auto get_extended_crack_faces = [&]() {
6180 auto get_faces_of_crack_front_verts = [&](
auto crack_faces_org) {
6181 ParallelComm *pcomm =
6186 if (!pcomm->rank()) {
6188 auto get_nodes = [&](
auto &&e) {
6191 "get connectivity");
6195 auto get_adj = [&](
auto &&e,
auto dim,
6196 auto t = moab::Interface::UNION) {
6208 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6211 auto front_block_nodes = get_nodes(front_block_edges);
6215 s = crack_faces.size();
6217 auto crack_face_nodes = get_nodes(crack_faces_org);
6218 auto crack_faces_edges =
6219 get_adj(crack_faces_org, 1, moab::Interface::UNION);
6222 front_block_edges = subtract(front_block_edges, crack_skin);
6223 auto crack_skin_nodes = get_nodes(crack_skin);
6224 crack_skin_nodes.merge(front_block_nodes);
6226 auto crack_skin_faces =
6227 get_adj(crack_skin, 2, moab::Interface::UNION);
6229 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
6231 crack_faces = crack_faces_org;
6232 for (
auto f : crack_skin_faces) {
6233 auto edges = intersect(
6234 get_adj(
Range(
f,
f), 1, moab::Interface::UNION), crack_skin);
6238 if (edges.size() == 2) {
6240 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6244 if (edges.size() == 2) {
6245 auto edge_conn = get_nodes(
Range(edges));
6246 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
6248 if (faces.size() == 2) {
6249 auto edge0_conn = get_nodes(
Range(edges[0], edges[0]));
6250 auto edge1_conn = get_nodes(
Range(edges[1], edges[1]));
6251 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
6257 moab::Interface::INTERSECT),
6262 if (node_edges.size()) {
6267 auto get_t_dir = [&](
auto e_conn) {
6268 auto other_node = subtract(e_conn,
edges_conn);
6272 t_dir(
i) -= t_v0(
i);
6278 get_t_dir(edge0_conn)(
i) + get_t_dir(edge1_conn)(
i);
6281 t_crack_surface_ave_dir(
i) = 0;
6282 for (
auto e : node_edges) {
6283 auto e_conn = get_nodes(
Range(e, e));
6284 auto t_dir = get_t_dir(e_conn);
6285 t_crack_surface_ave_dir(
i) += t_dir(
i);
6288 auto dot = t_ave_dir(
i) * t_crack_surface_ave_dir(
i);
6291 if (dot < -std::numeric_limits<double>::epsilon()) {
6292 crack_faces.insert(
f);
6295 crack_faces.insert(
f);
6299 }
else if (edges.size() == 3) {
6300 crack_faces.insert(
f);
6304 if (edges.size() == 1) {
6306 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6309 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6310 front_block_edges));
6311 if (edges.size() == 2) {
6312 crack_faces.insert(
f);
6318 crack_faces_org = crack_faces;
6320 }
while (s != crack_faces.size());
6326 return get_faces_of_crack_front_verts(get_crack_faces());
6331 get_extended_crack_faces());
6334 auto reconstruct_crack_faces = [&](
auto crack_faces) {
6335 ParallelComm *pcomm =
6341 Range new_crack_faces;
6342 if (!pcomm->rank()) {
6344 auto get_nodes = [&](
auto &&e) {
6347 "get connectivity");
6351 auto get_adj = [&](
auto &&e,
auto dim,
6352 auto t = moab::Interface::UNION) {
6360 auto get_test_on_crack_surface = [&]() {
6361 auto crack_faces_nodes =
6362 get_nodes(crack_faces);
6363 auto crack_faces_tets =
6364 get_adj(crack_faces_nodes, 3,
6365 moab::Interface::UNION);
6369 auto crack_faces_tets_nodes =
6370 get_nodes(crack_faces_tets);
6371 crack_faces_tets_nodes =
6372 subtract(crack_faces_tets_nodes, crack_faces_nodes);
6374 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
6375 moab::Interface::UNION));
6377 get_adj(crack_faces_tets, 2,
6378 moab::Interface::UNION);
6380 new_crack_faces.merge(crack_faces);
6382 return std::make_tuple(new_crack_faces, crack_faces_tets);
6385 auto carck_faces_test_edges = [&](
auto faces,
auto tets) {
6386 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
6387 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
6388 moab::Interface::UNION);
6389 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
6392 adj_faces_edges.merge(geometry_edges);
6393 adj_faces_edges.merge(front_block_edges);
6395 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
6396 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
6397 auto boundary_test_nodes_edges =
6398 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
6399 auto boundary_test_nodes_edges_nodes = subtract(
6400 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
6402 boundary_tets_edges =
6403 subtract(boundary_test_nodes_edges,
6404 get_adj(boundary_test_nodes_edges_nodes, 1,
6405 moab::Interface::UNION));
6412 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6413 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
6415 body_skin = intersect(body_skin, adj_tets_faces);
6416 body_skin_edges = subtract(
6417 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
6420 for (
auto e : body_skin_edges) {
6421 auto adj_tet = intersect(
6422 get_adj(
Range(e, e), 3, moab::Interface::INTERSECT), tets);
6423 if (adj_tet.size() == 1) {
6424 boundary_tets_edges.insert(e);
6428 return boundary_tets_edges;
6431 auto p = get_test_on_crack_surface();
6432 auto &[new_crack_faces, crack_faces_tets] = p;
6443 auto boundary_tets_edges =
6444 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
6446 boundary_tets_edges);
6448 auto resolve_surface = [&](
auto boundary_tets_edges,
6449 auto crack_faces_tets) {
6450 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
6451 auto crack_faces_tets_faces =
6452 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6454 Range all_removed_faces;
6455 Range all_removed_tets;
6459 while (size != crack_faces_tets.size()) {
6461 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6465 auto skin_skin_nodes = get_nodes(skin_skin);
6467 size = crack_faces_tets.size();
6469 <<
"Crack faces tets size " << crack_faces_tets.size()
6470 <<
" crack faces size " << crack_faces_tets_faces.size();
6471 auto skin_tets_nodes = subtract(
6472 get_nodes(skin_tets),
6473 boundary_tets_edges_nodes);
6475 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
6477 Range removed_nodes;
6478 Range tets_to_remove;
6479 Range faces_to_remove;
6480 for (
auto n : skin_tets_nodes) {
6482 intersect(get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT),
6484 if (tets.size() == 0) {
6488 auto hole_detetction = [&]() {
6490 get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT);
6495 if (adj_tets.size() == 0) {
6496 return std::make_pair(
6498 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6503 std::vector<Range> tets_groups;
6504 auto test_adj_tets = adj_tets;
6505 while (test_adj_tets.size()) {
6507 Range seed =
Range(test_adj_tets[0], test_adj_tets[0]);
6508 while (seed.size() != seed_size) {
6510 subtract(get_adj(seed, 2, moab::Interface::UNION),
6513 seed_size = seed.size();
6515 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
6518 tets_groups.push_back(seed);
6519 test_adj_tets = subtract(test_adj_tets, seed);
6521 if (tets_groups.size() == 1) {
6523 return std::make_pair(
6525 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6531 Range tets_to_remove;
6532 Range faces_to_remove;
6533 for (
auto &r : tets_groups) {
6534 auto f = get_adj(r, 2, moab::Interface::UNION);
6535 auto t = intersect(get_adj(
f, 3, moab::Interface::UNION),
6538 if (
f.size() > faces_to_remove.size() ||
6539 faces_to_remove.size() == 0) {
6540 faces_to_remove =
f;
6545 <<
"Hole detection: faces to remove "
6546 << faces_to_remove.size() <<
" tets to remove "
6547 << tets_to_remove.size();
6548 return std::make_pair(faces_to_remove, tets_to_remove);
6551 if (tets.size() < tets_to_remove.size() ||
6552 tets_to_remove.size() == 0) {
6554 auto [h_faces_to_remove, h_tets_to_remove] =
6556 faces_to_remove = h_faces_to_remove;
6557 tets_to_remove = h_tets_to_remove;
6565 all_removed_faces.merge(faces_to_remove);
6566 all_removed_tets.merge(tets_to_remove);
6569 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
6570 crack_faces_tets_faces =
6571 subtract(crack_faces_tets_faces, faces_to_remove);
6576 boost::lexical_cast<std::string>(counter) +
".vtk",
6579 "faces_to_remove_" +
6580 boost::lexical_cast<std::string>(counter) +
".vtk",
6584 boost::lexical_cast<std::string>(counter) +
".vtk",
6587 "crack_faces_tets_faces_" +
6588 boost::lexical_cast<std::string>(counter) +
".vtk",
6589 crack_faces_tets_faces);
6591 "crack_faces_tets_" +
6592 boost::lexical_cast<std::string>(counter) +
".vtk",
6598 auto cese_internal_faces = [&]() {
6601 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
6603 subtract(adj_faces, skin_tets);
6604 auto adj_tets = get_adj(adj_faces, 3,
6605 moab::Interface::UNION);
6608 subtract(crack_faces_tets,
6611 crack_faces_tets_faces =
6612 subtract(crack_faces_tets_faces, adj_faces);
6614 all_removed_faces.merge(adj_faces);
6615 all_removed_tets.merge(adj_tets);
6619 <<
"Remove internal faces size " << adj_faces.size()
6620 <<
" tets size " << adj_tets.size();
6624 auto case_only_one_free_edge = [&]() {
6627 for (
auto t :
Range(crack_faces_tets)) {
6629 auto adj_faces = get_adj(
6631 moab::Interface::UNION);
6632 auto crack_surface_edges =
6633 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6636 moab::Interface::UNION);
6639 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6640 crack_surface_edges);
6641 adj_edges = subtract(
6643 boundary_tets_edges);
6645 if (adj_edges.size() == 1) {
6647 subtract(crack_faces_tets,
6651 auto faces_to_remove =
6652 get_adj(adj_edges, 2, moab::Interface::UNION);
6655 crack_faces_tets_faces =
6656 subtract(crack_faces_tets_faces, faces_to_remove);
6658 all_removed_faces.merge(faces_to_remove);
6659 all_removed_tets.merge(
Range(
t,
t));
6661 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Remove free one edges ";
6665 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6666 crack_faces_tets_faces =
6667 subtract(crack_faces_tets_faces, all_removed_faces);
6672 auto cese_flat_tet = [&](
auto max_adj_edges) {
6679 auto body_skin_edges =
6680 get_adj(body_skin, 1, moab::Interface::UNION);
6682 for (
auto t :
Range(crack_faces_tets)) {
6684 auto adj_faces = get_adj(
6686 moab::Interface::UNION);
6687 auto crack_surface_edges =
6688 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6691 moab::Interface::UNION);
6694 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6695 crack_surface_edges);
6696 adj_edges = subtract(adj_edges, body_skin_edges);
6698 auto tet_edges = get_adj(
Range(
t,
t), 1,
6699 moab::Interface::UNION);
6701 tet_edges = subtract(tet_edges, adj_edges);
6703 for (
auto e : tet_edges) {
6704 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
6705 auto get_side = [&](
auto e) {
6706 int side, sense, offset;
6709 "get side number failed");
6712 auto get_side_ent = [&](
auto side) {
6719 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
6722 if (adj_edges.size() <= max_adj_edges) {
6725 Range faces_to_remove;
6726 for (
auto e : adj_edges) {
6727 auto edge_adj_faces =
6728 get_adj(
Range(e, e), 2, moab::Interface::UNION);
6729 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
6730 if (edge_adj_faces.size() != 2) {
6732 "Adj faces size is not 2 for edge " +
6733 boost::lexical_cast<std::string>(e));
6736 auto get_normal = [&](
auto f) {
6740 "get tri normal failed");
6743 auto t_n0 = get_normal(edge_adj_faces[0]);
6744 auto t_n1 = get_normal(edge_adj_faces[1]);
6745 auto get_sense = [&](
auto f) {
6746 int side, sense, offset;
6749 "get side number failed");
6752 auto sense0 = get_sense(edge_adj_faces[0]);
6753 auto sense1 = get_sense(edge_adj_faces[1]);
6758 auto dot_e = (sense0 * sense1) * t_n0(
i) * t_n1(
i);
6759 if (dot_e < dot || e == adj_edges[0]) {
6761 faces_to_remove = edge_adj_faces;
6765 all_removed_faces.merge(faces_to_remove);
6766 all_removed_tets.merge(
Range(
t,
t));
6769 <<
"Remove free edges on flat tet, with considered nb. of "
6771 << adj_edges.size();
6775 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6776 crack_faces_tets_faces =
6777 subtract(crack_faces_tets_faces, all_removed_faces);
6783 "Case only one free edge failed");
6784 for (
auto max_adj_edges : {0, 1, 2, 3}) {
6786 "Case only one free edge failed");
6789 "Case internal faces failed");
6793 "crack_faces_tets_faces_" +
6794 boost::lexical_cast<std::string>(counter) +
".vtk",
6795 crack_faces_tets_faces);
6797 "crack_faces_tets_" +
6798 boost::lexical_cast<std::string>(counter) +
".vtk",
6802 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
6803 all_removed_faces, all_removed_tets);
6806 auto [resolved_faces, resolved_tets, all_removed_faces,
6808 resolve_surface(boundary_tets_edges, crack_faces_tets);
6809 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
6817 crack_faces = resolved_faces;
6829 auto resolve_consisten_crack_extension = [&]() {
6831 auto crack_meshset =
6834 auto meshset = crack_meshset->getMeshset();
6838 Range old_crack_faces;
6841 auto extendeded_crack_faces = get_extended_crack_faces();
6842 auto reconstructed_crack_faces =
6843 subtract(reconstruct_crack_faces(extendeded_crack_faces),
6845 if (
nbCrackFaces >= reconstructed_crack_faces.size()) {
6847 <<
"No new crack faces to add, skipping adding to meshset";
6848 extendeded_crack_faces = subtract(
6849 extendeded_crack_faces, subtract(*
crackFaces, old_crack_faces));
6851 <<
"Number crack faces size (extended) "
6852 << extendeded_crack_faces.size();
6858 reconstructed_crack_faces);
6860 <<
"Number crack faces size (reconstructed) "
6861 << reconstructed_crack_faces.size();
6880 CHKERR resolve_consisten_crack_extension();
6893 verts = subtract(verts, conn);
6894 std::vector<double> coords(3 * verts.size());
6896 double def_coords[] = {0., 0., 0.};
6899 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
6900 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
6920 auto post_proc_norm_fe =
6921 boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
6923 auto post_proc_norm_rule_hook = [](int, int,
int p) ->
int {
6926 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
6928 post_proc_norm_fe->getUserPolynomialBase() =
6931 CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
6935 enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
6938 CHKERR VecZeroEntries(norms_vec);
6940 auto u_l2_ptr = boost::make_shared<MatrixDouble>();
6941 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
6942 post_proc_norm_fe->getOpPtrVector().push_back(
6944 post_proc_norm_fe->getOpPtrVector().push_back(
6946 post_proc_norm_fe->getOpPtrVector().push_back(
6948 post_proc_norm_fe->getOpPtrVector().push_back(
6950 post_proc_norm_fe->getOpPtrVector().push_back(
6954 auto piola_ptr = boost::make_shared<MatrixDouble>();
6955 post_proc_norm_fe->getOpPtrVector().push_back(
6957 post_proc_norm_fe->getOpPtrVector().push_back(
6961 post_proc_norm_fe->getOpPtrVector().push_back(
6964 TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
6966 *post_proc_norm_fe);
6967 TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
6969 CHKERR VecAssemblyBegin(norms_vec);
6970 CHKERR VecAssemblyEnd(norms_vec);
6971 const double *norms;
6972 CHKERR VecGetArrayRead(norms_vec, &norms);
6973 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
6974 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
6976 <<
"norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
6978 <<
"norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
6979 CHKERR VecRestoreArrayRead(norms_vec, &norms);
6993 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6994 if (
auto disp_bc = bc.second->dispBcPtr) {
6999 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7000 MOFEM_LOG(
"EP", Sev::noisy) <<
"Displacement BC: " << *disp_bc;
7002 std::vector<double> block_attributes(6, 0.);
7003 if (disp_bc->data.flag1 == 1) {
7004 block_attributes[0] = disp_bc->data.value1;
7005 block_attributes[3] = 1;
7007 if (disp_bc->data.flag2 == 1) {
7008 block_attributes[1] = disp_bc->data.value2;
7009 block_attributes[4] = 1;
7011 if (disp_bc->data.flag3 == 1) {
7012 block_attributes[2] = disp_bc->data.value3;
7013 block_attributes[5] = 1;
7015 auto faces = bc.second->bcEnts.subset_by_dimension(2);
7023 boost::make_shared<NormalDisplacementBcVec>();
7024 for (
auto bc : bc_mng->getBcMapByBlockName()) {
7025 auto block_name =
"(.*)NORMAL_DISPLACEMENT(.*)";
7026 std::regex reg_name(block_name);
7027 if (std::regex_match(bc.first, reg_name)) {
7031 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7033 block_name, bc.second->bcAttributes,
7034 bc.second->bcEnts.subset_by_dimension(2));
7039 boost::make_shared<AnalyticalDisplacementBcVec>();
7041 for (
auto bc : bc_mng->getBcMapByBlockName()) {
7042 auto block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
7043 std::regex reg_name(block_name);
7044 if (std::regex_match(bc.first, reg_name)) {
7048 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7050 block_name, bc.second->bcAttributes,
7051 bc.second->bcEnts.subset_by_dimension(2));
7055 auto ts_displacement =
7056 boost::make_shared<DynamicRelaxationTimeScale>(
"disp_history.txt");
7059 <<
"Add time scaling displacement BC: " << bc.blockName;
7062 ts_displacement,
"disp_history",
".txt", bc.blockName);
7065 auto ts_normal_displacement =
7066 boost::make_shared<DynamicRelaxationTimeScale>(
"normal_disp_history.txt");
7069 <<
"Add time scaling normal displacement BC: " << bc.blockName;
7072 ts_normal_displacement,
"normal_disp_history",
".txt",
7088 for (
auto bc : bc_mng->getBcMapByBlockName()) {
7089 if (
auto force_bc = bc.second->forceBcPtr) {
7094 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7095 MOFEM_LOG(
"EP", Sev::noisy) <<
"Force BC: " << *force_bc;
7097 std::vector<double> block_attributes(6, 0.);
7098 block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
7099 block_attributes[3] = 1;
7100 block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
7101 block_attributes[4] = 1;
7102 block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
7103 block_attributes[5] = 1;
7104 auto faces = bc.second->bcEnts.subset_by_dimension(2);
7112 for (
auto bc : bc_mng->getBcMapByBlockName()) {
7113 auto block_name =
"(.*)PRESSURE(.*)";
7114 std::regex reg_name(block_name);
7115 if (std::regex_match(bc.first, reg_name)) {
7120 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7122 block_name, bc.second->bcAttributes,
7123 bc.second->bcEnts.subset_by_dimension(2));
7128 boost::make_shared<AnalyticalTractionBcVec>();
7130 for (
auto bc : bc_mng->getBcMapByBlockName()) {
7131 auto block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
7132 std::regex reg_name(block_name);
7133 if (std::regex_match(bc.first, reg_name)) {
7137 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
7139 block_name, bc.second->bcAttributes,
7140 bc.second->bcEnts.subset_by_dimension(2));
7145 boost::make_shared<DynamicRelaxationTimeScale>(
"traction_history.txt");
7149 ts_traction,
"traction_history",
".txt", bc.blockName);
7153 boost::make_shared<DynamicRelaxationTimeScale>(
"pressure_history.txt");
7157 ts_pressure,
"pressure_history",
".txt", bc.blockName);
7166 auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
7167 const std::string block_name,
7168 const int nb_attributes) {
7171 std::regex((boost::format(
"(.*)%s(.*)") % block_name).str()))
7173 std::vector<double> block_attributes;
7174 CHKERR it->getAttributes(block_attributes);
7175 if (block_attributes.size() < nb_attributes) {
7177 "In block %s expected %d attributes, but given %ld",
7178 it->getName().c_str(), nb_attributes, block_attributes.size());
7181 auto get_block_ents = [&]() {
7187 auto Ents = get_block_ents();
7188 ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
7198 auto ts_pre_stretch =
7199 boost::make_shared<DynamicRelaxationTimeScale>(
"externalstrain_history.txt");
7202 <<
"Add time scaling external strain: " << ext_strain_block.blockName;
7205 ts_pre_stretch,
"externalstrain_history",
".txt", ext_strain_block.blockName);
7214 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
7217 CHKERR VecGetLocalSize(
v.second, &size);
7219 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
7220 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( " << low
7221 <<
" " << high <<
" ) ";
7253 MOFEM_LOG(
"EP", Sev::inform) <<
"Calculate crack area";
7255 for (
auto f : crack_faces) {
7258 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
7260 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
7262 if (success != MPI_SUCCESS) {
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Implementation of CGGUserPolynomialBase class.
Auxilary functions for Eshelbian plasticity.
Contains definition of EshelbianMonitor class.
FormsIntegrators< FaceElementForcesAndSourcesCore::UserDataOperator >::Assembly< A >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassVectorFace
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
static auto get_block_meshset(MoFEM::Interface &m_field, const int ms_id, const unsigned int cubit_bc_type)
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
static auto get_range_from_block_map(MoFEM::Interface &m_field, const std::string block_name, int dim)
static auto filter_owners(MoFEM::Interface &m_field, Range skin)
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
static auto get_crack_front_edges(MoFEM::Interface &m_field, Range crack_faces)
Eshelbian plasticity interface.
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
#define FTENSOR_INDEX(DIM, I)
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
constexpr int SPACE_DIM
[Define dimension]
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Tensor1< T, Tensor_Dim > normalize()
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
Order displacement.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
Add CUBIT meshset to manager.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark DOFs on block entities for boundary conditions.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
auto id_from_handle(const EntityHandle h)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
static constexpr int edges_conn[]
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
constexpr IntegrationType I
constexpr AssemblyType A
[Define dimension]
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
constexpr auto field_name
#define QUAD_2D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_3D_TABLE[]
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
FTensor::Index< 'm', 3 > m
CGG User Polynomial Base.
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static enum StretchSelector stretchSelector
boost::shared_ptr< Range > frontAdjEdges
MoFEMErrorCode createCrackSurfaceMeshset()
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
const std::string skeletonElement
static double inv_f_linear(const double v)
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
boost::shared_ptr< Range > contactFaces
static double dd_f_log_e_quadratic(const double v)
static double dd_f_log(const double v)
BitRefLevel bitAdjEnt
bit ref level for parent
static boost::function< double(const double)> inv_dd_f
MoFEM::Interface & mField
const std::string spatialL2Disp
static double inv_d_f_log(const double v)
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
static PetscBool l2UserBaseScale
SmartPetscObj< DM > dM
Coupled problem all fields.
static int internalStressInterpOrder
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
const std::string materialH1Positions
static int nbJIntegralContours
MoFEMErrorCode setBlockTagsOnSkin()
static PetscBool crackingOn
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
Remove all, but entities where kinematic constrains are applied.
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates, SmartPetscObj< Vec > ver_vec, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe)
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
static double griffithEnergy
Griffith energy.
MoFEMErrorCode calculateCrackArea(boost::shared_ptr< double > area_ptr)
const std::string elementVolumeName
static double dd_f_log_e(const double v)
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
CommInterface::EntitiesPetscVector vertexExchange
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
boost::shared_ptr< Range > maxMovedFaces
static PetscBool dynamicRelaxation
const std::string spatialH1Disp
MoFEMErrorCode solveElastic(TS ts, Vec x)
static double d_f_log(const double v)
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
static double crackingStartTime
static enum MaterialModel materialModel
MoFEMErrorCode getOptions()
const std::string piolaStress
MoFEMErrorCode setElasticElementToTs(DM dm)
static double inv_d_f_log_e(const double v)
int contactRefinementLevels
MoFEMErrorCode gettingNorms()
[Getting norms]
MoFEMErrorCode setVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
const std::string bubbleField
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation)
static double inv_f_log(const double v)
static PetscBool noStretch
MoFEMErrorCode setNewFrontCoordinates()
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
static double exponentBase
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
static double dd_f_linear(const double v)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
boost::shared_ptr< Range > skeletonFaces
boost::shared_ptr< PhysicalEquations > physicalEquations
const std::string rotAxis
static double inv_d_f_linear(const double v)
BitRefLevel bitAdjParentMask
bit ref level for parent parent
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x, int start_step, double start_time)
Solve problem using dynamic relaxation method.
static double inv_dd_f_log(const double v)
const std::string contactDisp
static std::string internalStressTagName
static enum SymmetrySelector symmetrySelector
CommInterface::EntitiesPetscVector edgeExchange
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
static boost::function< double(const double)> f
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
const std::string skinElement
static PetscBool internalStressVoigt
static double inv_dd_f_linear(const double v)
static double inv_dd_f_log_e(const double v)
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
static PetscBool setSingularity
static double d_f_log_e(const double v)
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
static double f_log_e_quadratic(const double v)
MoFEMErrorCode addCrackSurfaces(const bool debug=false)
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
BitRefLevel bitAdjParent
bit ref level for parent
MoFEMErrorCode postProcessResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, Vec var_vec=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
static double d_f_log_e_quadratic(const double v)
CommInterface::EntitiesPetscVector volumeExchange
MoFEMErrorCode saveOrgCoords()
const std::string naturalBcElement
static boost::function< double(const double)> dd_f
static double f_log_e(const double v)
static int addCrackMeshsetId
static double inv_f_log_e(const double v)
MoFEMErrorCode createExchangeVectors(Sev sev)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< Range > crackFaces
static boost::function< double(const double)> d_f
boost::shared_ptr< Range > frontVertices
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_d_f
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
static double d_f_linear(const double v)
const std::string hybridSpatialDisp
SmartPetscObj< Vec > solTSStep
static double f_log(const double v)
CommInterface::EntitiesPetscVector faceExchange
SmartPetscObj< DM > dmElastic
Elastic problem.
EshelbianCore(MoFEM::Interface &m_field)
boost::shared_ptr< Range > frontEdges
static boost::function< double(const double)> inv_f
const std::string stretchTensor
BitRefLevel bitAdjEntMask
bit ref level for parent parent
static double f_linear(const double v)
const std::string contactElement
AnalyticalDisplacementBc(std::string name, std::vector< double > attr, Range faces)
AnalyticalTractionBc(std::string name, std::vector< double > attr, Range faces)
BcRot(std::string name, std::vector< double > attr, Range faces)
ExternalStrain(std::string name, std::vector< double > attr, Range ents)
int operator()(int p_row, int p_col, int p_data) const
NormalDisplacementBc(std::string name, std::vector< double > attr, Range faces)
PressureBc(std::string name, std::vector< double > attr, Range faces)
boost::shared_ptr< Range > frontNodes
SetIntegrationAtFrontFace(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges)
boost::shared_ptr< Range > frontEdges
boost::function< int(int)> FunRule
SetIntegrationAtFrontFace(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges, FunRule fun_rule)
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
static std::map< long int, MatrixDouble > mapRefCoords
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
boost::shared_ptr< Range > frontNodes
static std::map< long int, MatrixDouble > mapRefCoords
boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cachePhi
SetIntegrationAtFrontVolume(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges, boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cache_phi=nullptr)
boost::shared_ptr< Range > frontEdges
TractionBc(std::string name, std::vector< double > attr, Range faces)
Set integration rule on element.
int operator()(int p_row, int p_col, int p_data) const
static auto setup(EshelbianCore *ep_ptr, TS ts, Vec x, bool set_ts_monitor)
static auto exp(A &&t_w_vee, B &&theta)
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Boundary condition manager for finite element problem setup.
static std::pair< std::string, std::string > extractStringFromBlockId(const std::string block_id, const std::string prb_name)
Extract block name and block name from block id.
static MoFEMErrorCode updateEntitiesPetscVector(moab::Interface &moab, EntitiesPetscVector &vec, Tag tag, UpdateGhosts update_gosts=defaultUpdateGhosts)
Exchange data between vector and data.
static Range getPartEntities(moab::Interface &moab, int part)
static EntitiesPetscVector createEntitiesPetscVector(MPI_Comm comm, moab::Interface &moab, int dim, const int nb_coeffs, Sev sev=Sev::verbose, int root_rank=0)
Create a ghost vector for exchanging data.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)> > > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
Structure for user loop methods on finite elements.
EntityHandle getFEEntityHandle() const
Get the entity handle of the current finite element.
default operator for TRI element
Field data structure for finite element approximation.
Definition of the force bc data structure.
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Constructor for operators working on finite element spaces.
structure to get information from mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
static boost::shared_ptr< ScalingMethod > get(boost::shared_ptr< ScalingMethod > ts, std::string file_prefix, std::string file_suffix, std::string block_name, Args &&...args)
Section manager is used to create indexes and sections.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Operator for broken loop side.
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients time derivative at integration pts for scalar field rank 0, i....
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Execute "this" element in the operator.
Post post-proc data at points from hash maps.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Template struct for dimension-specific finite element types.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
default operator for TET element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Apply rotation boundary condition.
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)
BoundaryEle::UserDataOperator BdyEleOp
ElementsAndOps< SPACE_DIM >::SideEle SideEle