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"
48 using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
52 using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
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_edges)
370 int order_col,
int order_data) {
373 constexpr bool debug =
false;
375 constexpr int numNodes = 4;
376 constexpr int numEdges = 6;
377 constexpr int refinementLevels = 4;
379 auto &m_field = fe_raw_ptr->
mField;
380 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
383 auto set_base_quadrature = [&]() {
385 int rule = 2 * order_data + 1;
396 auto &gauss_pts = fe_ptr->gaussPts;
397 gauss_pts.resize(4, nb_gauss_pts,
false);
398 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
399 &gauss_pts(0, 0), 1);
400 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
401 &gauss_pts(1, 0), 1);
402 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
403 &gauss_pts(2, 0), 1);
405 &gauss_pts(3, 0), 1);
406 auto &data = fe_ptr->dataOnElement[
H1];
407 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
410 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
411 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
420 CHKERR set_base_quadrature();
424 auto get_singular_nodes = [&]() {
427 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
429 std::bitset<numNodes> singular_nodes;
430 for (
auto nn = 0; nn != numNodes; ++nn) {
432 singular_nodes.set(nn);
434 singular_nodes.reset(nn);
437 return singular_nodes;
440 auto get_singular_edges = [&]() {
441 std::bitset<numEdges> singular_edges;
442 for (
int ee = 0; ee != numEdges; ee++) {
444 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
446 singular_edges.set(ee);
448 singular_edges.reset(ee);
451 return singular_edges;
454 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
456 fe_ptr->gaussPts.swap(ref_gauss_pts);
457 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
458 auto &data = fe_ptr->dataOnElement[
H1];
459 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
461 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
463 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
468 auto singular_nodes = get_singular_nodes();
469 if (singular_nodes.count()) {
470 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
472 CHKERR set_gauss_pts(it_map_ref_coords->second);
476 auto refine_quadrature = [&]() {
479 const int max_level = refinementLevels;
483 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
485 for (
int nn = 0; nn != 4; nn++)
486 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
487 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
491 Range tets(tet, tet);
494 tets, 1,
true, edges, moab::Interface::UNION);
499 Range nodes_at_front;
500 for (
int nn = 0; nn != numNodes; nn++) {
501 if (singular_nodes[nn]) {
503 CHKERR moab_ref.side_element(tet, 0, nn, ent);
504 nodes_at_front.insert(ent);
508 auto singular_edges = get_singular_edges();
511 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
512 for (
int ee = 0; ee != numEdges; ee++) {
513 if (singular_edges[ee]) {
515 CHKERR moab_ref.side_element(tet, 1, ee, ent);
516 CHKERR moab_ref.add_entities(meshset, &ent, 1);
522 for (
int ll = 0; ll != max_level; ll++) {
525 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
529 CHKERR moab_ref.get_adjacencies(
530 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
531 ref_edges = intersect(ref_edges, edges);
533 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
534 ref_edges = intersect(ref_edges, ents);
537 ->getEntitiesByTypeAndRefLevel(
539 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
543 ->updateMeshsetByEntitiesChildren(meshset,
545 meshset, MBEDGE,
true);
551 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
561 for (Range::iterator tit = tets.begin(); tit != tets.end();
565 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
566 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
569 auto &data = fe_ptr->dataOnElement[
H1];
570 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
571 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
573 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
575 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
576 double *tet_coords = &ref_coords(tt, 0);
579 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
580 for (
int dd = 0; dd != 3; dd++) {
581 ref_gauss_pts(dd, gg) =
582 shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
583 shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
584 shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
585 shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
587 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
591 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
597 CHKERR refine_quadrature();
607 using ForcesAndSourcesCore::dataOnElement;
610 using ForcesAndSourcesCore::ForcesAndSourcesCore;
625 boost::shared_ptr<Range> front_edges)
629 boost::shared_ptr<Range> front_edges,
634 int order_col,
int order_data) {
637 constexpr bool debug =
false;
639 constexpr int numNodes = 3;
640 constexpr int numEdges = 3;
641 constexpr int refinementLevels = 4;
643 auto &m_field = fe_raw_ptr->
mField;
644 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
647 auto set_base_quadrature = [&]() {
649 int rule =
funRule(order_data);
659 fe_ptr->gaussPts.resize(3, nb_gauss_pts,
false);
660 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
661 &fe_ptr->gaussPts(0, 0), 1);
662 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
663 &fe_ptr->gaussPts(1, 0), 1);
665 &fe_ptr->gaussPts(2, 0), 1);
666 auto &data = fe_ptr->dataOnElement[
H1];
667 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
670 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
671 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
673 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
676 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
685 CHKERR set_base_quadrature();
689 auto get_singular_nodes = [&]() {
692 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
694 std::bitset<numNodes> singular_nodes;
695 for (
auto nn = 0; nn != numNodes; ++nn) {
697 singular_nodes.set(nn);
699 singular_nodes.reset(nn);
702 return singular_nodes;
705 auto get_singular_edges = [&]() {
706 std::bitset<numEdges> singular_edges;
707 for (
int ee = 0; ee != numEdges; ee++) {
709 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
711 singular_edges.set(ee);
713 singular_edges.reset(ee);
716 return singular_edges;
719 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
721 fe_ptr->gaussPts.swap(ref_gauss_pts);
722 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
723 auto &data = fe_ptr->dataOnElement[
H1];
724 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
726 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
728 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
732 auto singular_nodes = get_singular_nodes();
733 if (singular_nodes.count()) {
734 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
736 CHKERR set_gauss_pts(it_map_ref_coords->second);
740 auto refine_quadrature = [&]() {
743 const int max_level = refinementLevels;
746 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
748 for (
int nn = 0; nn != numNodes; nn++)
749 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
751 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
755 Range tris(tri, tri);
758 tris, 1,
true, edges, moab::Interface::UNION);
763 Range nodes_at_front;
764 for (
int nn = 0; nn != numNodes; nn++) {
765 if (singular_nodes[nn]) {
767 CHKERR moab_ref.side_element(tri, 0, nn, ent);
768 nodes_at_front.insert(ent);
772 auto singular_edges = get_singular_edges();
775 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
776 for (
int ee = 0; ee != numEdges; ee++) {
777 if (singular_edges[ee]) {
779 CHKERR moab_ref.side_element(tri, 1, ee, ent);
780 CHKERR moab_ref.add_entities(meshset, &ent, 1);
786 for (
int ll = 0; ll != max_level; ll++) {
789 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
793 CHKERR moab_ref.get_adjacencies(
794 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
795 ref_edges = intersect(ref_edges, edges);
797 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
798 ref_edges = intersect(ref_edges, ents);
801 ->getEntitiesByTypeAndRefLevel(
803 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
807 ->updateMeshsetByEntitiesChildren(meshset,
809 meshset, MBEDGE,
true);
815 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
826 for (Range::iterator tit = tris.begin(); tit != tris.end();
830 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
831 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
834 auto &data = fe_ptr->dataOnElement[
H1];
835 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
836 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
838 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
840 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
841 double *tri_coords = &ref_coords(tt, 0);
844 auto det = t_normal.
l2();
845 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
846 for (
int dd = 0; dd != 2; dd++) {
847 ref_gauss_pts(dd, gg) =
848 shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
849 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
850 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
852 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
856 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
862 CHKERR refine_quadrature();
872 using ForcesAndSourcesCore::dataOnElement;
875 using ForcesAndSourcesCore::ForcesAndSourcesCore;
924 const char *list_rots[] = {
"small",
"moderate",
"large",
"no_h1"};
925 const char *list_symm[] = {
"symm",
"not_symm"};
926 const char *list_release[] = {
"griffith_force",
"griffith_skeleton"};
927 const char *list_stretches[] = {
"linear",
"log",
"log_quadratic"};
932 PetscInt choice_stretch = StretchSelector::LOG;
933 char analytical_expr_file_name[255] =
"analytical_expr.py";
935 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
937 CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
939 CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
941 CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
943 CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
945 CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
947 CHKERR PetscOptionsScalar(
"-viscosity_alpha_omega",
"rot viscosity",
"",
949 CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
951 CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
952 LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
954 CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
955 list_rots, NO_H1_CONFIGURATION + 1,
956 list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
957 CHKERR PetscOptionsEList(
"-symm",
"symmetric variant",
"", list_symm, 2,
958 list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
962 CHKERR PetscOptionsEList(
"-stretches",
"stretches",
"", list_stretches,
963 StretchSelector::STRETCH_SELECTOR_LAST,
964 list_stretches[choice_stretch], &choice_stretch,
967 CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
969 CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
971 CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
975 CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
979 CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
986 CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
988 CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
990 CHKERR PetscOptionsEList(
"-energy_release_variant",
"energy release variant",
991 "", list_release, 2, list_release[choice_release],
992 &choice_release, PETSC_NULLPTR);
993 CHKERR PetscOptionsInt(
"-nb_J_integral_levels",
"Number of J integarl levels",
997 char tag_name[255] =
"";
998 CHKERR PetscOptionsString(
"-internal_stress_tag_name",
999 "internal stress tag name",
"",
"", tag_name, 255,
1003 "-internal_stress_interp_order",
"internal stress interpolation order",
1005 CHKERR PetscOptionsBool(
"-internal_stress_voigt",
"Voigt index notation",
"",
1010 analytical_expr_file_name, 255, PETSC_NULLPTR);
1016 "Unsupported internal stress interpolation order %d",
1029 static_cast<EnergyReleaseSelector
>(choice_release);
1032 case StretchSelector::LINEAR:
1040 case StretchSelector::LOG:
1042 std::numeric_limits<float>::epsilon()) {
1058 case StretchSelector::LOG_QUADRATIC:
1064 "No logarithmic quadratic stretch for this case");
1080 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaU: -viscosity_alpha_u " <<
alphaU;
1081 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaW: -viscosity_alpha_w " <<
alphaW;
1083 <<
"alphaOmega: -viscosity_alpha_omega " <<
alphaOmega;
1087 MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
1095 MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
1097 <<
"Stretch: -stretches " << list_stretches[choice_stretch];
1123 <<
"Energy release variant: -energy_release_variant "
1126 <<
"Number of J integral levels: -nb_J_integral_levels "
1129#ifdef ENABLE_PYTHON_BINDING
1130 auto file_exists = [](std::string myfile) {
1131 std::ifstream file(myfile.c_str());
1138 if (file_exists(analytical_expr_file_name)) {
1139 MOFEM_LOG(
"EP", Sev::inform) << analytical_expr_file_name <<
" file found";
1143 analytical_expr_file_name);
1147 << analytical_expr_file_name <<
" file NOT found";
1160 auto get_tets = [&]() {
1166 auto get_tets_skin = [&]() {
1167 Range tets_skin_part;
1169 CHKERR skin.find_skin(0, get_tets(),
false, tets_skin_part);
1170 ParallelComm *pcomm =
1173 CHKERR pcomm->filter_pstatus(tets_skin_part,
1174 PSTATUS_SHARED | PSTATUS_MULTISHARED,
1175 PSTATUS_NOT, -1, &tets_skin);
1179 auto subtract_boundary_conditions = [&](
auto &&tets_skin) {
1185 tets_skin = subtract(tets_skin,
v.faces);
1189 tets_skin = subtract(tets_skin,
v.faces);
1194 tets_skin = subtract(tets_skin,
v.faces);
1200 auto add_blockset = [&](
auto block_name,
auto &&tets_skin) {
1203 tets_skin.merge(crack_faces);
1207 auto subtract_blockset = [&](
auto block_name,
auto &&tets_skin) {
1208 auto contact_range =
1210 tets_skin = subtract(tets_skin, contact_range);
1214 auto get_stress_trace_faces = [&](
auto &&tets_skin) {
1217 faces, moab::Interface::UNION);
1218 Range trace_faces = subtract(faces, tets_skin);
1222 auto tets = get_tets();
1226 auto trace_faces = get_stress_trace_faces(
1228 subtract_blockset(
"CONTACT",
1229 subtract_boundary_conditions(get_tets_skin()))
1236 boost::make_shared<Range>(subtract(trace_faces, *
contactFaces));
1256 auto add_broken_hdiv_field = [
this, meshset](
const std::string
field_name,
1262 auto get_side_map_hdiv = [&]() {
1265 std::pair<EntityType,
1280 get_side_map_hdiv(), MB_TAG_DENSE,
MF_ZERO);
1286 auto add_l2_field = [
this, meshset](
const std::string
field_name,
1287 const int order,
const int dim) {
1296 auto add_h1_field = [
this, meshset](
const std::string
field_name,
1297 const int order,
const int dim) {
1309 auto add_l2_field_by_range = [
this](
const std::string
field_name,
1310 const int order,
const int dim,
1311 const int field_dim,
Range &&r) {
1321 auto add_bubble_field = [
this, meshset](
const std::string
field_name,
1322 const int order,
const int dim) {
1328 auto field_order_table =
1329 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1330 auto get_cgg_bubble_order_zero = [](
int p) {
return 0; };
1331 auto get_cgg_bubble_order_tet = [](
int p) {
1334 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
1335 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
1336 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
1337 field_order_table[MBTET] = get_cgg_bubble_order_tet;
1344 auto add_user_l2_field = [
this, meshset](
const std::string
field_name,
1345 const int order,
const int dim) {
1351 auto field_order_table =
1352 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1353 auto zero_dofs = [](
int p) {
return 0; };
1355 field_order_table[MBVERTEX] = zero_dofs;
1356 field_order_table[MBEDGE] = zero_dofs;
1357 field_order_table[MBTRI] = zero_dofs;
1358 field_order_table[MBTET] = dof_l2_tet;
1376 auto get_hybridised_disp = [&]() {
1378 auto skin = subtract_boundary_conditions(get_tets_skin());
1380 faces.merge(intersect(bc.faces, skin));
1386 get_hybridised_disp());
1412 auto project_ho_geometry = [&](
auto field) {
1418 auto get_adj_front_edges = [&](
auto &front_edges) {
1419 Range front_crack_nodes;
1420 Range crack_front_edges_with_both_nodes_not_at_front;
1424 CHKERR moab.get_connectivity(front_edges, front_crack_nodes,
true);
1425 Range crack_front_edges;
1427 crack_front_edges, moab::Interface::UNION);
1428 crack_front_edges_with_both_nodes_not_at_front =
1429 subtract(crack_front_edges, front_edges);
1433 crack_front_edges_with_both_nodes_not_at_front =
send_type(
1434 mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
1436 return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
1437 boost::make_shared<Range>(
1438 crack_front_edges_with_both_nodes_not_at_front));
1445 auto [front_vertices, front_adj_edges] = get_adj_front_edges(*
frontEdges);
1456 (boost::format(
"crack_faces_%d.vtk") % rank).str(),
1459 (boost::format(
"front_edges_%d.vtk") % rank).str(),
1470 auto set_singular_dofs = [&](
auto &front_adj_edges,
auto &front_vertices) {
1478 MOFEM_LOG(
"EP", Sev::inform) <<
"Singularity eps " << beta;
1481 for (
auto edge : front_adj_edges) {
1484 CHKERR moab.get_connectivity(edge, conn, num_nodes,
false);
1486 CHKERR moab.get_coords(conn, num_nodes, coords);
1487 const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
1488 coords[5] - coords[2]};
1489 double dof[3] = {0, 0, 0};
1490 if (front_vertices.find(conn[0]) != front_vertices.end()) {
1491 for (
int dd = 0; dd != 3; dd++) {
1492 dof[dd] = -dir[dd] *
eps;
1494 }
else if (front_vertices.find(conn[1]) != front_vertices.end()) {
1495 for (
int dd = 0; dd != 3; dd++) {
1496 dof[dd] = +dir[dd] *
eps;
1501 const int idx = dit->get()->getEntDofIdx();
1503 dit->get()->getFieldData() = 0;
1505 dit->get()->getFieldData() = dof[idx];
1524 auto add_field_to_fe = [
this](
const std::string fe,
1589 auto set_fe_adjacency = [&](
auto fe_name) {
1592 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
1600 auto add_field_to_fe = [
this](
const std::string fe,
1612 Range natural_bc_elements;
1615 natural_bc_elements.merge(
v.faces);
1620 natural_bc_elements.merge(
v.faces);
1625 natural_bc_elements.merge(
v.faces);
1630 natural_bc_elements.merge(
v.faces);
1635 natural_bc_elements.merge(
v.faces);
1640 natural_bc_elements.merge(
v.faces);
1645 natural_bc_elements.merge(
v.faces);
1648 natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
1659 auto get_skin = [&](
auto &body_ents) {
1662 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
1667 Range boundary_ents;
1668 ParallelComm *pcomm =
1670 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
1671 PSTATUS_NOT, -1, &boundary_ents);
1672 return boundary_ents;
1806 auto remove_dofs_on_broken_skin = [&](
const std::string prb_name) {
1808 for (
int d : {0, 1, 2}) {
1809 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1811 ->getSideDofsOnBrokenSpaceEntities(
1822 CHKERR remove_dofs_on_broken_skin(
"ESHELBY_PLASTICITY");
1858 auto set_zero_block = [&]() {
1890 auto set_section = [&]() {
1892 PetscSection section;
1897 CHKERR PetscSectionDestroy(§ion);
1920BcDisp::BcDisp(std::string name, std::vector<double> attr,
Range faces)
1921 : blockName(name), faces(faces) {
1922 vals.resize(3,
false);
1923 flags.resize(3,
false);
1924 for (
int ii = 0; ii != 3; ++ii) {
1925 vals[ii] = attr[ii];
1926 flags[ii] =
static_cast<int>(attr[ii + 3]);
1929 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
1931 <<
"Add BCDisp vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1933 <<
"Add BCDisp flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1934 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " <<
faces.size();
1938 : blockName(name), faces(faces) {
1939 vals.resize(attr.size(),
false);
1940 for (
int ii = 0; ii != attr.size(); ++ii) {
1941 vals[ii] = attr[ii];
1947 : blockName(name), faces(faces) {
1948 vals.resize(3,
false);
1949 flags.resize(3,
false);
1950 for (
int ii = 0; ii != 3; ++ii) {
1951 vals[ii] = attr[ii];
1952 flags[ii] =
static_cast<int>(attr[ii + 3]);
1955 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
1957 <<
"Add BCForce vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1959 <<
"Add BCForce flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1960 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " <<
faces.size();
1964 std::vector<double> attr,
1966 : blockName(name), faces(faces) {
1969 if (attr.size() < 1) {
1971 "Wrong size of normal displacement BC");
1976 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc " << name;
1977 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc val " <<
val;
1979 <<
"Add NormalDisplacementBc nb. of faces " <<
faces.size();
1983 : blockName(name), faces(faces) {
1986 if (attr.size() < 1) {
1988 "Wrong size of normal displacement BC");
1993 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc " << name;
1994 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc val " <<
val;
1996 <<
"Add PressureBc nb. of faces " <<
faces.size();
2002 : blockName(name), ents(ents) {
2005 if (attr.size() < 2) {
2007 "Wrong size of external strain attribute");
2013 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain " << name;
2014 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain val " <<
val;
2015 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain bulk modulus K "
2018 <<
"Add ExternalStrain nb. of tets " <<
ents.size();
2024 std::vector<double> attr,
2026 : blockName(name), faces(faces) {
2029 if (attr.size() < 3) {
2031 "Wrong size of analytical displacement BC");
2034 flags.resize(3,
false);
2035 for (
int ii = 0; ii != 3; ++ii) {
2036 flags[ii] = attr[ii];
2039 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalDisplacementBc " << name;
2041 <<
"Add AnalyticalDisplacementBc flags " <<
flags[0] <<
" " <<
flags[1]
2044 <<
"Add AnalyticalDisplacementBc nb. of faces " <<
faces.size();
2048 std::vector<double> attr,
2050 : blockName(name), faces(faces) {
2053 if (attr.size() < 3) {
2055 "Wrong size of analytical traction BC");
2058 flags.resize(3,
false);
2059 for (
int ii = 0; ii != 3; ++ii) {
2060 flags[ii] = attr[ii];
2063 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalTractionBc " << name;
2065 <<
"Add AnalyticalTractionBc flags " <<
flags[0] <<
" " <<
flags[1]
2068 <<
"Add AnalyticalTractionBc nb. of faces " <<
faces.size();
2074 boost::shared_ptr<TractionFreeBc> &bc_ptr,
2075 const std::string contact_set_name) {
2080 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
2081 Range tets_skin_part;
2082 Skinner skin(&mField.get_moab());
2083 CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
2084 ParallelComm *pcomm =
2087 CHKERR pcomm->filter_pstatus(tets_skin_part,
2088 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2089 PSTATUS_NOT, -1, &tets_skin);
2092 for (
int dd = 0; dd != 3; ++dd)
2093 (*bc_ptr)[dd] = tets_skin;
2096 if (bcSpatialDispVecPtr)
2097 for (
auto &
v : *bcSpatialDispVecPtr) {
2099 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2101 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2103 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2107 if (bcSpatialRotationVecPtr)
2108 for (
auto &
v : *bcSpatialRotationVecPtr) {
2109 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2110 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2111 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2114 if (bcSpatialNormalDisplacementVecPtr)
2115 for (
auto &
v : *bcSpatialNormalDisplacementVecPtr) {
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 (bcSpatialAnalyticalDisplacementVecPtr)
2122 for (
auto &
v : *bcSpatialAnalyticalDisplacementVecPtr) {
2124 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2126 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2128 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2131 if (bcSpatialTractionVecPtr)
2132 for (
auto &
v : *bcSpatialTractionVecPtr) {
2133 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2134 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2135 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2138 if (bcSpatialAnalyticalTractionVecPtr)
2139 for (
auto &
v : *bcSpatialAnalyticalTractionVecPtr) {
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 (bcSpatialPressureVecPtr)
2146 for (
auto &
v : *bcSpatialPressureVecPtr) {
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);
2154 std::regex((boost::format(
"%s(.*)") % contact_set_name).str()))) {
2156 CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
2158 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
2159 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
2160 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
2177 return 2 * p_data + 1;
2183 return 2 * (p_data + 1);
2190 const int tag,
const bool do_rhs,
const bool do_lhs,
const bool calc_rates,
2192 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
2196 boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0,
MatrixDouble());
2197 fe->getUserPolynomialBase() =
2198 boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
2199 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2200 fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
2203 fe->getRuleHook = [](int, int, int) {
return -1; };
2204 fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges);
2210 dataAtPts->physicsPtr = physicalEquations;
2215 piolaStress, dataAtPts->getApproxPAtPts()));
2217 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
2219 piolaStress, dataAtPts->getDivPAtPts()));
2222 fe->getOpPtrVector().push_back(
2223 physicalEquations->returnOpCalculateStretchFromStress(
2224 dataAtPts, physicalEquations));
2226 fe->getOpPtrVector().push_back(
2228 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
2232 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
2234 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
2236 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
2240 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
2242 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
2247 spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
2250 fe->getOpPtrVector().push_back(
2252 stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
2253 fe->getOpPtrVector().push_back(
2255 stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
2259 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
2261 rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
2264 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2266 spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
2273 piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2276 bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2277 var_vec, MBMAXTYPE));
2279 rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
2281 piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
2283 spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
2286 fe->getOpPtrVector().push_back(
2287 physicalEquations->returnOpCalculateVarStretchFromStress(
2288 dataAtPts, physicalEquations));
2290 fe->getOpPtrVector().push_back(
2292 stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
2298 dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
2303 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
2304 tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
2311 const int tag,
const bool add_elastic,
const bool add_material,
2312 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
2313 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
2317 auto get_body_range = [
this](
auto name,
int dim) {
2318 std::map<int, Range> map;
2321 mField.getInterface<
MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2323 (boost::format(
"%s(.*)") % name).str()
2329 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2332 map[m_ptr->getMeshsetId()] = ents;
2338 auto rule_contact = [](int, int,
int o) {
return -1; };
2341 auto set_rule_contact = [refine](
2344 int order_col,
int order_data
2348 auto rule = 2 * order_data;
2354 fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2361 fe_rhs->getOpPtrVector().push_back(
2363 fe_rhs->getOpPtrVector().push_back(
2368 if (!internalStressTagName.empty()) {
2369 switch (internalStressInterpOrder) {
2371 fe_rhs->getOpPtrVector().push_back(
2375 fe_rhs->getOpPtrVector().push_back(
2380 "Unsupported internal stress interpolation order %d",
2381 internalStressInterpOrder);
2383 if (internalStressVoigt) {
2384 fe_rhs->getOpPtrVector().push_back(
2388 fe_rhs->getOpPtrVector().push_back(
2393 if (
auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
2394 stretchTensor, dataAtPts, externalStrainVecPtr, timeScaleMap)) {
2395 fe_rhs->getOpPtrVector().push_back(op);
2396 }
else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
2398 "OpSpatialPhysicalExternalStrain not implemented for this "
2402 fe_rhs->getOpPtrVector().push_back(
2403 physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
2406 fe_rhs->getOpPtrVector().push_back(
2408 fe_rhs->getOpPtrVector().push_back(
2410 fe_rhs->getOpPtrVector().push_back(
2413 auto set_hybridisation = [&](
auto &pip) {
2420 using SideEleOp = EleOnSide::UserDataOperator;
2421 using BdyEleOp = BoundaryEle::UserDataOperator;
2426 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2428 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2431 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2432 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2434 CHKERR EshelbianPlasticity::
2435 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2436 op_loop_skeleton_side->getOpPtrVector(), {L2},
2437 materialH1Positions, frontAdjEdges);
2441 auto broken_data_ptr =
2442 boost::make_shared<std::vector<BrokenBaseSideData>>();
2445 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2446 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2447 boost::make_shared<CGGUserPolynomialBase>();
2449 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2450 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2451 materialH1Positions, frontAdjEdges);
2452 op_loop_domain_side->getOpPtrVector().push_back(
2454 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
2455 op_loop_domain_side->getOpPtrVector().push_back(
2458 op_loop_domain_side->getOpPtrVector().push_back(
2462 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2464 GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
2466 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2467 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dHybrid(
2468 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
2469 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2470 op_loop_skeleton_side->getOpPtrVector().push_back(
2473 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dBroken(
2474 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
2477 pip.push_back(op_loop_skeleton_side);
2482 auto set_contact = [&](
auto &pip) {
2489 using SideEleOp = EleOnSide::UserDataOperator;
2490 using BdyEleOp = BoundaryEle::UserDataOperator;
2495 mField, contactElement,
SPACE_DIM - 1, Sev::noisy);
2497 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2498 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2499 CHKERR EshelbianPlasticity::
2500 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2501 op_loop_skeleton_side->getOpPtrVector(), {L2},
2502 materialH1Positions, frontAdjEdges);
2506 auto broken_data_ptr =
2507 boost::make_shared<std::vector<BrokenBaseSideData>>();
2510 auto contact_common_data_ptr =
2511 boost::make_shared<ContactOps::CommonData>();
2513 auto add_ops_domain_side = [&](
auto &pip) {
2517 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2518 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2519 boost::make_shared<CGGUserPolynomialBase>();
2521 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2522 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2523 materialH1Positions, frontAdjEdges);
2524 op_loop_domain_side->getOpPtrVector().push_back(
2527 op_loop_domain_side->getOpPtrVector().push_back(
2529 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2530 pip.push_back(op_loop_domain_side);
2534 auto add_ops_contact_rhs = [&](
auto &pip) {
2537 auto contact_sfd_map_range_ptr =
2538 boost::make_shared<std::map<int, Range>>(
2539 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2542 contactDisp, contact_common_data_ptr->contactDispPtr()));
2543 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2546 pip.push_back(
new OpTreeSearch(
2547 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2551 contactDisp, contact_common_data_ptr, contactTreeRhs,
2552 contact_sfd_map_range_ptr));
2555 broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
2561 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2562 CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
2565 pip.push_back(op_loop_skeleton_side);
2570 CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
2571 CHKERR set_contact(fe_rhs->getOpPtrVector());
2574 using BodyNaturalBC =
2576 Assembly<PETSC>::LinearForm<
GAUSS>;
2578 BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
2580 auto body_time_scale =
2581 boost::make_shared<DynamicRelaxationTimeScale>(
"body_force.txt");
2582 CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
2583 fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
2584 "BODY_FORCE", Sev::inform);
2588 fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2596 fe_lhs->getOpPtrVector().push_back(
2599 bubbleField, piolaStress, dataAtPts));
2600 fe_lhs->getOpPtrVector().push_back(
2604 fe_lhs->getOpPtrVector().push_back(
2605 physicalEquations->returnOpSpatialPhysical_du_du(
2606 stretchTensor, stretchTensor, dataAtPts, alphaU));
2608 stretchTensor, piolaStress, dataAtPts,
true));
2610 stretchTensor, bubbleField, dataAtPts,
true));
2612 stretchTensor, rotAxis, dataAtPts,
2613 symmetrySelector ==
SYMMETRIC ?
true :
false));
2617 spatialL2Disp, piolaStress, dataAtPts,
true));
2619 spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
2622 piolaStress, rotAxis, dataAtPts,
2623 symmetrySelector ==
SYMMETRIC ?
true :
false));
2625 bubbleField, rotAxis, dataAtPts,
2626 symmetrySelector ==
SYMMETRIC ?
true :
false));
2631 rotAxis, stretchTensor, dataAtPts,
false));
2634 rotAxis, piolaStress, dataAtPts,
false));
2636 rotAxis, bubbleField, dataAtPts,
false));
2639 rotAxis, rotAxis, dataAtPts, alphaOmega));
2641 auto set_hybridisation = [&](
auto &pip) {
2648 using SideEleOp = EleOnSide::UserDataOperator;
2649 using BdyEleOp = BoundaryEle::UserDataOperator;
2654 mField, skeletonElement,
SPACE_DIM - 1, Sev::noisy);
2656 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2659 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2660 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2661 CHKERR EshelbianPlasticity::
2662 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2663 op_loop_skeleton_side->getOpPtrVector(), {L2},
2664 materialH1Positions, frontAdjEdges);
2668 auto broken_data_ptr =
2669 boost::make_shared<std::vector<BrokenBaseSideData>>();
2672 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2673 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2674 boost::make_shared<CGGUserPolynomialBase>();
2676 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2677 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2678 materialH1Positions, frontAdjEdges);
2679 op_loop_domain_side->getOpPtrVector().push_back(
2682 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2684 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2685 op_loop_skeleton_side->getOpPtrVector().push_back(
2686 new OpC(hybridSpatialDisp, broken_data_ptr,
2687 boost::make_shared<double>(1.0),
true,
false));
2689 pip.push_back(op_loop_skeleton_side);
2694 auto set_contact = [&](
auto &pip) {
2701 using SideEleOp = EleOnSide::UserDataOperator;
2702 using BdyEleOp = BoundaryEle::UserDataOperator;
2707 mField, contactElement,
SPACE_DIM - 1, Sev::noisy);
2709 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2710 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2711 CHKERR EshelbianPlasticity::
2712 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2713 op_loop_skeleton_side->getOpPtrVector(), {L2},
2714 materialH1Positions, frontAdjEdges);
2718 auto broken_data_ptr =
2719 boost::make_shared<std::vector<BrokenBaseSideData>>();
2722 auto contact_common_data_ptr =
2723 boost::make_shared<ContactOps::CommonData>();
2725 auto add_ops_domain_side = [&](
auto &pip) {
2729 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2730 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2731 boost::make_shared<CGGUserPolynomialBase>();
2733 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2734 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2735 materialH1Positions, frontAdjEdges);
2736 op_loop_domain_side->getOpPtrVector().push_back(
2739 op_loop_domain_side->getOpPtrVector().push_back(
2741 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2742 pip.push_back(op_loop_domain_side);
2746 auto add_ops_contact_lhs = [&](
auto &pip) {
2749 contactDisp, contact_common_data_ptr->contactDispPtr()));
2750 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2753 pip.push_back(
new OpTreeSearch(
2754 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2759 auto contact_sfd_map_range_ptr =
2760 boost::make_shared<std::map<int, Range>>(
2761 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2764 contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
2765 contact_sfd_map_range_ptr));
2768 contactDisp, broken_data_ptr, contact_common_data_ptr,
2769 contactTreeRhs, contact_sfd_map_range_ptr));
2772 broken_data_ptr, contactDisp, contact_common_data_ptr,
2779 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2780 CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
2783 pip.push_back(op_loop_skeleton_side);
2788 CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
2789 CHKERR set_contact(fe_lhs->getOpPtrVector());
2799 const bool add_elastic,
const bool add_material,
2800 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
2801 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
2804 fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2805 fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2810 fe_rhs->getRuleHook = [](int, int, int) {
return -1; };
2811 fe_lhs->getRuleHook = [](int, int, int) {
return -1; };
2812 fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2813 fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2816 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2817 fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2819 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2820 fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2824 auto get_broken_op_side = [
this](
auto &pip) {
2827 using SideEleOp = EleOnSide::UserDataOperator;
2829 auto broken_data_ptr =
2830 boost::make_shared<std::vector<BrokenBaseSideData>>();
2833 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
2834 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2835 boost::make_shared<CGGUserPolynomialBase>();
2837 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2838 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2839 materialH1Positions, frontAdjEdges);
2840 op_loop_domain_side->getOpPtrVector().push_back(
2842 boost::shared_ptr<double> piola_scale_ptr(
new double);
2843 op_loop_domain_side->getOpPtrVector().push_back(
2844 physicalEquations->returnOpSetScale(piola_scale_ptr,
2845 physicalEquations));
2847 auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
2848 op_loop_domain_side->getOpPtrVector().push_back(
2850 piola_stress_mat_ptr));
2851 pip.push_back(op_loop_domain_side);
2852 return std::make_tuple(broken_data_ptr, piola_scale_ptr,
2853 piola_stress_mat_ptr);
2856 auto set_rhs = [&]() {
2859 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
2860 get_broken_op_side(fe_rhs->getOpPtrVector());
2862 fe_rhs->getOpPtrVector().push_back(
2863 new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr, timeScaleMap));
2865 broken_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
2868 broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
2870 fe_rhs->getOpPtrVector().push_back(
2872 piola_scale_ptr, timeScaleMap));
2873 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
2875 fe_rhs->getOpPtrVector().push_back(
2877 hybridSpatialDisp, hybrid_grad_ptr));
2879 hybridSpatialDisp, bcSpatialPressureVecPtr, piola_scale_ptr,
2880 hybrid_grad_ptr, timeScaleMap));
2882 hybridSpatialDisp, bcSpatialAnalyticalTractionVecPtr, piola_scale_ptr,
2885 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2886 fe_rhs->getOpPtrVector().push_back(
2890 hybridSpatialDisp, hybrid_ptr, piola_stress_mat_ptr,
2891 bcSpatialNormalDisplacementVecPtr, timeScaleMap));
2893 auto get_normal_disp_bc_faces = [&]() {
2896 return boost::make_shared<Range>(faces);
2901 using BdyEleOp = BoundaryEle::UserDataOperator;
2903 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2904 fe_rhs->getOpPtrVector().push_back(
new OpC_dBroken(
2905 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
2906 get_normal_disp_bc_faces()));
2911 auto set_lhs = [&]() {
2914 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
2915 get_broken_op_side(fe_lhs->getOpPtrVector());
2918 hybridSpatialDisp, bcSpatialNormalDisplacementVecPtr, timeScaleMap));
2920 hybridSpatialDisp, broken_data_ptr, bcSpatialNormalDisplacementVecPtr,
2923 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
2925 fe_rhs->getOpPtrVector().push_back(
2927 hybridSpatialDisp, hybrid_grad_ptr));
2929 hybridSpatialDisp, bcSpatialPressureVecPtr, hybrid_grad_ptr,
2932 auto get_normal_disp_bc_faces = [&]() {
2935 return boost::make_shared<Range>(faces);
2940 using BdyEleOp = BoundaryEle::UserDataOperator;
2942 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2943 fe_lhs->getOpPtrVector().push_back(
new OpC(
2944 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
2945 true,
true, get_normal_disp_bc_faces()));
2959 boost::shared_ptr<ContactTree> &fe_contact_tree
2965 auto get_body_range = [
this](
auto name,
int dim,
auto sev) {
2966 std::map<int, Range> map;
2969 mField.getInterface<
MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2971 (boost::format(
"%s(.*)") % name).str()
2977 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2980 map[m_ptr->getMeshsetId()] = ents;
2981 MOFEM_LOG(
"EPSYNC", sev) <<
"Meshset: " << m_ptr->getMeshsetId() <<
" "
2982 << ents.size() <<
" entities";
2989 auto get_map_skin = [
this](
auto &&map) {
2990 ParallelComm *pcomm =
2993 Skinner skin(&mField.get_moab());
2994 for (
auto &
m : map) {
2996 CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
2998 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2999 PSTATUS_NOT, -1,
nullptr),
3001 m.second.swap(skin_faces);
3011 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3013 auto calcs_side_traction = [&](
auto &pip) {
3017 using SideEleOp = EleOnSide::UserDataOperator;
3019 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
3020 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3021 boost::make_shared<CGGUserPolynomialBase>();
3022 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3023 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3024 materialH1Positions, frontAdjEdges);
3025 op_loop_domain_side->getOpPtrVector().push_back(
3027 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3028 boost::make_shared<double>(1.0)));
3029 pip.push_back(op_loop_domain_side);
3033 auto add_contact_three = [&]() {
3035 auto tree_moab_ptr = boost::make_shared<moab::Core>();
3036 fe_contact_tree = boost::make_shared<ContactTree>(
3037 mField, tree_moab_ptr, spaceOrder,
3038 get_body_range(
"CONTACT",
SPACE_DIM - 1, Sev::inform));
3039 fe_contact_tree->getOpPtrVector().push_back(
3041 contactDisp, contact_common_data_ptr->contactDispPtr()));
3042 CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
3043 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3044 fe_contact_tree->getOpPtrVector().push_back(
3046 fe_contact_tree->getOpPtrVector().push_back(
3047 new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
3051 CHKERR add_contact_three();
3061 CHKERR setContactElementRhsOps(contactTreeRhs);
3063 CHKERR setVolumeElementOps(tag,
true,
false, elasticFeRhs, elasticFeLhs);
3064 CHKERR setFaceElementOps(
true,
false, elasticBcRhs, elasticBcLhs);
3067 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
3069 auto get_op_contact_bc = [&]() {
3072 mField, contactElement,
SPACE_DIM - 1, Sev::noisy, adj_cache);
3073 return op_loop_side;
3081 boost::shared_ptr<FEMethod> null;
3083 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3116 bool set_ts_monitor) {
3118#ifdef ENABLE_PYTHON_BINDING
3119 auto setup_sdf = [&]() {
3120 boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
3122 auto file_exists = [](std::string myfile) {
3123 std::ifstream file(myfile.c_str());
3130 char sdf_file_name[255] =
"sdf.py";
3132 sdf_file_name, 255, PETSC_NULLPTR);
3134 if (file_exists(sdf_file_name)) {
3135 MOFEM_LOG(
"EP", Sev::inform) << sdf_file_name <<
" file found";
3136 sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
3137 CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
3138 ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
3139 MOFEM_LOG(
"EP", Sev::inform) <<
"SdfPython initialized";
3141 MOFEM_LOG(
"EP", Sev::warning) << sdf_file_name <<
" file NOT found";
3143 return sdf_python_ptr;
3145 auto sdf_python_ptr = setup_sdf();
3148 auto setup_ts_monitor = [&]() {
3149 boost::shared_ptr<TsCtx>
ts_ctx;
3151 if (set_ts_monitor) {
3153 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
3157 MOFEM_LOG(
"EP", Sev::inform) <<
"TS monitor setup";
3158 return std::make_tuple(
ts_ctx);
3161 auto setup_snes_monitor = [&]() {
3164 CHKERR TSGetSNES(ts, &snes);
3166 CHKERR SNESMonitorSet(snes,
3169 (
void *)(snes_ctx.get()), PETSC_NULLPTR);
3170 MOFEM_LOG(
"EP", Sev::inform) <<
"SNES monitor setup";
3174 auto setup_snes_conergence_test = [&]() {
3177 auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
3178 PetscReal snorm, PetscReal fnorm,
3179 SNESConvergedReason *reason,
void *cctx) {
3182 CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
3186 CHKERR SNESGetSolutionUpdate(snes, &x_update);
3187 CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
3255 auto setup_section = [&]() {
3256 PetscSection section_raw;
3259 CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
3260 for (
int ff = 0; ff != num_fields; ff++) {
3268 auto set_vector_on_mesh = [&]() {
3272 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
3273 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
3274 MOFEM_LOG(
"EP", Sev::inform) <<
"Vector set on mesh";
3278 auto setup_schur_block_solver = [&]() {
3279 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver";
3280 CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
3281 CHKERR TSSetFromOptions(ts);
3284 boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
3285 if constexpr (
A == AssemblyType::BLOCK_MAT) {
3290 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver done";
3297#ifdef ENABLE_PYTHON_BINDING
3298 return std::make_tuple(setup_sdf(), setup_ts_monitor(),
3299 setup_snes_monitor(), setup_snes_conergence_test(),
3300 setup_section(), set_vector_on_mesh(),
3301 setup_schur_block_solver());
3303 return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
3304 setup_snes_conergence_test(), setup_section(),
3305 set_vector_on_mesh(), setup_schur_block_solver());
3313 PetscBool debug_model = PETSC_FALSE;
3317 if (debug_model == PETSC_TRUE) {
3319 auto post_proc = [&](TS ts, PetscReal
t, Vec u, Vec u_t, Vec u_tt, Vec
F,
3324 CHKERR TSGetSNES(ts, &snes);
3326 CHKERR SNESGetIterationNumber(snes, &it);
3327 std::string file_name =
"snes_iteration_" + std::to_string(it) +
".h5m";
3328 CHKERR postProcessResults(1, file_name,
F, u_t);
3329 std::string file_skel_name =
3330 "snes_iteration_skel_" + std::to_string(it) +
".h5m";
3332 auto get_material_force_tag = [&]() {
3333 auto &moab = mField.get_moab();
3340 CHKERR calculateFaceMaterialForce(1, ts);
3341 CHKERR postProcessSkeletonResults(1, file_skel_name,
F,
3342 {get_material_force_tag()});
3346 ts_ctx_ptr->tsDebugHook = post_proc;
3351 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3353 CHKERR VecDuplicate(x, &xx);
3354 CHKERR VecZeroEntries(xx);
3355 CHKERR TS2SetSolution(ts, x, xx);
3358 CHKERR TSSetSolution(ts, x);
3361 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3362 {elasticFeLhs.get(), elasticFeRhs.get()});
3367 CHKERR TSSolve(ts, PETSC_NULLPTR);
3369 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3370 {elasticFeLhs.get(), elasticFeRhs.get()});
3373 CHKERR TSGetSNES(ts, &snes);
3374 int lin_solver_iterations;
3375 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
3377 <<
"Number of linear solver iterations " << lin_solver_iterations;
3379 PetscBool test_cook_flg = PETSC_FALSE;
3382 if (test_cook_flg) {
3383 constexpr int expected_lin_solver_iterations = 11;
3384 if (lin_solver_iterations > expected_lin_solver_iterations)
3387 "Expected number of iterations is different than expected %d > %d",
3388 lin_solver_iterations, expected_lin_solver_iterations);
3391 PetscBool test_sslv116_flag = PETSC_FALSE;
3393 &test_sslv116_flag, PETSC_NULLPTR);
3395 if (test_sslv116_flag) {
3396 double max_val = 0.0;
3397 double min_val = 0.0;
3398 auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3400 auto ent_type = ent_ptr->getEntType();
3401 if (ent_type == MBVERTEX) {
3402 max_val = std::max(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], max_val);
3403 min_val = std::min(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], min_val);
3408 field_min_max, spatialH1Disp);
3410 double global_max_val = 0.0;
3411 double global_min_val = 0.0;
3412 MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
3414 MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
3417 <<
"Max " << spatialH1Disp <<
" value: " << global_max_val;
3419 <<
"Min " << spatialH1Disp <<
" value: " << global_min_val;
3421 double ref_max_val = 0.00767;
3422 double ref_min_val = -0.00329;
3423 if (std::abs(global_max_val - ref_max_val) > 1e-5) {
3425 "Incorrect max value of the displacement field: %f != %f",
3426 global_max_val, ref_max_val);
3428 if (std::abs(global_min_val - ref_min_val) > 4e-5) {
3430 "Incorrect min value of the displacement field: %f != %f",
3431 global_min_val, ref_min_val);
3445 double final_time = 1;
3446 double delta_time = 0.1;
3448 PetscBool ts_h1_update = PETSC_FALSE;
3450 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Dynamic Relaxation Options",
3453 CHKERR PetscOptionsScalar(
"-dynamic_final_time",
3454 "dynamic relaxation final time",
"", final_time,
3455 &final_time, PETSC_NULLPTR);
3456 CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
3457 "dynamic relaxation final time",
"", delta_time,
3458 &delta_time, PETSC_NULLPTR);
3459 CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
3460 max_it, &max_it, PETSC_NULLPTR);
3461 CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
3462 ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
3466 auto setup_ts_monitor = [&]() {
3467 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*
this);
3470 auto monitor_ptr = setup_ts_monitor();
3472 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3473 {elasticFeLhs.get(), elasticFeRhs.get()});
3477 double ts_delta_time;
3478 CHKERR TSGetTimeStep(ts, &ts_delta_time);
3490 monitor_ptr->ts_u = PETSC_NULLPTR;
3491 monitor_ptr->ts_t = dynamicTime;
3492 monitor_ptr->ts_step = dynamicStep;
3494 MOFEM_LOG(
"EP", Sev::inform) <<
"IT " << dynamicStep <<
" Time "
3495 << dynamicTime <<
" delta time " << delta_time;
3496 dynamicTime += delta_time;
3499 for (; dynamicTime < final_time; dynamicTime += delta_time) {
3500 MOFEM_LOG(
"EP", Sev::inform) <<
"IT " << dynamicStep <<
" Time "
3501 << dynamicTime <<
" delta time " << delta_time;
3503 CHKERR TSSetStepNumber(ts, 0);
3505 CHKERR TSSetTimeStep(ts, ts_delta_time);
3506 if (!ts_h1_update) {
3509 CHKERR TSSolve(ts, PETSC_NULLPTR);
3510 if (!ts_h1_update) {
3514 monitor_ptr->ts_u = PETSC_NULLPTR;
3515 monitor_ptr->ts_t = dynamicTime;
3516 monitor_ptr->ts_step = dynamicStep;
3520 if (dynamicStep > max_it)
3525 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3526 {elasticFeLhs.get(), elasticFeRhs.get()});
3534 auto set_block = [&](
auto name,
int dim) {
3535 std::map<int, Range> map;
3536 auto set_tag_impl = [&](
auto name) {
3541 std::regex((boost::format(
"%s(.*)") % name).str())
3544 for (
auto bc : bcs) {
3546 CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
3548 map[bc->getMeshsetId()] = r;
3553 CHKERR set_tag_impl(name);
3555 return std::make_pair(name, map);
3558 auto set_skin = [&](
auto &&map) {
3559 for (
auto &
m : map.second) {
3566 auto set_tag = [&](
auto &&map) {
3568 auto name = map.first;
3569 int def_val[] = {-1};
3571 mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER,
th,
3572 MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
3574 for (
auto &
m : map.second) {
3582 listTagsToTransfer.push_back(set_tag(set_skin(set_block(
"BODY", 3))));
3583 listTagsToTransfer.push_back(set_tag(set_skin(set_block(
"MAT_ELASTIC", 3))));
3584 listTagsToTransfer.push_back(
3585 set_tag(set_skin(set_block(
"MAT_NEOHOOKEAN", 3))));
3586 listTagsToTransfer.push_back(set_tag(set_block(
"CONTACT", 2)));
3593 Vec f_residual, Vec var_vector,
3594 std::vector<Tag> tags_to_transfer) {
3600 rval = mField.get_moab().tag_get_handle(
"CRACK",
th);
3601 if (
rval == MB_SUCCESS) {
3602 CHKERR mField.get_moab().tag_delete(
th);
3604 int def_val[] = {0};
3605 CHKERR mField.get_moab().tag_get_handle(
3606 "CRACK", 1, MB_TYPE_INTEGER,
th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
3608 CHKERR mField.get_moab().tag_clear_data(
th, *crackFaces, mark);
3609 tags_to_transfer.push_back(
th);
3611 auto get_tag = [&](
auto name,
auto dim) {
3612 auto &mob = mField.get_moab();
3614 double def_val[] = {0., 0., 0.};
3615 CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3616 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3621 tags_to_transfer.push_back(get_tag(
"MaterialForce", 3));
3626 for (
auto t : listTagsToTransfer) {
3627 tags_to_transfer.push_back(
t);
3635 struct exclude_sdf {
3636 exclude_sdf(
Range &&r) : map(r) {}
3637 bool operator()(
FEMethod *fe_method_ptr) {
3639 if (map.find(ent) != map.end()) {
3649 contactTreeRhs->exeTestHook =
3654 auto get_post_proc = [&](
auto &post_proc_mesh,
auto sense) {
3656 auto post_proc_ptr =
3657 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3658 mField, post_proc_mesh);
3659 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3660 post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
3663 auto domain_ops = [&](
auto &fe,
int sense) {
3665 fe.getUserPolynomialBase() =
3667 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3668 fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
3670 auto piola_scale_ptr = boost::make_shared<double>(1.0);
3671 fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
3672 piola_scale_ptr, physicalEquations));
3674 piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
3676 bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
3679 fe.getOpPtrVector().push_back(
3680 physicalEquations->returnOpCalculateStretchFromStress(
3681 dataAtPts, physicalEquations));
3683 fe.getOpPtrVector().push_back(
3685 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
3690 piolaStress, dataAtPts->getVarPiolaPts(),
3691 boost::make_shared<double>(1), vec));
3693 bubbleField, dataAtPts->getVarPiolaPts(),
3694 boost::make_shared<double>(1), vec, MBMAXTYPE));
3696 rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
3698 fe.getOpPtrVector().push_back(
3699 physicalEquations->returnOpCalculateVarStretchFromStress(
3700 dataAtPts, physicalEquations));
3702 fe.getOpPtrVector().push_back(
3704 stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
3709 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
3711 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
3714 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
3716 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
3718 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
3720 fe.getOpPtrVector().push_back(
3724 fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
3725 tag,
true,
false, dataAtPts, physicalEquations));
3727 physicalEquations->returnOpCalculateEnergy(dataAtPts,
nullptr)) {
3728 fe.getOpPtrVector().push_back(op);
3737 struct OpSidePPMap :
public OpPPMap {
3738 OpSidePPMap(moab::Interface &post_proc_mesh,
3739 std::vector<EntityHandle> &map_gauss_pts,
3740 DataMapVec data_map_scalar, DataMapMat data_map_vec,
3741 DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
3743 :
OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
3744 data_map_vec, data_map_mat, data_symm_map_mat),
3751 if (tagSense != OpPPMap::getSkeletonSense())
3763 vec_fields[
"SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
3764 vec_fields[
"SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
3765 vec_fields[
"Omega"] = dataAtPts->getRotAxisAtPts();
3766 vec_fields[
"ContactDisplacement"] = dataAtPts->getContactL2AtPts();
3767 vec_fields[
"AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
3768 vec_fields[
"X"] = dataAtPts->getLargeXH1AtPts();
3770 vec_fields[
"EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
3774 vec_fields[
"VarOmega"] = dataAtPts->getVarRotAxisPts();
3775 vec_fields[
"VarSpatialDisplacementL2"] =
3776 boost::make_shared<MatrixDouble>();
3778 spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], vec, MBTET));
3782 vec_fields[
"ResSpatialDisplacementL2"] =
3783 boost::make_shared<MatrixDouble>();
3785 spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], vec, MBTET));
3786 vec_fields[
"ResOmega"] = boost::make_shared<MatrixDouble>();
3788 rotAxis, vec_fields[
"ResOmega"], vec, MBTET));
3792 mat_fields[
"PiolaStress"] = dataAtPts->getApproxPAtPts();
3794 mat_fields[
"VarPiolaStress"] = dataAtPts->getVarPiolaPts();
3798 mat_fields[
"ResPiolaStress"] = boost::make_shared<MatrixDouble>();
3800 piolaStress, mat_fields[
"ResPiolaStress"],
3801 boost::make_shared<double>(1), vec));
3803 bubbleField, mat_fields[
"ResPiolaStress"],
3804 boost::make_shared<double>(1), vec, MBMAXTYPE));
3806 if (!internalStressTagName.empty()) {
3807 mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
3808 switch (internalStressInterpOrder) {
3810 fe.getOpPtrVector().push_back(
3814 fe.getOpPtrVector().push_back(
3819 "Unsupported internal stress interpolation order %d",
3820 internalStressInterpOrder);
3825 mat_fields_symm[
"LogSpatialStretch"] =
3826 dataAtPts->getLogStretchTensorAtPts();
3827 mat_fields_symm[
"SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
3829 mat_fields_symm[
"VarLogSpatialStretch"] =
3830 dataAtPts->getVarLogStreachPts();
3835 mat_fields_symm[
"ResLogSpatialStretch"] =
3836 boost::make_shared<MatrixDouble>();
3837 fe.getOpPtrVector().push_back(
3839 stretchTensor, mat_fields_symm[
"ResLogSpatialStretch"], vec,
3844 fe.getOpPtrVector().push_back(
3848 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3865 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3871 post_proc_ptr->getOpPtrVector().push_back(
3873 dataAtPts->getContactL2AtPts()));
3875 auto X_h1_ptr = boost::make_shared<MatrixDouble>();
3877 post_proc_ptr->getOpPtrVector().push_back(
3879 dataAtPts->getLargeXH1AtPts()));
3884 domain_ops(*(op_loop_side->getSideFEPtr()), sense);
3885 post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
3887 return post_proc_ptr;
3891 auto calcs_side_traction_and_displacements = [&](
auto &post_proc_ptr,
3894 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3898 using SideEleOp = EleOnSide::UserDataOperator;
3900 mField, elementVolumeName,
SPACE_DIM, Sev::noisy);
3901 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3903 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3904 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3905 materialH1Positions, frontAdjEdges);
3906 op_loop_domain_side->getOpPtrVector().push_back(
3908 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3909 boost::make_shared<double>(1.0)));
3910 pip.push_back(op_loop_domain_side);
3912 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3915 contactDisp, contact_common_data_ptr->contactDispPtr()));
3916 pip.push_back(
new OpTreeSearch(
3917 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
3919 &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
3926 pip.push_back(op_this);
3927 auto contact_residual = boost::make_shared<MatrixDouble>();
3928 op_this->getOpPtrVector().push_back(
3930 contactDisp, contact_residual,
3933 op_this->getOpPtrVector().push_back(
3937 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3941 {{
"res_contact", contact_residual}},
3955 auto post_proc_mesh = boost::make_shared<moab::Core>();
3956 auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
3957 auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
3958 CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
3959 post_proc_ptr->getOpPtrVector());
3965 CHKERR mField.get_moab().get_adjacencies(own_tets,
SPACE_DIM - 1,
true,
3966 own_faces, moab::Interface::UNION);
3968 auto get_post_negative = [&](
auto &&ents) {
3969 auto crack_faces_pos = ents;
3970 auto crack_faces_neg = crack_faces_pos;
3971 auto skin =
get_skin(mField, own_tets);
3972 auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
3973 for (
auto f : crack_on_proc_skin) {
3975 CHKERR mField.get_moab().get_adjacencies(&f, 1,
SPACE_DIM,
false, tet);
3976 tet = intersect(tet, own_tets);
3977 int side_number, sense, offset;
3978 CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
3981 crack_faces_neg.erase(f);
3983 crack_faces_pos.erase(f);
3986 return std::make_pair(crack_faces_pos, crack_faces_neg);
3989 auto get_crack_faces = [&](
auto crack_faces) {
3990 auto get_adj = [&](
auto e,
auto dim) {
3992 CHKERR mField.get_moab().get_adjacencies(e, dim,
true, adj,
3993 moab::Interface::UNION);
3996 auto tets = get_adj(crack_faces, 3);
3997 auto faces = subtract(get_adj(tets, 2), crack_faces);
3998 tets = subtract(tets, get_adj(faces, 3));
3999 return subtract(crack_faces, get_adj(tets, 2));
4002 auto [crack_faces_pos, crack_faces_neg] =
4003 get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
4005 auto only_crack_faces = [&](
FEMethod *fe_method_ptr) {
4006 auto ent = fe_method_ptr->getFEEntityHandle();
4007 if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
4013 auto only_negative_crack_faces = [&](
FEMethod *fe_method_ptr) {
4014 auto ent = fe_method_ptr->getFEEntityHandle();
4015 if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
4021 auto get_adj_front = [&]() {
4022 auto skeleton_faces = *skeletonFaces;
4024 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2,
true, adj_front,
4025 moab::Interface::UNION);
4027 adj_front = intersect(adj_front, skeleton_faces);
4028 adj_front = subtract(adj_front, *crackFaces);
4029 adj_front = intersect(own_faces, adj_front);
4033 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4034 post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
4036 auto post_proc_begin =
4040 post_proc_ptr->exeTestHook = only_crack_faces;
4041 post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
4043 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4045 post_proc_negative_sense_ptr, 0,
4046 mField.get_comm_size());
4048 constexpr bool debug =
false;
4051 auto [adj_front_pos, adj_front_neg] =
4054 auto only_front_faces_pos = [adj_front_pos](
FEMethod *fe_method_ptr) {
4055 auto ent = fe_method_ptr->getFEEntityHandle();
4056 if (adj_front_pos.find(ent) == adj_front_pos.end()) {
4062 auto only_front_faces_neg = [adj_front_neg](
FEMethod *fe_method_ptr) {
4063 auto ent = fe_method_ptr->getFEEntityHandle();
4064 if (adj_front_neg.find(ent) == adj_front_neg.end()) {
4070 post_proc_ptr->exeTestHook = only_front_faces_pos;
4072 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4073 post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
4075 post_proc_negative_sense_ptr, 0,
4076 mField.get_comm_size());
4081 CHKERR post_proc_end.writeFile(file.c_str());
4088 std::vector<Tag> tags_to_transfer) {
4093 auto post_proc_mesh = boost::make_shared<moab::Core>();
4094 auto post_proc_ptr =
4095 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
4097 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
4101 auto hybrid_disp = boost::make_shared<MatrixDouble>();
4102 post_proc_ptr->getOpPtrVector().push_back(
4104 post_proc_ptr->getOpPtrVector().push_back(
4108 auto op_loop_domain_side =
4111 post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4114 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4115 boost::make_shared<CGGUserPolynomialBase>();
4116 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4117 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4119 op_loop_domain_side->getOpPtrVector().push_back(
4122 op_loop_domain_side->getOpPtrVector().push_back(
4125 op_loop_domain_side->getOpPtrVector().push_back(
4129 op_loop_domain_side->getOpPtrVector().push_back(
4133 op_loop_domain_side->getOpPtrVector().push_back(
4141 vec_fields[
"HybridDisplacement"] = hybrid_disp;
4142 vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
4144 mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
4145 mat_fields[
"HybridDisplacementGradient"] =
4148 mat_fields_symm[
"LogSpatialStretch"] =
dataAtPts->getLogStretchTensorAtPts();
4150 post_proc_ptr->getOpPtrVector().push_back(
4154 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4169 auto hybrid_res = boost::make_shared<MatrixDouble>();
4170 post_proc_ptr->getOpPtrVector().push_back(
4175 post_proc_ptr->getOpPtrVector().push_back(
4179 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4183 {{
"res_hybrid", hybrid_res}},
4194 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4196 auto post_proc_begin =
4203 CHKERR post_proc_end.writeFile(file.c_str());
4211 constexpr bool debug =
false;
4213 auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
4214 std::vector<Tag> tags;
4215 tags.reserve(names.size());
4216 auto create_and_clean = [&]() {
4218 for (
auto n : names) {
4219 tags.push_back(
Tag());
4220 auto &tag = tags.back();
4222 auto rval = moab.tag_get_handle(
n.first.c_str(), tag);
4223 if (
rval == MB_SUCCESS) {
4224 moab.tag_delete(tag);
4226 double def_val[] = {0., 0., 0.};
4227 CHKERR moab.tag_get_handle(
n.first.c_str(),
n.second, MB_TYPE_DOUBLE,
4228 tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4236 enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
4238 auto tags = get_tags_vec({{
"MaterialForce", 3},
4240 {
"GriffithForce", 1},
4241 {
"FacePressure", 1}});
4243 auto calculate_material_forces = [&]() {
4249 auto get_face_material_force_fe = [&]() {
4251 auto fe_ptr = boost::make_shared<FaceEle>(
mField);
4252 fe_ptr->getRuleHook = [](int, int, int) {
return -1; };
4253 fe_ptr->setRuleHook =
4257 EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
4261 fe_ptr->getOpPtrVector().push_back(
4264 auto op_loop_domain_side =
4267 fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4271 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4272 boost::make_shared<CGGUserPolynomialBase>();
4274 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4275 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4277 op_loop_domain_side->getOpPtrVector().push_back(
4280 op_loop_domain_side->getOpPtrVector().push_back(
4292 op_loop_domain_side->getOpPtrVector().push_back(
4296 op_loop_domain_side->getOpPtrVector().push_back(
4300 op_loop_domain_side->getOpPtrVector().push_back(
4305 op_loop_domain_side->getOpPtrVector().push_back(
4311 auto integrate_face_material_force_fe = [&](
auto &&face_energy_fe) {
4319 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
4322 CHKERR VecGetLocalSize(
v.second, &size);
4324 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
4325 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( "
4326 << low <<
" " << high <<
" ) ";
4330 CHKERR print_loc_size(face_exchange,
"material face_exchange",
4334 mField.
get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
4341 "front_skin_faces_material_force_" +
4350 CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
4355 auto calculate_front_material_force = [&]() {
4359 auto get_conn = [&](
auto e) {
4362 "get connectivity");
4366 auto get_conn_range = [&](
auto e) {
4369 "get connectivity");
4373 auto get_adj = [&](
auto e,
auto dim) {
4380 auto get_adj_range = [&](
auto e,
auto dim) {
4383 moab::Interface::UNION),
4388 auto get_material_force = [&](
auto r,
auto th) {
4393 return material_forces;
4398 auto crack_edges = get_adj_range(*
crackFaces, 1);
4399 auto front_nodes = get_conn_range(*
frontEdges);
4406 Range all_skin_faces;
4407 Range all_front_faces;
4410 auto calculate_edge_direction = [&](
auto e) {
4415 "get connectivity");
4416 std::array<double, 6> coords;
4421 &coords[0], &coords[1], &coords[2]};
4423 &coords[3], &coords[4], &coords[5]};
4426 t_dir(
i) = t_p1(
i) - t_p0(
i);
4431 auto calculate_force_through_node = [&]() {
4443 auto body_skin_conn = get_conn_range(body_skin);
4446 for (
auto n : front_nodes) {
4447 auto adj_tets = get_adj(
n, 3);
4449 auto conn = get_conn_range(adj_tets);
4450 adj_tets = get_adj_range(conn, 3);
4453 auto material_forces = get_material_force(skin_faces, tags[0]);
4457 all_skin_faces.merge(skin_faces);
4461 auto calculate_node_material_force = [&]() {
4463 getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
4465 for (
auto face : skin_faces) {
4468 t_face_force_tmp(
I) = t_face_T(
I);
4471 auto face_tets = intersect(get_adj(face, 3), adj_tets);
4473 if (face_tets.empty()) {
4477 if (face_tets.size() != 1) {
4479 "face_tets.size() != 1");
4482 int side_number, sense, offset;
4486 "moab side number");
4487 t_face_force_tmp(
I) *= sense;
4488 t_node_force(
I) += t_face_force_tmp(
I);
4491 return t_node_force;
4494 auto calculate_crack_area_growth_direction =
4495 [&](
auto n,
auto &t_node_force) {
4498 auto boundary_node = intersect(
Range(
n,
n), body_skin_conn);
4499 if (boundary_node.size()) {
4500 auto faces = intersect(get_adj(
n, 2), body_skin);
4501 for (
auto f : faces) {
4504 f, &t_normal_face(0));
4505 t_project(
I) += t_normal_face(
I);
4507 t_project.normalize();
4514 if (boundary_node.size()) {
4515 t_Q(
I,
J) -= t_project(
I) * t_project(
J);
4518 auto adj_faces = intersect(get_adj(
n, 2), *
crackFaces);
4519 if (adj_faces.empty()) {
4521 intersect(get_adj(
n, 1), unite(*
frontEdges, body_edges));
4523 for (
auto e : adj_edges) {
4524 auto t_dir = calculate_edge_direction(e);
4530 t_node_force_tmp(
I) = t_node_force(
I);
4532 t_area_dir(
I) = -t_node_force_tmp(
I);
4533 t_area_dir(
I) *=
l / 2;
4538 auto front_edges = get_adj(
n, 1);
4540 for (
auto f : adj_faces) {
4545 std::array<double, 9> coords;
4551 coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
4552 auto n_it = std::find(conn, conn + num_nodes,
n);
4553 auto n_index = std::distance(conn, n_it);
4556 t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
4557 t_d_normal(0, n_index, 2),
4559 t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
4560 t_d_normal(1, n_index, 2),
4562 t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
4563 t_d_normal(2, n_index, 2)};
4566 t_projected_hessian(
I,
J) =
4567 t_Q(
I,
K) * (t_face_hessian(
K,
L) * t_Q(
L,
J));
4570 t_face_normal(
I) * t_projected_hessian(
I,
K) / 2.;
4576 auto t_node_force = calculate_node_material_force();
4580 &
n, 1, &t_node_force(0)),
4583 auto get_area_dir = [&]() {
4584 auto adj_faces = get_adj_range(adj_tets, 2);
4585 auto t_area_dir = calculate_crack_area_growth_direction(
4588 auto adj_edges = intersect(get_adj_range(adj_tets, 1),
4590 auto seed_n = get_conn_range(adj_edges);
4592 skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
4593 seed_n = subtract(seed_n, skin_adj_edges);
4595 for (
auto sn : seed_n) {
4596 auto t_area_dir_sn = calculate_crack_area_growth_direction(
4598 t_area_dir(
I) += t_area_dir_sn(
I);
4600 for (
auto sn : skin_adj_edges) {
4601 auto t_area_dir_sn = calculate_crack_area_growth_direction(
4603 t_area_dir(
I) += t_area_dir_sn(
I) / 2;
4610 auto t_area_dir = get_area_dir();
4616 auto griffith = -t_node_force(
I) * t_area_dir(
I) /
4617 (t_area_dir(
K) * t_area_dir(
K));
4625 auto ave_node_force = [&](
auto th) {
4630 auto conn = get_conn(e);
4631 auto data = get_material_force(conn,
th);
4632 auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
4634 for (
auto n : conn) {
4635 t_edge(
I) += t_node(
I);
4638 t_edge(
I) /= conn.size();
4641 calculate_edge_direction(e);
4659 auto ave_node_griffith_energy = [&](
auto th) {
4664 tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
4667 &e, 1, &t_edge_area_dir(0));
4668 double griffith_energy = -t_edge_force(
I) * t_edge_area_dir(
I) /
4669 (t_edge_area_dir(
K) * t_edge_area_dir(
K));
4675 CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
4676 CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
4677 CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
4682 CHKERR calculate_force_through_node();
4686 auto adj_faces = get_adj(e, 2);
4687 auto crack_face = intersect(get_adj(e, 2), *
crackFaces);
4691 all_front_faces.merge(adj_faces);
4697 &e, 1, &t_edge_force(0));
4699 calculate_edge_direction(e);
4708 for (
auto f : adj_faces) {
4712 int side_number, sense, offset;
4715 auto dot = -sense * t_cross(
I) * t_normal(
I);
4717 tags[ExhangeTags::GRIFFITHFORCE], &
f, 1, &dot),
4725 CHKERR TSGetStepNumber(ts, &ts_step);
4727 "front_edges_material_force_" +
4728 std::to_string(ts_step) +
".vtk",
4731 "front_skin_faces_material_force_" +
4732 std::to_string(ts_step) +
".vtk",
4735 "front_faces_material_force_" +
4736 std::to_string(ts_step) +
".vtk",
4745 mField.
get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
4747 mField.
get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
4754 auto print_results = [&]() {
4757 auto get_conn_range = [&](
auto e) {
4760 "get connectivity");
4764 auto get_tag_data = [&](
auto &ents,
auto tag,
auto dim) {
4765 std::vector<double> data(ents.size() * dim);
4772 auto at_nodes = [&]() {
4775 auto material_force =
4776 get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
4777 auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
4778 auto griffith_force =
4779 get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
4780 std::vector<double> coords(conn.size() * 3);
4784 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Material force at nodes";
4785 for (
size_t i = 0;
i < conn.size(); ++
i) {
4787 <<
"Node " << conn[
i] <<
" coords "
4788 << coords[
i * 3 + 0] <<
" " << coords[
i * 3 + 1] <<
" "
4789 << coords[
i * 3 + 2] <<
" material force "
4790 << material_force[
i * 3 + 0] <<
" "
4791 << material_force[
i * 3 + 1] <<
" "
4792 << material_force[
i * 3 + 2] <<
" area growth "
4793 << area_growth[
i * 3 + 0] <<
" " << area_growth[
i * 3 + 1]
4794 <<
" " << area_growth[
i * 3 + 2] <<
" griffith force "
4795 << griffith_force[
i];
4805 CHKERR calculate_material_forces();
4806 CHKERR calculate_front_material_force();
4813 bool set_orientation) {
4816 constexpr bool debug =
false;
4817 constexpr auto sev = Sev::verbose;
4822 Range body_skin_edges;
4824 moab::Interface::UNION);
4825 Range boundary_skin_verts;
4827 boundary_skin_verts,
true);
4830 Range geometry_edges_verts;
4832 geometry_edges_verts,
true);
4833 Range crack_faces_verts;
4836 Range crack_faces_edges;
4838 *
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
4839 Range crack_faces_tets;
4841 *
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
4847 moab::Interface::UNION);
4848 Range front_verts_edges;
4850 front_verts, 1,
true, front_verts_edges, moab::Interface::UNION);
4852 auto get_tags_vec = [&](
auto tag_name,
int dim) {
4853 std::vector<Tag> tags(1);
4858 auto create_and_clean = [&]() {
4861 auto rval = moab.tag_get_handle(tag_name, tags[0]);
4862 if (
rval == MB_SUCCESS) {
4863 moab.tag_delete(tags[0]);
4865 double def_val[] = {0., 0., 0.};
4866 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
4867 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4876 auto get_adj_front = [&](
bool subtract_crack) {
4879 adj_front, moab::Interface::UNION);
4881 adj_front = subtract(adj_front, *
crackFaces);
4887 auto th_front_position = get_tags_vec(
"FrontPosition", 3);
4888 auto th_max_face_energy = get_tags_vec(
"MaxFaceEnergy", 1);
4892 auto get_crack_adj_tets = [&](
auto r) {
4893 Range crack_faces_conn;
4895 Range crack_faces_conn_tets;
4897 true, crack_faces_conn_tets,
4898 moab::Interface::UNION);
4899 return crack_faces_conn_tets;
4902 auto get_layers_for_sides = [&](
auto &side) {
4903 std::vector<Range> layers;
4907 auto get_adj = [&](
auto &r,
int dim) {
4910 moab::Interface::UNION);
4914 auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
4919 Range front_faces = get_adj(front_nodes, 2);
4920 front_faces = subtract(front_faces, *
crackFaces);
4921 auto front_tets = get_tets(front_nodes);
4922 auto front_side = intersect(side, front_tets);
4923 layers.push_back(front_side);
4926 adj_faces = intersect(adj_faces, front_faces);
4927 auto adj_faces_tets = get_tets(adj_faces);
4928 adj_faces_tets = intersect(adj_faces_tets, front_tets);
4929 layers.push_back(unite(layers.back(), adj_faces_tets));
4930 if (layers.back().size() == layers[layers.size() - 2].size()) {
4941 auto layers_top = get_layers_for_sides(sides_pair.first);
4942 auto layers_bottom = get_layers_for_sides(sides_pair.second);
4954 MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
4956 for (
auto &r : layers_top) {
4957 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
4960 "layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
4965 for (
auto &r : layers_bottom) {
4966 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
4969 "layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
4975 auto get_cross = [&](
auto t_dir,
auto f) {
4987 auto get_sense = [&](
auto f,
auto e) {
4988 int side, sense, offset;
4991 return std::make_tuple(side, sense, offset);
4994 auto calculate_edge_direction = [&](
auto e,
auto normalize =
true) {
4998 std::array<double, 6> coords;
5001 &coords[0], &coords[1], &coords[2]};
5003 &coords[3], &coords[4], &coords[5]};
5006 t_dir(
i) = t_p1(
i) - t_p0(
i);
5012 auto evaluate_face_energy_and_set_orientation = [&](
auto front_edges,
5019 Tag th_material_force;
5021 case GRIFFITH_FORCE:
5022 case GRIFFITH_SKELETON:
5029 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5030 "Unknown energy release selector");
5038 auto find_maximal_face_energy = [&](
auto front_edges,
auto front_faces,
5039 auto &edge_face_max_energy_map) {
5048 for (
auto e : front_edges) {
5050 double griffith_force;
5056 faces = subtract(intersect(faces, front_faces), body_skin);
5057 std::vector<double> face_energy(faces.size());
5059 face_energy.data());
5060 auto max_energy_it =
5061 std::max_element(face_energy.begin(), face_energy.end());
5063 max_energy_it != face_energy.end() ? *max_energy_it : 0;
5065 edge_face_max_energy_map[e] =
5066 std::make_tuple(faces[max_energy_it - face_energy.begin()],
5067 griffith_force,
static_cast<double>(0));
5069 <<
"Edge " << e <<
" griffith force " << griffith_force
5070 <<
" max face energy " << max_energy <<
" factor "
5071 << max_energy / griffith_force;
5073 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
5095 auto calculate_face_orientation = [&](
auto &edge_face_max_energy_map) {
5098 auto up_down_face = [&](
5100 auto &face_angle_map_up,
5101 auto &face_angle_map_down
5106 for (
auto &
m : edge_face_max_energy_map) {
5108 auto [max_face, energy, opt_angle] =
m.second;
5112 faces = intersect(faces, front_faces);
5116 moab::Interface::UNION);
5117 if (adj_tets.size()) {
5122 moab::Interface::UNION);
5123 if (adj_tets.size()) {
5125 Range adj_tets_faces;
5128 adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
5129 moab::Interface::UNION);
5130 adj_tets_faces = intersect(adj_tets_faces, faces);
5135 get_cross(calculate_edge_direction(e,
true), max_face);
5136 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
5137 t_cross_max(
i) *= sense_max;
5139 for (
auto t : adj_tets) {
5140 Range adj_tets_faces;
5142 &
t, 1,
SPACE_DIM - 1,
false, adj_tets_faces);
5143 adj_tets_faces = intersect(adj_tets_faces, faces);
5145 subtract(adj_tets_faces,
Range(max_face, max_face));
5147 if (adj_tets_faces.size() == 1) {
5151 auto t_cross = get_cross(calculate_edge_direction(e,
true),
5153 auto [side, sense, offset] =
5154 get_sense(adj_tets_faces[0], e);
5155 t_cross(
i) *= sense;
5156 double dot = t_cross(
i) * t_cross_max(
i);
5157 auto angle = std::acos(dot);
5161 th_face_energy, adj_tets_faces, &face_energy);
5163 auto [side_face, sense_face, offset_face] =
5164 get_sense(
t, max_face);
5166 if (sense_face > 0) {
5167 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
5171 face_angle_map_down[e] = std::make_tuple(
5172 face_energy, -angle, adj_tets_faces[0]);
5183 auto calc_optimal_angle = [&](
5185 auto &face_angle_map_up,
5186 auto &face_angle_map_down
5191 for (
auto &
m : edge_face_max_energy_map) {
5193 auto &[max_face, e0,
a0] =
m.second;
5195 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
5197 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
5198 face_angle_map_down.find(e) == face_angle_map_down.end()) {
5203 case GRIFFITH_FORCE:
5204 case GRIFFITH_SKELETON: {
5206 Tag th_material_force;
5211 th_material_force, &e, 1, &t_material_force(0));
5212 auto material_force_magnitude = t_material_force.
l2();
5213 if (material_force_magnitude <
5214 std::numeric_limits<double>::epsilon()) {
5219 auto t_edge_dir = calculate_edge_direction(e,
true);
5220 auto t_cross_max = get_cross(t_edge_dir, max_face);
5221 auto [side, sense, offset] = get_sense(max_face, e);
5222 t_cross_max(sense) *= sense;
5229 t_cross_max.normalize();
5232 t_material_force(
J) * t_cross_max(
K);
5233 a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
5236 <<
"Optimal angle " <<
a0 <<
" energy " << e0;
5242 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5243 "Unknown energy release selector");
5253 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5255 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5256 face_angle_map_down;
5257 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
5258 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
5262 auto th_angle = get_tags_vec(
"Angle", 1);
5264 for (
auto &
m : face_angle_map_up) {
5265 auto [e,
a, face] =
m.second;
5270 for (
auto &
m : face_angle_map_down) {
5271 auto [e,
a, face] =
m.second;
5276 Range max_energy_faces;
5277 for (
auto &
m : edge_face_max_energy_map) {
5278 auto [face, e, angle] =
m.second;
5279 max_energy_faces.insert(face);
5295 auto get_conn = [&](
auto e) {
5302 auto get_adj = [&](
auto e,
auto dim) {
5305 e, dim,
false, adj, moab::Interface::UNION),
5310 auto get_coords = [&](
auto v) {
5318 auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
5321 auto t_edge_dir = calculate_edge_direction(e,
true);
5322 auto [side, sense, offset] = get_sense(
f, e);
5323 t_edge_dir(
i) *= sense;
5324 t_edge_dir.normalize();
5325 t_edge_dir(
i) *= angle;
5330 t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
5331 return std::make_tuple(t_normal, t_rotated_normal);
5334 auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
5335 auto &t_move,
auto gamma) {
5336 auto index = adj_vertex_tets_verts.index(
v);
5338 for (
auto ii : {0, 1, 2}) {
5339 coords[3 * index + ii] += gamma * t_move(ii);
5346 auto tets_quality = [&](
auto quality,
auto &adj_vertex_tets_verts,
5347 auto &adj_vertex_tets,
auto &coords) {
5348 for (
auto t : adj_vertex_tets) {
5352 std::array<double, 12> tet_coords;
5353 for (
auto n = 0;
n != 4; ++
n) {
5354 auto index = adj_vertex_tets_verts.index(conn[
n]);
5358 for (
auto ii = 0; ii != 3; ++ii) {
5359 tet_coords[3 *
n + ii] = coords[3 * index + ii];
5363 if (!std::isnormal(q))
5365 quality = std::min(quality, q);
5371 auto calculate_free_face_node_displacement =
5372 [&](
auto &edge_face_max_energy_map) {
5374 auto get_vertex_edges = [&](
auto vertex) {
5381 vertex_edges = subtract(vertex_edges, front_verts_edges);
5383 if (boundary_skin_verts.size() &&
5384 boundary_skin_verts.find(vertex[0]) !=
5385 boundary_skin_verts.end()) {
5386 MOFEM_LOG(
"EP", sev) <<
"Boundary vertex";
5387 vertex_edges = intersect(vertex_edges, body_skin_edges);
5389 if (geometry_edges_verts.size() &&
5390 geometry_edges_verts.find(vertex[0]) !=
5391 geometry_edges_verts.end()) {
5392 MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
5393 vertex_edges = intersect(vertex_edges, geometry_edges);
5395 if (crack_faces_verts.size() &&
5396 crack_faces_verts.find(vertex[0]) !=
5397 crack_faces_verts.end()) {
5398 MOFEM_LOG(
"EP", sev) <<
"Crack face vertex";
5399 vertex_edges = intersect(vertex_edges, crack_faces_edges);
5406 return vertex_edges;
5411 using Bundle = std::vector<
5417 std::map<EntityHandle, Bundle> edge_bundle_map;
5419 for (
auto &
m : edge_face_max_energy_map) {
5421 auto edge =
m.first;
5422 auto &[max_face, energy, opt_angle] =
m.second;
5425 auto [t_normal, t_rotated_normal] =
5426 get_rotated_normal(edge, max_face, opt_angle);
5428 auto front_vertex = get_conn(
Range(
m.first,
m.first));
5429 auto adj_tets = get_adj(
Range(max_face, max_face), 3);
5430 auto adj_tets_faces = get_adj(adj_tets, 2);
5431 auto adj_front_faces = subtract(
5432 intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
5434 if (adj_front_faces.size() > 3)
5436 "adj_front_faces.size()>3");
5440 &t_material_force(0));
5441 std::vector<double> griffith_energy(adj_front_faces.size());
5443 th_face_energy, adj_front_faces, griffith_energy.data());
5446 auto set_edge_bundle = [&](
auto min_gamma) {
5447 for (
auto rotated_f : adj_front_faces) {
5449 double rotated_face_energy =
5450 griffith_energy[adj_front_faces.index(rotated_f)];
5452 auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
5454 if (vertex.size() != 1) {
5456 "Wrong number of vertex to move");
5458 auto front_vertex_edges_vertex = get_conn(
5459 intersect(get_adj(front_vertex, 1), crack_faces_edges));
5461 vertex, front_vertex_edges_vertex);
5462 if (vertex.empty()) {
5466 auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
5469 subtract(body_skin_edges, crack_faces_edges));
5472 for (;
c < 10; ++
c) {
5474 subtract(get_adj(faces, 1), seen_front_edges);
5475 if (front_edges.size() == 0) {
5478 auto front_connected_edges =
5479 intersect(front_edges, whole_front);
5480 if (front_connected_edges.size()) {
5481 seen_front_edges.merge(front_connected_edges);
5484 faces.merge(get_adj(front_edges, 2));
5491 double rotated_face_cardinality = face_cardinality(
5497 rotated_face_cardinality = std::max(rotated_face_cardinality,
5500 auto t_vertex_coords = get_coords(vertex);
5501 auto vertex_edges = get_vertex_edges(vertex);
5510 for (
auto e_used_to_move_detection : vertex_edges) {
5511 auto edge_conn = get_conn(
Range(e_used_to_move_detection,
5512 e_used_to_move_detection));
5513 edge_conn = subtract(edge_conn, vertex);
5523 t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
5527 (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5529 (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5532 constexpr double eps =
5533 std::numeric_limits<double>::epsilon();
5534 if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
5537 t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
5539 auto check_rotated_face_directoon = [&]() {
5541 t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
5544 (t_material_force(
i) / t_material_force.
l2()) *
5546 return -dot > 0 ? true :
false;
5549 if (check_rotated_face_directoon()) {
5552 <<
"Crack edge " << edge <<
" moved face "
5554 <<
" edge: " << e_used_to_move_detection
5555 <<
" face direction/energy " << rotated_face_energy
5556 <<
" face cardinality " << rotated_face_cardinality
5557 <<
" gamma: " << gamma;
5559 auto &bundle = edge_bundle_map[edge];
5560 bundle.emplace_back(rotated_f, e_used_to_move_detection,
5561 vertex[0], t_move, 1,
5562 rotated_face_cardinality, gamma);
5569 set_edge_bundle(std::numeric_limits<double>::epsilon());
5570 if (edge_bundle_map[edge].empty()) {
5571 set_edge_bundle(-1.);
5575 return edge_bundle_map;
5578 auto get_sort_by_energy = [&](
auto &edge_face_max_energy_map) {
5579 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
5582 for (
auto &
m : edge_face_max_energy_map) {
5584 auto &[max_face, energy, opt_angle] =
m.second;
5585 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
5588 return sort_by_energy;
5591 auto set_tag = [&](
auto &&adj_edges_map,
auto &&sort_by_energy) {
5594 Tag th_face_pressure;
5596 mField.
get_moab().tag_get_handle(
"FacePressure", th_face_pressure),
5598 auto get_face_pressure = [&](
auto face) {
5607 <<
"Number of edges to check " << sort_by_energy.size();
5609 enum face_energy { POSITIVE, NEGATIVE };
5610 constexpr bool skip_negative =
true;
5612 for (
auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
5616 for (
auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
5619 auto energy = it->first;
5620 auto [max_edge, max_face, opt_angle] = it->second;
5622 auto face_pressure = get_face_pressure(max_face);
5623 if (skip_negative) {
5624 if (fe == face_energy::POSITIVE) {
5625 if (face_pressure < 0) {
5627 <<
"Skip negative face " << max_face <<
" with energy "
5628 << energy <<
" and pressure " << face_pressure;
5635 <<
"Check face " << max_face <<
" edge " << max_edge
5636 <<
" energy " << energy <<
" optimal angle " << opt_angle
5637 <<
" face pressure " << face_pressure;
5639 auto jt = adj_edges_map.find(max_edge);
5640 if (jt == adj_edges_map.end()) {
5642 <<
"Edge " << max_edge <<
" not found in adj_edges_map";
5645 auto &bundle = jt->second;
5647 auto find_max_in_bundle_impl = [&](
auto edge,
auto &bundle,
5654 double max_quality = -2;
5655 double max_quality_evaluated = -2;
5656 double min_cardinality = std::numeric_limits<double>::max();
5660 for (
auto &b : bundle) {
5661 auto &[face, move_edge, vertex, t_move, quality, cardinality,
5664 auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
5665 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
5666 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
5668 adj_vertex_tets_verts, coords.data()),
5671 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
5672 quality = tets_quality(quality, adj_vertex_tets_verts,
5673 adj_vertex_tets, coords);
5675 auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
5679 return ((edge_gamma < 0) ? (q / 2) : q) / pow(
c, 2);
5683 if (eval_quality(quality, cardinality, edge_gamma) >=
5684 max_quality_evaluated) {
5685 max_quality = quality;
5686 min_cardinality = cardinality;
5687 vertex_max = vertex;
5689 move_edge_max = move_edge;
5690 t_move_last(
i) = t_move(
i);
5691 max_quality_evaluated =
5692 eval_quality(max_quality, min_cardinality, edge_gamma);
5696 return std::make_tuple(vertex_max, face_max, t_move_last,
5697 max_quality, min_cardinality);
5700 auto find_max_in_bundle = [&](
auto edge,
auto &bundle) {
5701 auto b_org_bundle = bundle;
5702 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
5703 auto &[vertex_max, face_max, t_move_last, max_quality,
5705 if (max_quality < 0) {
5706 for (
double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
5707 bundle = b_org_bundle;
5708 r = find_max_in_bundle_impl(edge, bundle, gamma);
5709 auto &[vertex_max, face_max, t_move_last, max_quality,
5712 <<
"Back tracking: gamma " << gamma <<
" edge " << edge
5713 <<
" quality " << max_quality <<
" cardinality "
5715 if (max_quality > 0.01) {
5717 t_move_last(
I) *= gamma;
5728 auto set_tag_to_vertex_and_face = [&](
auto &&r,
auto &quality) {
5730 auto &[
v,
f, t_move, q, cardinality] = r;
5732 if ((q > 0 && std::isnormal(q)) && energy > 0) {
5735 <<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" "
5736 << max_edge <<
" move " << t_move <<
" energy " << energy
5737 <<
" quality " << q <<
" cardinality " << cardinality;
5748 double quality = -2;
5749 CHKERR set_tag_to_vertex_and_face(
5751 find_max_in_bundle(max_edge, bundle),
5757 if (quality > 0 && std::isnormal(quality) && energy > 0) {
5759 <<
"Crack face set with quality: " << quality;
5772 MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
5773 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
5774 edge_face_max_energy_map;
5775 CHKERR find_maximal_face_energy(front_edges, front_faces,
5776 edge_face_max_energy_map);
5777 CHKERR calculate_face_orientation(edge_face_max_energy_map);
5779 MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
5782 calculate_free_face_node_displacement(edge_face_max_energy_map),
5783 get_sort_by_energy(edge_face_max_energy_map)
5791 CHKERR evaluate_face_energy_and_set_orientation(
5792 *
frontEdges, get_adj_front(
false), sides_pair, th_front_position);
5810 auto get_max_moved_faces = [&]() {
5811 Range max_moved_faces;
5812 auto adj_front = get_adj_front(
false);
5813 std::vector<double> face_energy(adj_front.size());
5815 face_energy.data());
5816 for (
int i = 0;
i != adj_front.size(); ++
i) {
5817 if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
5818 max_moved_faces.insert(adj_front[
i]);
5822 return boost::make_shared<Range>(max_moved_faces);
5833 "max_moved_faces_" +
5848 Tag th_front_position;
5850 mField.
get_moab().tag_get_handle(
"FrontPosition", th_front_position);
5855 std::vector<double> coords(3 * verts.size());
5857 std::vector<double> pos(3 * verts.size());
5859 for (
int i = 0;
i != 3 * verts.size(); ++
i) {
5860 coords[
i] += pos[
i];
5863 double zero[] = {0., 0., 0.};
5868 constexpr bool debug =
false;
5873 "set_coords_faces_" +
5884 constexpr bool potential_crack_debug =
false;
5885 if constexpr (potential_crack_debug) {
5888 Range crack_front_verts;
5893 Range crack_front_faces;
5895 true, crack_front_faces,
5896 moab::Interface::UNION);
5897 crack_front_faces = intersect(crack_front_faces, add_ents);
5904 auto get_crack_faces = [&]() {
5912 auto get_extended_crack_faces = [&]() {
5913 auto get_faces_of_crack_front_verts = [&](
auto crack_faces_org) {
5914 ParallelComm *pcomm =
5919 if (!pcomm->rank()) {
5921 auto get_nodes = [&](
auto &&e) {
5924 "get connectivity");
5928 auto get_adj = [&](
auto &&e,
auto dim,
5929 auto t = moab::Interface::UNION) {
5941 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
5944 auto front_block_nodes = get_nodes(front_block_edges);
5948 s = crack_faces.size();
5950 auto crack_face_nodes = get_nodes(crack_faces_org);
5951 auto crack_faces_edges =
5952 get_adj(crack_faces_org, 1, moab::Interface::UNION);
5955 front_block_edges = subtract(front_block_edges, crack_skin);
5956 auto crack_skin_nodes = get_nodes(crack_skin);
5957 crack_skin_nodes.merge(front_block_nodes);
5959 auto crack_skin_faces =
5960 get_adj(crack_skin, 2, moab::Interface::UNION);
5962 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
5964 crack_faces = crack_faces_org;
5965 for (
auto f : crack_skin_faces) {
5966 auto edges = intersect(
5967 get_adj(
Range(
f,
f), 1, moab::Interface::UNION), crack_skin);
5971 if (edges.size() == 2) {
5973 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
5977 if (edges.size() == 2) {
5978 auto edge_conn = get_nodes(
Range(edges));
5979 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
5981 if (faces.size() == 2) {
5982 auto edge0_conn = get_nodes(
Range(edges[0], edges[0]));
5983 auto edge1_conn = get_nodes(
Range(edges[1], edges[1]));
5984 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
5990 moab::Interface::INTERSECT),
5995 if (node_edges.size()) {
6000 auto get_t_dir = [&](
auto e_conn) {
6001 auto other_node = subtract(e_conn,
edges_conn);
6005 t_dir(
i) -= t_v0(
i);
6011 get_t_dir(edge0_conn)(
i) + get_t_dir(edge1_conn)(
i);
6014 t_crack_surface_ave_dir(
i) = 0;
6015 for (
auto e : node_edges) {
6016 auto e_conn = get_nodes(
Range(e, e));
6017 auto t_dir = get_t_dir(e_conn);
6018 t_crack_surface_ave_dir(
i) += t_dir(
i);
6021 auto dot = t_ave_dir(
i) * t_crack_surface_ave_dir(
i);
6024 if (dot < -std::numeric_limits<double>::epsilon()) {
6025 crack_faces.insert(
f);
6028 crack_faces.insert(
f);
6032 }
else if (edges.size() == 3) {
6033 crack_faces.insert(
f);
6037 if (edges.size() == 1) {
6039 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6042 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6043 front_block_edges));
6044 if (edges.size() == 2) {
6045 crack_faces.insert(
f);
6051 crack_faces_org = crack_faces;
6053 }
while (s != crack_faces.size());
6059 return get_faces_of_crack_front_verts(get_crack_faces());
6064 get_extended_crack_faces());
6067 auto reconstruct_crack_faces = [&](
auto crack_faces) {
6068 ParallelComm *pcomm =
6074 Range new_crack_faces;
6075 if (!pcomm->rank()) {
6077 auto get_nodes = [&](
auto &&e) {
6080 "get connectivity");
6084 auto get_adj = [&](
auto &&e,
auto dim,
6085 auto t = moab::Interface::UNION) {
6093 auto get_test_on_crack_surface = [&]() {
6094 auto crack_faces_nodes =
6095 get_nodes(crack_faces);
6096 auto crack_faces_tets =
6097 get_adj(crack_faces_nodes, 3,
6098 moab::Interface::UNION);
6102 auto crack_faces_tets_nodes =
6103 get_nodes(crack_faces_tets);
6104 crack_faces_tets_nodes =
6105 subtract(crack_faces_tets_nodes, crack_faces_nodes);
6107 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
6108 moab::Interface::UNION));
6110 get_adj(crack_faces_tets, 2,
6111 moab::Interface::UNION);
6113 new_crack_faces.merge(crack_faces);
6115 return std::make_tuple(new_crack_faces, crack_faces_tets);
6118 auto carck_faces_test_edges = [&](
auto faces,
auto tets) {
6119 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
6120 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
6121 moab::Interface::UNION);
6122 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
6125 adj_faces_edges.merge(geometry_edges);
6126 adj_faces_edges.merge(front_block_edges);
6128 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
6129 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
6130 auto boundary_test_nodes_edges =
6131 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
6132 auto boundary_test_nodes_edges_nodes = subtract(
6133 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
6135 boundary_tets_edges =
6136 subtract(boundary_test_nodes_edges,
6137 get_adj(boundary_test_nodes_edges_nodes, 1,
6138 moab::Interface::UNION));
6145 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6146 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
6148 body_skin = intersect(body_skin, adj_tets_faces);
6149 body_skin_edges = subtract(
6150 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
6153 for (
auto e : body_skin_edges) {
6154 auto adj_tet = intersect(
6155 get_adj(
Range(e, e), 3, moab::Interface::INTERSECT), tets);
6156 if (adj_tet.size() == 1) {
6157 boundary_tets_edges.insert(e);
6161 return boundary_tets_edges;
6164 auto p = get_test_on_crack_surface();
6165 auto &[new_crack_faces, crack_faces_tets] = p;
6176 auto boundary_tets_edges =
6177 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
6179 boundary_tets_edges);
6181 auto resolve_surface = [&](
auto boundary_tets_edges,
6182 auto crack_faces_tets) {
6183 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
6184 auto crack_faces_tets_faces =
6185 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6187 Range all_removed_faces;
6188 Range all_removed_tets;
6192 while (size != crack_faces_tets.size()) {
6194 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6198 auto skin_skin_nodes = get_nodes(skin_skin);
6200 size = crack_faces_tets.size();
6202 <<
"Crack faces tets size " << crack_faces_tets.size()
6203 <<
" crack faces size " << crack_faces_tets_faces.size();
6204 auto skin_tets_nodes = subtract(
6205 get_nodes(skin_tets),
6206 boundary_tets_edges_nodes);
6208 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
6210 Range removed_nodes;
6211 Range tets_to_remove;
6212 Range faces_to_remove;
6213 for (
auto n : skin_tets_nodes) {
6215 intersect(get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT),
6217 if (tets.size() == 0) {
6221 auto hole_detetction = [&]() {
6223 get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT);
6228 if (adj_tets.size() == 0) {
6229 return std::make_pair(
6231 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6236 std::vector<Range> tets_groups;
6237 auto test_adj_tets = adj_tets;
6238 while (test_adj_tets.size()) {
6240 Range seed =
Range(test_adj_tets[0], test_adj_tets[0]);
6241 while (seed.size() != seed_size) {
6243 subtract(get_adj(seed, 2, moab::Interface::UNION),
6246 seed_size = seed.size();
6248 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
6251 tets_groups.push_back(seed);
6252 test_adj_tets = subtract(test_adj_tets, seed);
6254 if (tets_groups.size() == 1) {
6256 return std::make_pair(
6258 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6264 Range tets_to_remove;
6265 Range faces_to_remove;
6266 for (
auto &r : tets_groups) {
6267 auto f = get_adj(r, 2, moab::Interface::UNION);
6268 auto t = intersect(get_adj(
f, 3, moab::Interface::UNION),
6271 if (
f.size() > faces_to_remove.size() ||
6272 faces_to_remove.size() == 0) {
6273 faces_to_remove =
f;
6278 <<
"Hole detection: faces to remove "
6279 << faces_to_remove.size() <<
" tets to remove "
6280 << tets_to_remove.size();
6281 return std::make_pair(faces_to_remove, tets_to_remove);
6284 if (tets.size() < tets_to_remove.size() ||
6285 tets_to_remove.size() == 0) {
6287 auto [h_faces_to_remove, h_tets_to_remove] =
6289 faces_to_remove = h_faces_to_remove;
6290 tets_to_remove = h_tets_to_remove;
6298 all_removed_faces.merge(faces_to_remove);
6299 all_removed_tets.merge(tets_to_remove);
6302 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
6303 crack_faces_tets_faces =
6304 subtract(crack_faces_tets_faces, faces_to_remove);
6309 boost::lexical_cast<std::string>(counter) +
".vtk",
6312 "faces_to_remove_" +
6313 boost::lexical_cast<std::string>(counter) +
".vtk",
6317 boost::lexical_cast<std::string>(counter) +
".vtk",
6320 "crack_faces_tets_faces_" +
6321 boost::lexical_cast<std::string>(counter) +
".vtk",
6322 crack_faces_tets_faces);
6324 "crack_faces_tets_" +
6325 boost::lexical_cast<std::string>(counter) +
".vtk",
6331 auto cese_internal_faces = [&]() {
6334 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
6336 subtract(adj_faces, skin_tets);
6337 auto adj_tets = get_adj(adj_faces, 3,
6338 moab::Interface::UNION);
6341 subtract(crack_faces_tets,
6344 crack_faces_tets_faces =
6345 subtract(crack_faces_tets_faces, adj_faces);
6347 all_removed_faces.merge(adj_faces);
6348 all_removed_tets.merge(adj_tets);
6352 <<
"Remove internal faces size " << adj_faces.size()
6353 <<
" tets size " << adj_tets.size();
6357 auto case_only_one_free_edge = [&]() {
6360 for (
auto t :
Range(crack_faces_tets)) {
6362 auto adj_faces = get_adj(
6364 moab::Interface::UNION);
6365 auto crack_surface_edges =
6366 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6369 moab::Interface::UNION);
6372 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6373 crack_surface_edges);
6374 adj_edges = subtract(
6376 boundary_tets_edges);
6378 if (adj_edges.size() == 1) {
6380 subtract(crack_faces_tets,
6384 auto faces_to_remove =
6385 get_adj(adj_edges, 2, moab::Interface::UNION);
6388 crack_faces_tets_faces =
6389 subtract(crack_faces_tets_faces, faces_to_remove);
6391 all_removed_faces.merge(faces_to_remove);
6392 all_removed_tets.merge(
Range(
t,
t));
6394 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Remove free one edges ";
6398 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6399 crack_faces_tets_faces =
6400 subtract(crack_faces_tets_faces, all_removed_faces);
6405 auto cese_flat_tet = [&](
auto max_adj_edges) {
6412 auto body_skin_edges =
6413 get_adj(body_skin, 1, moab::Interface::UNION);
6415 for (
auto t :
Range(crack_faces_tets)) {
6417 auto adj_faces = get_adj(
6419 moab::Interface::UNION);
6420 auto crack_surface_edges =
6421 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6424 moab::Interface::UNION);
6427 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6428 crack_surface_edges);
6429 adj_edges = subtract(adj_edges, body_skin_edges);
6431 auto tet_edges = get_adj(
Range(
t,
t), 1,
6432 moab::Interface::UNION);
6434 tet_edges = subtract(tet_edges, adj_edges);
6436 for (
auto e : tet_edges) {
6437 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
6438 auto get_side = [&](
auto e) {
6439 int side, sense, offset;
6442 "get side number failed");
6445 auto get_side_ent = [&](
auto side) {
6452 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
6455 if (adj_edges.size() <= max_adj_edges) {
6458 Range faces_to_remove;
6459 for (
auto e : adj_edges) {
6460 auto edge_adj_faces =
6461 get_adj(
Range(e, e), 2, moab::Interface::UNION);
6462 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
6463 if (edge_adj_faces.size() != 2) {
6465 "Adj faces size is not 2 for edge " +
6466 boost::lexical_cast<std::string>(e));
6469 auto get_normal = [&](
auto f) {
6473 "get tri normal failed");
6476 auto t_n0 = get_normal(edge_adj_faces[0]);
6477 auto t_n1 = get_normal(edge_adj_faces[1]);
6478 auto get_sense = [&](
auto f) {
6479 int side, sense, offset;
6482 "get side number failed");
6485 auto sense0 = get_sense(edge_adj_faces[0]);
6486 auto sense1 = get_sense(edge_adj_faces[1]);
6491 auto dot_e = (sense0 * sense1) * t_n0(
i) * t_n1(
i);
6492 if (dot_e < dot || e == adj_edges[0]) {
6494 faces_to_remove = edge_adj_faces;
6498 all_removed_faces.merge(faces_to_remove);
6499 all_removed_tets.merge(
Range(
t,
t));
6502 <<
"Remove free edges on flat tet, with considered nb. of "
6504 << adj_edges.size();
6508 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6509 crack_faces_tets_faces =
6510 subtract(crack_faces_tets_faces, all_removed_faces);
6516 "Case only one free edge failed");
6517 for (
auto max_adj_edges : {0, 1, 2, 3}) {
6519 "Case only one free edge failed");
6522 "Case internal faces failed");
6526 "crack_faces_tets_faces_" +
6527 boost::lexical_cast<std::string>(counter) +
".vtk",
6528 crack_faces_tets_faces);
6530 "crack_faces_tets_" +
6531 boost::lexical_cast<std::string>(counter) +
".vtk",
6535 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
6536 all_removed_faces, all_removed_tets);
6539 auto [resolved_faces, resolved_tets, all_removed_faces,
6541 resolve_surface(boundary_tets_edges, crack_faces_tets);
6542 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
6550 crack_faces = resolved_faces;
6562 auto resolve_consisten_crack_extension = [&]() {
6564 auto crack_meshset =
6567 auto meshset = crack_meshset->getMeshset();
6571 Range old_crack_faces;
6574 auto extendeded_crack_faces = get_extended_crack_faces();
6575 auto reconstructed_crack_faces =
6576 subtract(reconstruct_crack_faces(extendeded_crack_faces),
6578 if (
nbCrackFaces >= reconstructed_crack_faces.size()) {
6580 <<
"No new crack faces to add, skipping adding to meshset";
6581 extendeded_crack_faces = subtract(
6582 extendeded_crack_faces, subtract(*
crackFaces, old_crack_faces));
6584 <<
"Number crack faces size (extended) "
6585 << extendeded_crack_faces.size();
6591 reconstructed_crack_faces);
6593 <<
"Number crack faces size (reconstructed) "
6594 << reconstructed_crack_faces.size();
6613 CHKERR resolve_consisten_crack_extension();
6626 verts = subtract(verts, conn);
6627 std::vector<double> coords(3 * verts.size());
6629 double def_coords[] = {0., 0., 0.};
6632 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
6633 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
6653 auto post_proc_norm_fe =
6654 boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
6656 auto post_proc_norm_rule_hook = [](int, int,
int p) ->
int {
6659 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
6661 post_proc_norm_fe->getUserPolynomialBase() =
6664 CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
6668 enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
6671 CHKERR VecZeroEntries(norms_vec);
6673 auto u_l2_ptr = boost::make_shared<MatrixDouble>();
6674 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
6675 post_proc_norm_fe->getOpPtrVector().push_back(
6677 post_proc_norm_fe->getOpPtrVector().push_back(
6679 post_proc_norm_fe->getOpPtrVector().push_back(
6681 post_proc_norm_fe->getOpPtrVector().push_back(
6683 post_proc_norm_fe->getOpPtrVector().push_back(
6687 auto piola_ptr = boost::make_shared<MatrixDouble>();
6688 post_proc_norm_fe->getOpPtrVector().push_back(
6690 post_proc_norm_fe->getOpPtrVector().push_back(
6694 post_proc_norm_fe->getOpPtrVector().push_back(
6697 TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
6699 *post_proc_norm_fe);
6700 TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
6702 CHKERR VecAssemblyBegin(norms_vec);
6703 CHKERR VecAssemblyEnd(norms_vec);
6704 const double *norms;
6705 CHKERR VecGetArrayRead(norms_vec, &norms);
6706 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
6707 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
6709 <<
"norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
6711 <<
"norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
6712 CHKERR VecRestoreArrayRead(norms_vec, &norms);
6726 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6727 if (
auto disp_bc = bc.second->dispBcPtr) {
6732 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6733 MOFEM_LOG(
"EP", Sev::noisy) <<
"Displacement BC: " << *disp_bc;
6735 std::vector<double> block_attributes(6, 0.);
6736 if (disp_bc->data.flag1 == 1) {
6737 block_attributes[0] = disp_bc->data.value1;
6738 block_attributes[3] = 1;
6740 if (disp_bc->data.flag2 == 1) {
6741 block_attributes[1] = disp_bc->data.value2;
6742 block_attributes[4] = 1;
6744 if (disp_bc->data.flag3 == 1) {
6745 block_attributes[2] = disp_bc->data.value3;
6746 block_attributes[5] = 1;
6748 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6756 boost::make_shared<NormalDisplacementBcVec>();
6757 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6758 auto block_name =
"(.*)NORMAL_DISPLACEMENT(.*)";
6759 std::regex reg_name(block_name);
6760 if (std::regex_match(bc.first, reg_name)) {
6764 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6766 block_name, bc.second->bcAttributes,
6767 bc.second->bcEnts.subset_by_dimension(2));
6772 boost::make_shared<AnalyticalDisplacementBcVec>();
6774 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6775 auto block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
6776 std::regex reg_name(block_name);
6777 if (std::regex_match(bc.first, reg_name)) {
6781 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6783 block_name, bc.second->bcAttributes,
6784 bc.second->bcEnts.subset_by_dimension(2));
6788 auto ts_displacement =
6789 boost::make_shared<DynamicRelaxationTimeScale>(
"disp_history.txt");
6792 <<
"Add time scaling displacement BC: " << bc.blockName;
6795 ts_displacement,
"disp_history",
".txt", bc.blockName);
6798 auto ts_normal_displacement =
6799 boost::make_shared<DynamicRelaxationTimeScale>(
"normal_disp_history.txt");
6802 <<
"Add time scaling normal displacement BC: " << bc.blockName;
6805 ts_normal_displacement,
"normal_disp_history",
".txt",
6821 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6822 if (
auto force_bc = bc.second->forceBcPtr) {
6827 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6828 MOFEM_LOG(
"EP", Sev::noisy) <<
"Force BC: " << *force_bc;
6830 std::vector<double> block_attributes(6, 0.);
6831 block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
6832 block_attributes[3] = 1;
6833 block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
6834 block_attributes[4] = 1;
6835 block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
6836 block_attributes[5] = 1;
6837 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6845 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6846 auto block_name =
"(.*)PRESSURE(.*)";
6847 std::regex reg_name(block_name);
6848 if (std::regex_match(bc.first, reg_name)) {
6853 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6855 block_name, bc.second->bcAttributes,
6856 bc.second->bcEnts.subset_by_dimension(2));
6861 boost::make_shared<AnalyticalTractionBcVec>();
6863 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6864 auto block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
6865 std::regex reg_name(block_name);
6866 if (std::regex_match(bc.first, reg_name)) {
6870 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6872 block_name, bc.second->bcAttributes,
6873 bc.second->bcEnts.subset_by_dimension(2));
6878 boost::make_shared<DynamicRelaxationTimeScale>(
"traction_history.txt");
6882 ts_traction,
"traction_history",
".txt", bc.blockName);
6886 boost::make_shared<DynamicRelaxationTimeScale>(
"pressure_history.txt");
6890 ts_pressure,
"pressure_history",
".txt", bc.blockName);
6899 auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
6900 const std::string block_name,
6901 const int nb_attributes) {
6904 std::regex((boost::format(
"%s(.*)") % block_name).str()))
6906 std::vector<double> block_attributes;
6907 CHKERR it->getAttributes(block_attributes);
6908 if (block_attributes.size() < nb_attributes) {
6910 "In block %s expected %d attributes, but given %ld",
6911 it->getName().c_str(), nb_attributes, block_attributes.size());
6914 auto get_block_ents = [&]() {
6920 auto Ents = get_block_ents();
6921 ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
6930 auto ts_pre_stretch =
6931 boost::make_shared<DynamicRelaxationTimeScale>(
"externalstrain_history.txt");
6934 <<
"Add time scaling external strain: " << ext_strain_block.blockName;
6937 ts_pre_stretch,
"externalstrain_history",
".txt", ext_strain_block.blockName);
6946 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
6949 CHKERR VecGetLocalSize(
v.second, &size);
6951 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
6952 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( " << low
6953 <<
" " << high <<
" ) ";
6985 MOFEM_LOG(
"EP", Sev::inform) <<
"Calculate crack area";
6987 for (
auto f : crack_faces) {
6990 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
6992 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
6994 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.
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)
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 ...
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)
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(MFIELD, NAME, ENT, IT)
loop over all dofs from a moFEM field and particular field
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
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ElementsAndOps< SPACE_DIM >::SideEle SideEle
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
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
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
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x)
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
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)
static int nbJIntegralLevels
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
SetIntegrationAtFrontVolume(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges)
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
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.
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)
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