15#include <boost/math/constants/constants.hpp>
18#ifdef ENABLE_PYTHON_BINDING
19 #include <boost/python.hpp>
20 #include <boost/python/def.hpp>
21 #include <boost/python/numpy.hpp>
22namespace bp = boost::python;
23namespace np = boost::python::numpy;
25 #pragma message "With ENABLE_PYTHON_BINDING"
29 #pragma message "Without ENABLE_PYTHON_BINDING"
44 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
45 using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
48 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
49 using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
55 const EntityType type) {
59 auto dim = CN::Dimension(type);
61 std::vector<int> sendcounts(pcomm->size());
62 std::vector<int> displs(pcomm->size());
63 std::vector<int> sendbuf(r.size());
64 if (pcomm->rank() == 0) {
65 for (
auto p = 1; p != pcomm->size(); p++) {
67 ->getPartEntities(m_field.
get_moab(), p)
70 CHKERR m_field.
get_moab().get_adjacencies(part_ents, dim,
true, faces,
71 moab::Interface::UNION);
72 faces = intersect(faces, r);
73 sendcounts[p] = faces.size();
74 displs[p] = sendbuf.size();
75 for (
auto f : faces) {
77 sendbuf.push_back(
id);
83 MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
85 std::vector<int> recvbuf(recv_data);
86 MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
87 recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
89 if (pcomm->rank() > 0) {
91 for (
auto &f : recvbuf) {
101 const std::string block_name,
int dim) {
107 std::regex((boost::format(
"%s(.*)") % block_name).str())
111 for (
auto bc : bcs) {
123 const std::string block_name,
int dim) {
124 std::map<std::string, Range> r;
129 std::regex((boost::format(
"%s(.*)") % block_name).str())
133 for (
auto bc : bcs) {
138 r[bc->getName()] = faces;
145 const unsigned int cubit_bc_type) {
148 CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
152static auto save_range(moab::Interface &moab,
const std::string name,
153 const Range r, std::vector<Tag> tags = {}) {
156 CHKERR moab.add_entities(*out_meshset, r);
158 CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1,
159 tags.data(), tags.size());
161 MOFEM_LOG(
"SELF", Sev::warning) <<
"Empty range for " << name;
168 ParallelComm *pcomm =
171 PSTATUS_SHARED | PSTATUS_MULTISHARED,
172 PSTATUS_NOT, -1, &boundary_ents),
174 return boundary_ents;
179 ParallelComm *pcomm =
181 CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
190 CHK_MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents),
"find_skin");
196 ParallelComm *pcomm =
199 Range crack_skin_without_bdy;
200 if (pcomm->rank() == 0) {
202 CHKERR moab.get_adjacencies(crack_faces, 1,
true, crack_edges,
203 moab::Interface::UNION);
204 auto crack_skin =
get_skin(m_field, crack_faces);
208 "get_entities_by_dimension");
209 auto body_skin =
get_skin(m_field, body_ents);
210 Range body_skin_edges;
211 CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1,
true, body_skin_edges,
212 moab::Interface::UNION),
214 crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
216 for (
auto &
m : front_edges_map) {
217 auto add_front = subtract(
m.second, crack_edges);
218 auto i = intersect(
m.second, crack_edges);
220 crack_skin_without_bdy.merge(add_front);
224 CHKERR moab.get_adjacencies(i_skin, 1,
true, adj_i_skin,
225 moab::Interface::UNION);
226 adj_i_skin = subtract(intersect(adj_i_skin,
m.second), crack_edges);
227 crack_skin_without_bdy.merge(adj_i_skin);
231 return send_type(m_field, crack_skin_without_bdy, MBEDGE);
237 ParallelComm *pcomm =
240 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface";
242 if (!pcomm->rank()) {
244 auto impl = [&](
auto &saids) {
249 auto get_adj = [&](
auto e,
auto dim) {
252 e, dim,
true, adj, moab::Interface::UNION),
257 auto get_conn = [&](
auto e) {
264 constexpr bool debug =
false;
269 auto body_skin_edges = get_adj(body_skin, 1);
272 subtract(
get_skin(m_field, crack_faces), body_skin_edges);
273 auto crack_skin_conn = get_conn(crack_skin);
274 auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
275 auto crack_edges = get_adj(crack_faces, 1);
276 crack_edges = subtract(crack_edges, crack_skin);
277 auto all_tets = get_adj(crack_edges, 3);
278 crack_edges = subtract(crack_edges, crack_skin_conn_edges);
279 auto crack_conn = get_conn(crack_edges);
280 all_tets.merge(get_adj(crack_conn, 3));
289 if (crack_faces.size()) {
290 auto grow = [&](
auto r) {
291 auto crack_faces_conn = get_conn(crack_faces);
294 while (size_r != r.size() && r.size() > 0) {
296 CHKERR moab.get_connectivity(r,
v,
true);
297 v = subtract(
v, crack_faces_conn);
300 moab::Interface::UNION);
301 r = intersect(r, all_tets);
310 Range all_tets_ord = all_tets;
311 while (all_tets.size()) {
312 Range faces = get_adj(unite(saids.first, saids.second), 2);
313 faces = subtract(crack_faces, faces);
316 auto fit = faces.begin();
317 for (; fit != faces.end(); ++fit) {
318 tets = intersect(get_adj(
Range(*fit, *fit), 3), all_tets);
319 if (tets.size() == 2) {
326 saids.first.insert(tets[0]);
327 saids.first = grow(saids.first);
328 all_tets = subtract(all_tets, saids.first);
329 if (tets.size() == 2) {
330 saids.second.insert(tets[1]);
331 saids.second = grow(saids.second);
332 all_tets = subtract(all_tets, saids.second);
340 saids.first = subtract(all_tets_ord, saids.second);
341 saids.second = subtract(all_tets_ord, saids.first);
347 std::pair<Range, Range> saids;
348 if (crack_faces.size())
353 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface <- done";
355 return std::pair<Range, Range>();
363 boost::shared_ptr<Range> front_edges)
367 int order_col,
int order_data) {
370 constexpr bool debug =
false;
372 constexpr int numNodes = 4;
373 constexpr int numEdges = 6;
374 constexpr int refinementLevels = 4;
376 auto &m_field = fe_raw_ptr->
mField;
377 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
380 auto set_base_quadrature = [&]() {
382 int rule = 2 * order_data + 1;
393 auto &gauss_pts = fe_ptr->gaussPts;
394 gauss_pts.resize(4, nb_gauss_pts,
false);
395 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
396 &gauss_pts(0, 0), 1);
397 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
398 &gauss_pts(1, 0), 1);
399 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
400 &gauss_pts(2, 0), 1);
402 &gauss_pts(3, 0), 1);
403 auto &data = fe_ptr->dataOnElement[
H1];
404 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
407 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
408 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
417 CHKERR set_base_quadrature();
421 auto get_singular_nodes = [&]() {
424 CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
426 std::bitset<numNodes> singular_nodes;
427 for (
auto nn = 0; nn != numNodes; ++nn) {
429 singular_nodes.set(nn);
431 singular_nodes.reset(nn);
434 return singular_nodes;
437 auto get_singular_edges = [&]() {
438 std::bitset<numEdges> singular_edges;
439 for (
int ee = 0; ee != numEdges; ee++) {
443 singular_edges.set(ee);
445 singular_edges.reset(ee);
448 return singular_edges;
451 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
453 fe_ptr->gaussPts.swap(ref_gauss_pts);
454 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
455 auto &data = fe_ptr->dataOnElement[
H1];
456 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
458 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
460 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
465 auto singular_nodes = get_singular_nodes();
466 if (singular_nodes.count()) {
467 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
469 CHKERR set_gauss_pts(it_map_ref_coords->second);
473 auto refine_quadrature = [&]() {
476 const int max_level = refinementLevels;
480 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
482 for (
int nn = 0; nn != 4; nn++)
483 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
484 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
488 Range tets(tet, tet);
491 tets, 1,
true, edges, moab::Interface::UNION);
496 Range nodes_at_front;
497 for (
int nn = 0; nn != numNodes; nn++) {
498 if (singular_nodes[nn]) {
500 CHKERR moab_ref.side_element(tet, 0, nn, ent);
501 nodes_at_front.insert(ent);
505 auto singular_edges = get_singular_edges();
508 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
509 for (
int ee = 0; ee != numEdges; ee++) {
510 if (singular_edges[ee]) {
512 CHKERR moab_ref.side_element(tet, 1, ee, ent);
513 CHKERR moab_ref.add_entities(meshset, &ent, 1);
519 for (
int ll = 0; ll != max_level; ll++) {
522 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
526 CHKERR moab_ref.get_adjacencies(
527 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
528 ref_edges = intersect(ref_edges, edges);
530 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
531 ref_edges = intersect(ref_edges, ents);
534 ->getEntitiesByTypeAndRefLevel(
536 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
540 ->updateMeshsetByEntitiesChildren(meshset,
542 meshset, MBEDGE,
true);
548 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
556 MatrixDouble ref_coords(tets.size(), 12,
false);
558 for (Range::iterator tit = tets.begin(); tit != tets.end();
562 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
563 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
566 auto &data = fe_ptr->dataOnElement[
H1];
567 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
568 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
569 MatrixDouble &shape_n =
570 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
572 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
573 double *tet_coords = &ref_coords(tt, 0);
576 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
577 for (
int dd = 0; dd != 3; dd++) {
578 ref_gauss_pts(dd, gg) =
579 shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
580 shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
581 shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
582 shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
584 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
588 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
594 CHKERR refine_quadrature();
604 using ForcesAndSourcesCore::dataOnElement;
607 using ForcesAndSourcesCore::ForcesAndSourcesCore;
622 boost::shared_ptr<Range> front_edges)
626 boost::shared_ptr<Range> front_edges,
631 int order_col,
int order_data) {
634 constexpr bool debug =
false;
636 constexpr int numNodes = 3;
637 constexpr int numEdges = 3;
638 constexpr int refinementLevels = 4;
640 auto &m_field = fe_raw_ptr->
mField;
641 auto fe_ptr =
static_cast<Fe *
>(fe_raw_ptr);
644 auto set_base_quadrature = [&]() {
646 int rule =
funRule(order_data);
656 fe_ptr->gaussPts.resize(3, nb_gauss_pts,
false);
657 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
658 &fe_ptr->gaussPts(0, 0), 1);
659 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
660 &fe_ptr->gaussPts(1, 0), 1);
662 &fe_ptr->gaussPts(2, 0), 1);
663 auto &data = fe_ptr->dataOnElement[
H1];
664 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
667 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
668 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
670 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
673 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
682 CHKERR set_base_quadrature();
686 auto get_singular_nodes = [&]() {
689 CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
691 std::bitset<numNodes> singular_nodes;
692 for (
auto nn = 0; nn != numNodes; ++nn) {
694 singular_nodes.set(nn);
696 singular_nodes.reset(nn);
699 return singular_nodes;
702 auto get_singular_edges = [&]() {
703 std::bitset<numEdges> singular_edges;
704 for (
int ee = 0; ee != numEdges; ee++) {
708 singular_edges.set(ee);
710 singular_edges.reset(ee);
713 return singular_edges;
716 auto set_gauss_pts = [&](
auto &ref_gauss_pts) {
718 fe_ptr->gaussPts.swap(ref_gauss_pts);
719 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
720 auto &data = fe_ptr->dataOnElement[
H1];
721 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
723 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
725 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
729 auto singular_nodes = get_singular_nodes();
730 if (singular_nodes.count()) {
731 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
733 CHKERR set_gauss_pts(it_map_ref_coords->second);
737 auto refine_quadrature = [&]() {
740 const int max_level = refinementLevels;
743 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
745 for (
int nn = 0; nn != numNodes; nn++)
746 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
748 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
752 Range tris(tri, tri);
755 tris, 1,
true, edges, moab::Interface::UNION);
760 Range nodes_at_front;
761 for (
int nn = 0; nn != numNodes; nn++) {
762 if (singular_nodes[nn]) {
764 CHKERR moab_ref.side_element(tri, 0, nn, ent);
765 nodes_at_front.insert(ent);
769 auto singular_edges = get_singular_edges();
772 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
773 for (
int ee = 0; ee != numEdges; ee++) {
774 if (singular_edges[ee]) {
776 CHKERR moab_ref.side_element(tri, 1, ee, ent);
777 CHKERR moab_ref.add_entities(meshset, &ent, 1);
783 for (
int ll = 0; ll != max_level; ll++) {
786 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
790 CHKERR moab_ref.get_adjacencies(
791 nodes_at_front, 1,
true, ref_edges, moab::Interface::UNION);
792 ref_edges = intersect(ref_edges, edges);
794 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
795 ref_edges = intersect(ref_edges, ents);
798 ->getEntitiesByTypeAndRefLevel(
800 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
804 ->updateMeshsetByEntitiesChildren(meshset,
806 meshset, MBEDGE,
true);
812 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
821 MatrixDouble ref_coords(tris.size(), 9,
false);
823 for (Range::iterator tit = tris.begin(); tit != tris.end();
827 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
828 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
831 auto &data = fe_ptr->dataOnElement[
H1];
832 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
833 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
834 MatrixDouble &shape_n =
835 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
837 for (
size_t tt = 0; tt != ref_coords.size1(); tt++) {
838 double *tri_coords = &ref_coords(tt, 0);
841 auto det = t_normal.
l2();
842 for (
size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
843 for (
int dd = 0; dd != 2; dd++) {
844 ref_gauss_pts(dd, gg) =
845 shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
846 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
847 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
849 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
853 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
859 CHKERR refine_quadrature();
869 using ForcesAndSourcesCore::dataOnElement;
872 using ForcesAndSourcesCore::ForcesAndSourcesCore;
921 const char *list_rots[] = {
"small",
"moderate",
"large",
"no_h1"};
922 const char *list_symm[] = {
"symm",
"not_symm"};
923 const char *list_release[] = {
"griffith_force",
"griffith_skelton"};
924 const char *list_stretches[] = {
"linear",
"log",
"log_quadratic"};
930 char analytical_expr_file_name[255] =
"analytical_expr.py";
932 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
934 CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
936 CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
938 CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
940 CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
942 CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
944 CHKERR PetscOptionsScalar(
"-viscosity_alpha_omega",
"rot viscosity",
"",
946 CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
948 CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
949 LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
951 CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
953 list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
954 CHKERR PetscOptionsEList(
"-symm",
"symmetric variant",
"", list_symm, 2,
955 list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
959 CHKERR PetscOptionsEList(
"-stretches",
"stretches",
"", list_stretches,
961 list_stretches[choice_stretch], &choice_stretch,
964 CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
966 CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
968 CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
972 CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
976 CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
983 CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
985 CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
987 CHKERR PetscOptionsEList(
"-energy_release_variant",
"energy release variant",
988 "", list_release, 2, list_release[choice_release],
989 &choice_release, PETSC_NULLPTR);
990 CHKERR PetscOptionsInt(
"-nb_J_integral_levels",
"Number of J integarl levels",
994 char tag_name[255] =
"";
995 CHKERR PetscOptionsString(
"-internal_stress_tag_name",
996 "internal stress tag name",
"",
"", tag_name, 255,
1000 "-internal_stress_interp_order",
"internal stress interpolation order",
1002 CHKERR PetscOptionsBool(
"-internal_stress_voigt",
"Voigt index notation",
"",
1007 analytical_expr_file_name, 255, PETSC_NULLPTR);
1013 "Unsupported internal stress interpolation order %d",
1039 std::numeric_limits<float>::epsilon()) {
1061 "No logarithmic quadratic stretch for this case");
1077 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaU: -viscosity_alpha_u " <<
alphaU;
1078 MOFEM_LOG(
"EP", Sev::inform) <<
"alphaW: -viscosity_alpha_w " <<
alphaW;
1080 <<
"alphaOmega: -viscosity_alpha_omega " <<
alphaOmega;
1084 MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
1092 MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
1094 <<
"Stretch: -stretches " << list_stretches[choice_stretch];
1120 <<
"Energy release variant: -energy_release_variant "
1123 <<
"Number of J integral levels: -nb_J_integral_levels "
1126#ifdef ENABLE_PYTHON_BINDING
1127 auto file_exists = [](std::string myfile) {
1128 std::ifstream file(myfile.c_str());
1135 if (file_exists(analytical_expr_file_name)) {
1136 MOFEM_LOG(
"EP", Sev::inform) << analytical_expr_file_name <<
" file found";
1140 analytical_expr_file_name);
1144 << analytical_expr_file_name <<
" file NOT found";
1157 auto get_tets = [&]() {
1163 auto get_tets_skin = [&]() {
1164 Range tets_skin_part;
1166 CHKERR skin.find_skin(0, get_tets(),
false, tets_skin_part);
1167 ParallelComm *pcomm =
1170 CHKERR pcomm->filter_pstatus(tets_skin_part,
1171 PSTATUS_SHARED | PSTATUS_MULTISHARED,
1172 PSTATUS_NOT, -1, &tets_skin);
1176 auto subtract_boundary_conditions = [&](
auto &&tets_skin) {
1182 tets_skin = subtract(tets_skin,
v.faces);
1186 tets_skin = subtract(tets_skin,
v.faces);
1191 tets_skin = subtract(tets_skin,
v.faces);
1197 auto add_blockset = [&](
auto block_name,
auto &&tets_skin) {
1200 tets_skin.merge(crack_faces);
1204 auto subtract_blockset = [&](
auto block_name,
auto &&tets_skin) {
1205 auto contact_range =
1207 tets_skin = subtract(tets_skin, contact_range);
1211 auto get_stress_trace_faces = [&](
auto &&tets_skin) {
1214 faces, moab::Interface::UNION);
1215 Range trace_faces = subtract(faces, tets_skin);
1219 auto tets = get_tets();
1223 auto trace_faces = get_stress_trace_faces(
1225 subtract_blockset(
"CONTACT",
1226 subtract_boundary_conditions(get_tets_skin()))
1233 boost::make_shared<Range>(subtract(trace_faces, *
contactFaces));
1253 auto add_broken_hdiv_field = [
this, meshset](
const std::string
field_name,
1259 auto get_side_map_hdiv = [&]() {
1262 std::pair<EntityType,
1277 get_side_map_hdiv(), MB_TAG_DENSE,
MF_ZERO);
1283 auto add_l2_field = [
this, meshset](
const std::string
field_name,
1284 const int order,
const int dim) {
1293 auto add_h1_field = [
this, meshset](
const std::string
field_name,
1294 const int order,
const int dim) {
1306 auto add_l2_field_by_range = [
this](
const std::string
field_name,
1307 const int order,
const int dim,
1308 const int field_dim,
Range &&r) {
1318 auto add_bubble_field = [
this, meshset](
const std::string
field_name,
1319 const int order,
const int dim) {
1325 auto field_order_table =
1326 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1327 auto get_cgg_bubble_order_zero = [](
int p) {
return 0; };
1328 auto get_cgg_bubble_order_tet = [](
int p) {
1331 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
1332 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
1333 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
1334 field_order_table[MBTET] = get_cgg_bubble_order_tet;
1341 auto add_user_l2_field = [
this, meshset](
const std::string
field_name,
1342 const int order,
const int dim) {
1348 auto field_order_table =
1349 const_cast<Field *
>(field_ptr)->getFieldOrderTable();
1350 auto zero_dofs = [](
int p) {
return 0; };
1352 field_order_table[MBVERTEX] = zero_dofs;
1353 field_order_table[MBEDGE] = zero_dofs;
1354 field_order_table[MBTRI] = zero_dofs;
1355 field_order_table[MBTET] = dof_l2_tet;
1373 auto get_hybridised_disp = [&]() {
1375 auto skin = subtract_boundary_conditions(get_tets_skin());
1377 faces.merge(intersect(bc.faces, skin));
1383 get_hybridised_disp());
1409 auto project_ho_geometry = [&](
auto field) {
1415 auto get_adj_front_edges = [&](
auto &front_edges) {
1416 Range front_crack_nodes;
1417 Range crack_front_edges_with_both_nodes_not_at_front;
1421 CHKERR moab.get_connectivity(front_edges, front_crack_nodes,
true);
1422 Range crack_front_edges;
1424 crack_front_edges, moab::Interface::UNION);
1426 intersect(subtract(crack_front_edges, front_edges), meshset_ents);
1427 Range crack_front_edges_nodes;
1428 CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
1430 crack_front_edges_nodes =
1431 subtract(crack_front_edges_nodes, front_crack_nodes);
1432 CHKERR moab.get_adjacencies(
1433 crack_front_edges_nodes,
SPACE_DIM - 2,
false,
1434 crack_front_edges_with_both_nodes_not_at_front,
1435 moab::Interface::UNION);
1436 crack_front_edges_with_both_nodes_not_at_front = intersect(
1437 crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
1438 crack_front_edges_with_both_nodes_not_at_front = intersect(
1439 crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
1443 crack_front_edges_with_both_nodes_not_at_front =
send_type(
1444 mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
1446 return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
1447 boost::make_shared<Range>(
1448 crack_front_edges_with_both_nodes_not_at_front));
1455 auto [front_vertices, front_adj_edges] = get_adj_front_edges(*
frontEdges);
1466 (boost::format(
"crack_faces_%d.vtk") % rank).str(),
1469 (boost::format(
"front_edges_%d.vtk") % rank).str(),
1480 auto set_singular_dofs = [&](
auto &front_adj_edges,
auto &front_vertices) {
1488 MOFEM_LOG(
"EP", Sev::inform) <<
"Singularity eps " << beta;
1491 for (
auto edge : front_adj_edges) {
1494 CHKERR moab.get_connectivity(edge, conn, num_nodes,
false);
1496 CHKERR moab.get_coords(conn, num_nodes, coords);
1497 const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
1498 coords[5] - coords[2]};
1499 double dof[3] = {0, 0, 0};
1500 if (front_vertices.find(conn[0]) != front_vertices.end()) {
1501 for (
int dd = 0; dd != 3; dd++) {
1502 dof[dd] = -dir[dd] *
eps;
1504 }
else if (front_vertices.find(conn[1]) != front_vertices.end()) {
1505 for (
int dd = 0; dd != 3; dd++) {
1506 dof[dd] = +dir[dd] *
eps;
1511 const int idx = dit->get()->getEntDofIdx();
1513 dit->get()->getFieldData() = 0;
1515 dit->get()->getFieldData() = dof[idx];
1534 auto add_field_to_fe = [
this](
const std::string fe,
1599 auto set_fe_adjacency = [&](
auto fe_name) {
1602 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
1610 auto add_field_to_fe = [
this](
const std::string fe,
1622 Range natural_bc_elements;
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);
1650 natural_bc_elements.merge(
v.faces);
1655 natural_bc_elements.merge(
v.faces);
1658 natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
1669 auto get_skin = [&](
auto &body_ents) {
1672 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
1677 Range boundary_ents;
1678 ParallelComm *pcomm =
1680 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
1681 PSTATUS_NOT, -1, &boundary_ents);
1682 return boundary_ents;
1816 auto remove_dofs_on_broken_skin = [&](
const std::string prb_name) {
1818 for (
int d : {0, 1, 2}) {
1819 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1821 ->getSideDofsOnBrokenSpaceEntities(
1832 CHKERR remove_dofs_on_broken_skin(
"ESHELBY_PLASTICITY");
1868 auto set_zero_block = [&]() {
1900 auto set_section = [&]() {
1902 PetscSection section;
1907 CHKERR PetscSectionDestroy(§ion);
1931 : blockName(name), faces(faces) {
1932 vals.resize(3,
false);
1933 flags.resize(3,
false);
1934 for (
int ii = 0; ii != 3; ++ii) {
1935 vals[ii] = attr[ii];
1936 flags[ii] =
static_cast<int>(attr[ii + 3]);
1939 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
1941 <<
"Add BCDisp vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1943 <<
"Add BCDisp flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1944 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " <<
faces.size();
1948 : blockName(name), faces(faces) {
1949 vals.resize(attr.size(),
false);
1950 for (
int ii = 0; ii != attr.size(); ++ii) {
1951 vals[ii] = attr[ii];
1957 : blockName(name), faces(faces) {
1958 vals.resize(3,
false);
1959 flags.resize(3,
false);
1960 for (
int ii = 0; ii != 3; ++ii) {
1961 vals[ii] = attr[ii];
1962 flags[ii] =
static_cast<int>(attr[ii + 3]);
1965 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
1967 <<
"Add BCForce vals " <<
vals[0] <<
" " <<
vals[1] <<
" " <<
vals[2];
1969 <<
"Add BCForce flags " <<
flags[0] <<
" " <<
flags[1] <<
" " <<
flags[2];
1970 MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " <<
faces.size();
1974 std::vector<double> attr,
1976 : blockName(name), faces(faces) {
1979 if (attr.size() < 1) {
1980 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1981 "Wrong size of normal displacement BC");
1986 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc " << name;
1987 MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc val " <<
val;
1989 <<
"Add NormalDisplacementBc nb. of faces " <<
faces.size();
1993 : blockName(name), faces(faces) {
1996 if (attr.size() < 1) {
1997 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1998 "Wrong size of normal displacement BC");
2003 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc " << name;
2004 MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc val " <<
val;
2006 <<
"Add PressureBc nb. of faces " <<
faces.size();
2012 : blockName(name), ents(ents) {
2015 if (attr.size() < 2) {
2016 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
2017 "Wrong size of external strain attribute");
2023 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain " << name;
2024 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain val " <<
val;
2025 MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain bulk modulus K "
2028 <<
"Add ExternalStrain nb. of tets " <<
ents.size();
2034 std::vector<double> attr,
2036 : blockName(name), faces(faces) {
2039 if (attr.size() < 3) {
2040 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
2041 "Wrong size of analytical displacement BC");
2044 flags.resize(3,
false);
2045 for (
int ii = 0; ii != 3; ++ii) {
2046 flags[ii] = attr[ii];
2049 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalDisplacementBc " << name;
2051 <<
"Add AnalyticalDisplacementBc flags " <<
flags[0] <<
" " <<
flags[1]
2054 <<
"Add AnalyticalDisplacementBc nb. of faces " <<
faces.size();
2058 std::vector<double> attr,
2060 : blockName(name), faces(faces) {
2063 if (attr.size() < 3) {
2064 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
2065 "Wrong size of analytical traction BC");
2068 flags.resize(3,
false);
2069 for (
int ii = 0; ii != 3; ++ii) {
2070 flags[ii] = attr[ii];
2073 MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalTractionBc " << name;
2075 <<
"Add AnalyticalTractionBc flags " <<
flags[0] <<
" " <<
flags[1]
2078 <<
"Add AnalyticalTractionBc nb. of faces " <<
faces.size();
2084 boost::shared_ptr<TractionFreeBc> &bc_ptr,
2085 const std::string contact_set_name) {
2091 Range tets_skin_part;
2093 CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
2094 ParallelComm *pcomm =
2097 CHKERR pcomm->filter_pstatus(tets_skin_part,
2098 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2099 PSTATUS_NOT, -1, &tets_skin);
2102 for (
int dd = 0; dd != 3; ++dd)
2103 (*bc_ptr)[dd] = tets_skin;
2109 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2111 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2113 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2119 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2120 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2121 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2126 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2127 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2128 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2134 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2136 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2138 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2143 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2144 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2145 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2150 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2151 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2152 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2157 (*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
2158 (*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
2159 (*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
2164 std::regex((boost::format(
"%s(.*)") % contact_set_name).str()))) {
2168 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
2169 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
2170 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
2187 return 2 * p_data + 1;
2193 return 2 * (p_data + 1);
2198 boost::shared_ptr<CachePhi> cache_phi_otr)
2202 boost::typeindex::type_index type_index,
2210 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
2215 int nb_gauss_pts = pts.size2();
2216 if (!nb_gauss_pts) {
2220 if (pts.size1() < 3) {
2222 "Wrong dimension of pts, should be at least 3 rows with "
2226 const auto base = this->cTx->bAse;
2229 switch (this->cTx->sPace) {
2234 data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
false);
2237 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
2238 data.
dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3,
false);
2242 CHKERR getValueL2AinsworthBase(pts);
2258 "Wrong base, should be USER_BASE");
2266 const int nb_gauss_pts = pts.size2();
2268 auto check_cache = [
this](
int order,
int nb_gauss_pts) ->
bool {
2277 if (check_cache(
order, nb_gauss_pts)) {
2280 phi.resize(mat.size1(), mat.size2(),
false);
2285 shapeFun.resize(nb_gauss_pts, 4,
false);
2287 &pts(2, 0), nb_gauss_pts);
2289 double diff_shape_fun[12];
2295 phi.resize(nb_gauss_pts, 9 * nb_base_functions,
false);
2318 const int tag,
const bool do_rhs,
const bool do_lhs,
const bool calc_rates,
2320 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
2324 boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0,
MatrixDouble());
2325 fe->getUserPolynomialBase() =
2326 boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
2331 fe->getRuleHook = [](int, int, int) {
return -1; };
2350 fe->getOpPtrVector().push_back(
2354 fe->getOpPtrVector().push_back(
2378 fe->getOpPtrVector().push_back(
2381 fe->getOpPtrVector().push_back(
2392 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
2405 var_vec, MBMAXTYPE));
2414 fe->getOpPtrVector().push_back(
2418 fe->getOpPtrVector().push_back(
2439 const int tag,
const bool add_elastic,
const bool add_material,
2440 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
2441 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
2445 auto get_body_range = [
this](
auto name,
int dim) {
2446 std::map<int, Range> map;
2451 (boost::format(
"%s(.*)") % name).str()
2460 map[m_ptr->getMeshsetId()] = ents;
2466 auto rule_contact = [](int, int,
int o) {
return -1; };
2469 auto set_rule_contact = [refine](
2472 int order_col,
int order_data
2476 auto rule = 2 * order_data;
2482 fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
2489 fe_rhs->getOpPtrVector().push_back(
2491 fe_rhs->getOpPtrVector().push_back(
2499 fe_rhs->getOpPtrVector().push_back(
2503 fe_rhs->getOpPtrVector().push_back(
2508 "Unsupported internal stress interpolation order %d",
2512 fe_rhs->getOpPtrVector().push_back(
2516 fe_rhs->getOpPtrVector().push_back(
2523 fe_rhs->getOpPtrVector().push_back(op);
2526 "OpSpatialPhysicalExternalStrain not implemented for this "
2530 fe_rhs->getOpPtrVector().push_back(
2534 fe_rhs->getOpPtrVector().push_back(
2536 fe_rhs->getOpPtrVector().push_back(
2538 fe_rhs->getOpPtrVector().push_back(
2541 auto set_hybridisation = [&](
auto &pip) {
2548 using SideEleOp = EleOnSide::UserDataOperator;
2549 using BdyEleOp = BoundaryEle::UserDataOperator;
2556 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2559 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2562 CHKERR EshelbianPlasticity::
2563 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2564 op_loop_skeleton_side->getOpPtrVector(), {L2},
2569 auto broken_data_ptr =
2570 boost::make_shared<std::vector<BrokenBaseSideData>>();
2574 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2575 boost::make_shared<CGGUserPolynomialBase>();
2578 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2580 op_loop_domain_side->getOpPtrVector().push_back(
2582 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
2583 op_loop_domain_side->getOpPtrVector().push_back(
2586 op_loop_domain_side->getOpPtrVector().push_back(
2590 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2592 GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
2594 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2595 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dHybrid(
2597 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2598 op_loop_skeleton_side->getOpPtrVector().push_back(
2601 op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dBroken(
2602 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
2605 pip.push_back(op_loop_skeleton_side);
2610 auto set_contact = [&](
auto &pip) {
2617 using SideEleOp = EleOnSide::UserDataOperator;
2618 using BdyEleOp = BoundaryEle::UserDataOperator;
2625 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2626 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2627 CHKERR EshelbianPlasticity::
2628 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2629 op_loop_skeleton_side->getOpPtrVector(), {L2},
2634 auto broken_data_ptr =
2635 boost::make_shared<std::vector<BrokenBaseSideData>>();
2638 auto contact_common_data_ptr =
2639 boost::make_shared<ContactOps::CommonData>();
2641 auto add_ops_domain_side = [&](
auto &pip) {
2646 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2647 boost::make_shared<CGGUserPolynomialBase>();
2650 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2652 op_loop_domain_side->getOpPtrVector().push_back(
2655 op_loop_domain_side->getOpPtrVector().push_back(
2657 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2658 pip.push_back(op_loop_domain_side);
2662 auto add_ops_contact_rhs = [&](
auto &pip) {
2665 auto contact_sfd_map_range_ptr =
2666 boost::make_shared<std::map<int, Range>>(
2667 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2670 contactDisp, contact_common_data_ptr->contactDispPtr()));
2671 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2680 contact_sfd_map_range_ptr));
2689 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2690 CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
2693 pip.push_back(op_loop_skeleton_side);
2698 CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
2699 CHKERR set_contact(fe_rhs->getOpPtrVector());
2702 using BodyNaturalBC =
2704 Assembly<PETSC>::LinearForm<
GAUSS>;
2706 BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
2708 auto body_time_scale =
2709 boost::make_shared<DynamicRelaxationTimeScale>(
"body_force.txt");
2710 CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
2712 "BODY_FORCE", Sev::inform);
2716 fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
2724 fe_lhs->getOpPtrVector().push_back(
2728 fe_lhs->getOpPtrVector().push_back(
2732 fe_lhs->getOpPtrVector().push_back(
2773 auto set_hybridisation = [&](
auto &pip) {
2780 using SideEleOp = EleOnSide::UserDataOperator;
2781 using BdyEleOp = BoundaryEle::UserDataOperator;
2788 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2791 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2793 CHKERR EshelbianPlasticity::
2794 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2795 op_loop_skeleton_side->getOpPtrVector(), {L2},
2800 auto broken_data_ptr =
2801 boost::make_shared<std::vector<BrokenBaseSideData>>();
2805 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2806 boost::make_shared<CGGUserPolynomialBase>();
2809 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2811 op_loop_domain_side->getOpPtrVector().push_back(
2814 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2816 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2817 op_loop_skeleton_side->getOpPtrVector().push_back(
2819 boost::make_shared<double>(1.0),
true,
false));
2821 pip.push_back(op_loop_skeleton_side);
2826 auto set_contact = [&](
auto &pip) {
2833 using SideEleOp = EleOnSide::UserDataOperator;
2834 using BdyEleOp = BoundaryEle::UserDataOperator;
2841 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2842 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2843 CHKERR EshelbianPlasticity::
2844 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2845 op_loop_skeleton_side->getOpPtrVector(), {L2},
2850 auto broken_data_ptr =
2851 boost::make_shared<std::vector<BrokenBaseSideData>>();
2854 auto contact_common_data_ptr =
2855 boost::make_shared<ContactOps::CommonData>();
2857 auto add_ops_domain_side = [&](
auto &pip) {
2862 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2863 boost::make_shared<CGGUserPolynomialBase>();
2866 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2868 op_loop_domain_side->getOpPtrVector().push_back(
2871 op_loop_domain_side->getOpPtrVector().push_back(
2873 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2874 pip.push_back(op_loop_domain_side);
2878 auto add_ops_contact_lhs = [&](
auto &pip) {
2881 contactDisp, contact_common_data_ptr->contactDispPtr()));
2882 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2891 auto contact_sfd_map_range_ptr =
2892 boost::make_shared<std::map<int, Range>>(
2893 get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
2897 contact_sfd_map_range_ptr));
2900 contactDisp, broken_data_ptr, contact_common_data_ptr,
2904 broken_data_ptr,
contactDisp, contact_common_data_ptr,
2911 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2912 CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
2915 pip.push_back(op_loop_skeleton_side);
2920 CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
2921 CHKERR set_contact(fe_lhs->getOpPtrVector());
2931 const bool add_elastic,
const bool add_material,
2932 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
2933 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
2936 fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
2937 fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
2942 fe_rhs->getRuleHook = [](int, int, int) {
return -1; };
2943 fe_lhs->getRuleHook = [](int, int, int) {
return -1; };
2956 auto get_broken_op_side = [
this](
auto &pip) {
2959 using SideEleOp = EleOnSide::UserDataOperator;
2961 auto broken_data_ptr =
2962 boost::make_shared<std::vector<BrokenBaseSideData>>();
2966 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2967 boost::make_shared<CGGUserPolynomialBase>();
2970 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2972 op_loop_domain_side->getOpPtrVector().push_back(
2974 boost::shared_ptr<double> piola_scale_ptr(
new double);
2975 op_loop_domain_side->getOpPtrVector().push_back(
2979 auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
2980 op_loop_domain_side->getOpPtrVector().push_back(
2982 piola_stress_mat_ptr));
2983 pip.push_back(op_loop_domain_side);
2984 return std::make_tuple(broken_data_ptr, piola_scale_ptr,
2985 piola_stress_mat_ptr);
2988 auto set_rhs = [&]() {
2991 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
2992 get_broken_op_side(fe_rhs->getOpPtrVector());
2994 fe_rhs->getOpPtrVector().push_back(
3002 fe_rhs->getOpPtrVector().push_back(
3005 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
3007 fe_rhs->getOpPtrVector().push_back(
3017 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
3018 fe_rhs->getOpPtrVector().push_back(
3025 auto get_normal_disp_bc_faces = [&]() {
3028 return boost::make_shared<Range>(faces);
3033 using BdyEleOp = BoundaryEle::UserDataOperator;
3035 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
3036 fe_rhs->getOpPtrVector().push_back(
new OpC_dBroken(
3037 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
3038 get_normal_disp_bc_faces()));
3043 auto set_lhs = [&]() {
3046 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
3047 get_broken_op_side(fe_lhs->getOpPtrVector());
3055 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
3057 fe_rhs->getOpPtrVector().push_back(
3064 auto get_normal_disp_bc_faces = [&]() {
3067 return boost::make_shared<Range>(faces);
3072 using BdyEleOp = BoundaryEle::UserDataOperator;
3074 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
3075 fe_lhs->getOpPtrVector().push_back(
new OpC(
3077 true,
true, get_normal_disp_bc_faces()));
3091 boost::shared_ptr<ContactTree> &fe_contact_tree
3097 auto get_body_range = [
this](
auto name,
int dim,
auto sev) {
3098 std::map<int, Range> map;
3103 (boost::format(
"%s(.*)") % name).str()
3112 map[m_ptr->getMeshsetId()] = ents;
3113 MOFEM_LOG(
"EPSYNC", sev) <<
"Meshset: " << m_ptr->getMeshsetId() <<
" "
3114 << ents.size() <<
" entities";
3121 auto get_map_skin = [
this](
auto &&map) {
3122 ParallelComm *pcomm =
3126 for (
auto &
m : map) {
3128 CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
3130 PSTATUS_SHARED | PSTATUS_MULTISHARED,
3131 PSTATUS_NOT, -1,
nullptr),
3133 m.second.swap(skin_faces);
3143 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3145 auto calcs_side_traction = [&](
auto &pip) {
3149 using SideEleOp = EleOnSide::UserDataOperator;
3152 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3153 boost::make_shared<CGGUserPolynomialBase>();
3155 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3157 op_loop_domain_side->getOpPtrVector().push_back(
3159 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3160 boost::make_shared<double>(1.0)));
3161 pip.push_back(op_loop_domain_side);
3165 auto add_contact_three = [&]() {
3167 auto tree_moab_ptr = boost::make_shared<moab::Core>();
3168 fe_contact_tree = boost::make_shared<ContactTree>(
3170 get_body_range(
"CONTACT",
SPACE_DIM - 1, Sev::inform));
3171 fe_contact_tree->getOpPtrVector().push_back(
3173 contactDisp, contact_common_data_ptr->contactDispPtr()));
3174 CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
3175 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3176 fe_contact_tree->getOpPtrVector().push_back(
3178 fe_contact_tree->getOpPtrVector().push_back(
3179 new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
3183 CHKERR add_contact_three();
3199 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
3201 auto get_op_contact_bc = [&]() {
3205 return op_loop_side;
3213 boost::shared_ptr<FEMethod> null;
3215 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
3247 bool set_ts_monitor) {
3249#ifdef ENABLE_PYTHON_BINDING
3250 auto setup_sdf = [&]() {
3251 boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
3253 auto file_exists = [](std::string myfile) {
3254 std::ifstream file(myfile.c_str());
3261 char sdf_file_name[255] =
"sdf.py";
3263 sdf_file_name, 255, PETSC_NULLPTR);
3265 if (file_exists(sdf_file_name)) {
3266 MOFEM_LOG(
"EP", Sev::inform) << sdf_file_name <<
" file found";
3267 sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
3268 CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
3269 ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
3270 MOFEM_LOG(
"EP", Sev::inform) <<
"SdfPython initialized";
3272 MOFEM_LOG(
"EP", Sev::warning) << sdf_file_name <<
" file NOT found";
3274 return sdf_python_ptr;
3276 auto sdf_python_ptr = setup_sdf();
3279 auto setup_ts_monitor = [&]() {
3280 boost::shared_ptr<TsCtx>
ts_ctx;
3282 if (set_ts_monitor) {
3284 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
3288 MOFEM_LOG(
"EP", Sev::inform) <<
"TS monitor setup";
3289 return std::make_tuple(
ts_ctx);
3292 auto setup_snes_monitor = [&]() {
3295 CHKERR TSGetSNES(ts, &snes);
3297 CHKERR SNESMonitorSet(snes,
3300 (
void *)(snes_ctx.get()), PETSC_NULLPTR);
3301 MOFEM_LOG(
"EP", Sev::inform) <<
"SNES monitor setup";
3305 auto setup_snes_conergence_test = [&]() {
3308 auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
3309 PetscReal snorm, PetscReal fnorm,
3310 SNESConvergedReason *reason,
void *cctx) {
3313 CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
3317 CHKERR SNESGetSolutionUpdate(snes, &x_update);
3318 CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
3386 auto setup_section = [&]() {
3387 PetscSection section_raw;
3390 CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
3391 for (
int ff = 0; ff != num_fields; ff++) {
3399 auto set_vector_on_mesh = [&]() {
3403 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
3404 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
3405 MOFEM_LOG(
"EP", Sev::inform) <<
"Vector set on mesh";
3409 auto setup_schur_block_solver = [&]() {
3410 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver";
3411 CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
3412 CHKERR TSSetFromOptions(ts);
3415 boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
3416 if constexpr (A == AssemblyType::BLOCK_MAT) {
3421 MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver done";
3428#ifdef ENABLE_PYTHON_BINDING
3429 return std::make_tuple(setup_sdf(), setup_ts_monitor(),
3430 setup_snes_monitor(), setup_snes_conergence_test(),
3431 setup_section(), set_vector_on_mesh(),
3432 setup_schur_block_solver());
3434 return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
3435 setup_snes_conergence_test(), setup_section(),
3436 set_vector_on_mesh(), setup_schur_block_solver());
3444 PetscBool debug_model = PETSC_FALSE;
3448 if (debug_model == PETSC_TRUE) {
3450 auto post_proc = [&](TS ts, PetscReal
t, Vec u, Vec u_t, Vec u_tt, Vec
F,
3455 CHKERR TSGetSNES(ts, &snes);
3457 CHKERR SNESGetIterationNumber(snes, &it);
3458 std::string file_name =
"snes_iteration_" + std::to_string(it) +
".h5m";
3460 std::string file_skel_name =
3461 "snes_iteration_skel_" + std::to_string(it) +
".h5m";
3463 auto get_material_force_tag = [&]() {
3473 {get_material_force_tag()});
3477 ts_ctx_ptr->tsDebugHook = post_proc;
3482 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
3484 CHKERR VecDuplicate(x, &xx);
3485 CHKERR VecZeroEntries(xx);
3486 CHKERR TS2SetSolution(ts, x, xx);
3489 CHKERR TSSetSolution(ts, x);
3492 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3498 CHKERR TSSolve(ts, PETSC_NULLPTR);
3500 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3504 CHKERR TSGetSNES(ts, &snes);
3505 int lin_solver_iterations;
3506 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
3508 <<
"Number of linear solver iterations " << lin_solver_iterations;
3510 PetscBool test_cook_flg = PETSC_FALSE;
3513 if (test_cook_flg) {
3514 constexpr int expected_lin_solver_iterations = 11;
3515 if (lin_solver_iterations > expected_lin_solver_iterations)
3518 "Expected number of iterations is different than expected %d > %d",
3519 lin_solver_iterations, expected_lin_solver_iterations);
3522 PetscBool test_sslv116_flag = PETSC_FALSE;
3524 &test_sslv116_flag, PETSC_NULLPTR);
3526 if (test_sslv116_flag) {
3527 double max_val = 0.0;
3528 double min_val = 0.0;
3529 auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3531 auto ent_type = ent_ptr->getEntType();
3532 if (ent_type == MBVERTEX) {
3533 max_val = std::max(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], max_val);
3534 min_val = std::min(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], min_val);
3541 double global_max_val = 0.0;
3542 double global_min_val = 0.0;
3543 MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
3545 MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
3552 double ref_max_val = 0.00767;
3553 double ref_min_val = -0.00329;
3554 if (std::abs(global_max_val - ref_max_val) > 1e-5) {
3556 "Incorrect max value of the displacement field: %f != %f",
3557 global_max_val, ref_max_val);
3559 if (std::abs(global_min_val - ref_min_val) > 4e-5) {
3561 "Incorrect min value of the displacement field: %f != %f",
3562 global_min_val, ref_min_val);
3576 double final_time = 1;
3577 double delta_time = 0.1;
3579 PetscBool ts_h1_update = PETSC_FALSE;
3581 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Dynamic Relaxation Options",
3584 CHKERR PetscOptionsScalar(
"-dynamic_final_time",
3585 "dynamic relaxation final time",
"", final_time,
3586 &final_time, PETSC_NULLPTR);
3587 CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
3588 "dynamic relaxation final time",
"", delta_time,
3589 &delta_time, PETSC_NULLPTR);
3590 CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
3591 max_it, &max_it, PETSC_NULLPTR);
3592 CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
3593 ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
3597 auto setup_ts_monitor = [&]() {
3598 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*
this);
3601 auto monitor_ptr = setup_ts_monitor();
3603 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3608 double ts_delta_time;
3609 CHKERR TSGetTimeStep(ts, &ts_delta_time);
3621 monitor_ptr->ts_u = PETSC_NULLPTR;
3634 CHKERR TSSetStepNumber(ts, 0);
3636 CHKERR TSSetTimeStep(ts, ts_delta_time);
3637 if (!ts_h1_update) {
3640 CHKERR TSSolve(ts, PETSC_NULLPTR);
3641 if (!ts_h1_update) {
3645 monitor_ptr->ts_u = PETSC_NULLPTR;
3656 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3665 auto set_block = [&](
auto name,
int dim) {
3666 std::map<int, Range> map;
3667 auto set_tag_impl = [&](
auto name) {
3672 std::regex((boost::format(
"%s(.*)") % name).str())
3675 for (
auto bc : bcs) {
3679 map[bc->getMeshsetId()] = r;
3684 CHKERR set_tag_impl(name);
3686 return std::make_pair(name, map);
3689 auto set_skin = [&](
auto &&map) {
3690 for (
auto &
m : map.second) {
3697 auto set_tag = [&](
auto &&map) {
3699 auto name = map.first;
3700 int def_val[] = {-1};
3703 MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
3705 for (
auto &
m : map.second) {
3716 set_tag(set_skin(set_block(
"MAT_NEOHOOKEAN", 3))));
3724 Vec f_residual, Vec var_vector,
3725 std::vector<Tag> tags_to_transfer) {
3732 if (
rval == MB_SUCCESS) {
3735 int def_val[] = {0};
3737 "CRACK", 1, MB_TYPE_INTEGER,
th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
3740 tags_to_transfer.push_back(
th);
3742 auto get_tag = [&](
auto name,
auto dim) {
3745 double def_val[] = {0., 0., 0.};
3746 CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3747 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3752 tags_to_transfer.push_back(get_tag(
"MaterialForce", 3));
3758 tags_to_transfer.push_back(
t);
3766 struct exclude_sdf {
3767 exclude_sdf(
Range &&r) : map(r) {}
3768 bool operator()(
FEMethod *fe_method_ptr) {
3770 if (map.find(ent) != map.end()) {
3785 auto get_post_proc = [&](
auto &post_proc_mesh,
auto sense) {
3787 auto post_proc_ptr =
3788 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3794 auto domain_ops = [&](
auto &fe,
int sense) {
3796 fe.getUserPolynomialBase() =
3801 auto piola_scale_ptr = boost::make_shared<double>(1.0);
3810 fe.getOpPtrVector().push_back(
3814 fe.getOpPtrVector().push_back(
3822 boost::make_shared<double>(1), vec));
3825 boost::make_shared<double>(1), vec, MBMAXTYPE));
3829 fe.getOpPtrVector().push_back(
3833 fe.getOpPtrVector().push_back(
3851 fe.getOpPtrVector().push_back(
3859 fe.getOpPtrVector().push_back(op);
3866 VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
3868 struct OpSidePPMap :
public OpPPMap {
3869 OpSidePPMap(moab::Interface &post_proc_mesh,
3870 std::vector<EntityHandle> &map_gauss_pts,
3871 DataMapVec data_map_scalar, DataMapMat data_map_vec,
3872 DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
3874 :
OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
3875 data_map_vec, data_map_mat, data_symm_map_mat),
3882 if (tagSense != OpPPMap::getSkeletonSense())
3894 vec_fields[
"SpatialDisplacementL2"] =
dataAtPts->getSmallWL2AtPts();
3895 vec_fields[
"SpatialDisplacementH1"] =
dataAtPts->getSmallWH1AtPts();
3896 vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
3897 vec_fields[
"ContactDisplacement"] =
dataAtPts->getContactL2AtPts();
3898 vec_fields[
"AngularMomentum"] =
dataAtPts->getLeviKirchhoffAtPts();
3899 vec_fields[
"X"] =
dataAtPts->getLargeXH1AtPts();
3901 vec_fields[
"EiegnLogStreach"] =
dataAtPts->getEigenValsAtPts();
3905 vec_fields[
"VarOmega"] =
dataAtPts->getVarRotAxisPts();
3906 vec_fields[
"VarSpatialDisplacementL2"] =
3907 boost::make_shared<MatrixDouble>();
3909 spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], vec, MBTET));
3913 vec_fields[
"ResSpatialDisplacementL2"] =
3914 boost::make_shared<MatrixDouble>();
3916 spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], vec, MBTET));
3917 vec_fields[
"ResOmega"] = boost::make_shared<MatrixDouble>();
3919 rotAxis, vec_fields[
"ResOmega"], vec, MBTET));
3923 mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
3925 mat_fields[
"VarPiolaStress"] =
dataAtPts->getVarPiolaPts();
3929 mat_fields[
"ResPiolaStress"] = boost::make_shared<MatrixDouble>();
3932 boost::make_shared<double>(1), vec));
3935 boost::make_shared<double>(1), vec, MBMAXTYPE));
3941 fe.getOpPtrVector().push_back(
3945 fe.getOpPtrVector().push_back(
3950 "Unsupported internal stress interpolation order %d",
3956 mat_fields_symm[
"LogSpatialStretch"] =
3958 mat_fields_symm[
"SpatialStretch"] =
dataAtPts->getStretchTensorAtPts();
3960 mat_fields_symm[
"VarLogSpatialStretch"] =
3966 mat_fields_symm[
"ResLogSpatialStretch"] =
3967 boost::make_shared<MatrixDouble>();
3968 fe.getOpPtrVector().push_back(
3970 stretchTensor, mat_fields_symm[
"ResLogSpatialStretch"], vec,
3975 fe.getOpPtrVector().push_back(
3979 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3996 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4002 post_proc_ptr->getOpPtrVector().push_back(
4006 auto X_h1_ptr = boost::make_shared<MatrixDouble>();
4008 post_proc_ptr->getOpPtrVector().push_back(
4015 domain_ops(*(op_loop_side->getSideFEPtr()), sense);
4016 post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
4018 return post_proc_ptr;
4022 auto calcs_side_traction_and_displacements = [&](
auto &post_proc_ptr,
4025 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
4029 using SideEleOp = EleOnSide::UserDataOperator;
4032 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4035 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4037 op_loop_domain_side->getOpPtrVector().push_back(
4039 piolaStress, contact_common_data_ptr->contactTractionPtr(),
4040 boost::make_shared<double>(1.0)));
4041 pip.push_back(op_loop_domain_side);
4043 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
4046 contactDisp, contact_common_data_ptr->contactDispPtr()));
4050 &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
4057 pip.push_back(op_this);
4058 auto contact_residual = boost::make_shared<MatrixDouble>();
4059 op_this->getOpPtrVector().push_back(
4064 op_this->getOpPtrVector().push_back(
4068 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4072 {{
"res_contact", contact_residual}},
4086 auto post_proc_mesh = boost::make_shared<moab::Core>();
4087 auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
4088 auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
4089 CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
4090 post_proc_ptr->getOpPtrVector());
4096 CHKERR mField.get_moab().get_adjacencies(own_tets,
SPACE_DIM - 1,
true,
4097 own_faces, moab::Interface::UNION);
4099 auto get_post_negative = [&](
auto &&ents) {
4100 auto crack_faces_pos = ents;
4101 auto crack_faces_neg = crack_faces_pos;
4102 auto skin =
get_skin(mField, own_tets);
4103 auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
4104 for (
auto f : crack_on_proc_skin) {
4106 CHKERR mField.get_moab().get_adjacencies(&f, 1,
SPACE_DIM,
false, tet);
4107 tet = intersect(tet, own_tets);
4108 int side_number, sense, offset;
4109 CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
4112 crack_faces_neg.erase(f);
4114 crack_faces_pos.erase(f);
4117 return std::make_pair(crack_faces_pos, crack_faces_neg);
4120 auto get_crack_faces = [&](
auto crack_faces) {
4121 auto get_adj = [&](
auto e,
auto dim) {
4123 CHKERR mField.get_moab().get_adjacencies(e, dim,
true, adj,
4124 moab::Interface::UNION);
4127 auto tets = get_adj(crack_faces, 3);
4128 auto faces = subtract(get_adj(tets, 2), crack_faces);
4129 tets = subtract(tets, get_adj(faces, 3));
4130 return subtract(crack_faces, get_adj(tets, 2));
4133 auto [crack_faces_pos, crack_faces_neg] =
4134 get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
4136 auto only_crack_faces = [&](
FEMethod *fe_method_ptr) {
4137 auto ent = fe_method_ptr->getFEEntityHandle();
4138 if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
4144 auto only_negative_crack_faces = [&](
FEMethod *fe_method_ptr) {
4145 auto ent = fe_method_ptr->getFEEntityHandle();
4146 if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
4152 auto get_adj_front = [&]() {
4153 auto skeleton_faces = *skeletonFaces;
4155 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2,
true, adj_front,
4156 moab::Interface::UNION);
4158 adj_front = intersect(adj_front, skeleton_faces);
4159 adj_front = subtract(adj_front, *crackFaces);
4160 adj_front = intersect(own_faces, adj_front);
4164 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4165 post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
4167 auto post_proc_begin =
4171 post_proc_ptr->exeTestHook = only_crack_faces;
4172 post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
4174 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4176 post_proc_negative_sense_ptr, 0,
4177 mField.get_comm_size());
4179 constexpr bool debug =
false;
4182 auto [adj_front_pos, adj_front_neg] =
4185 auto only_front_faces_pos = [adj_front_pos](
FEMethod *fe_method_ptr) {
4186 auto ent = fe_method_ptr->getFEEntityHandle();
4187 if (adj_front_pos.find(ent) == adj_front_pos.end()) {
4193 auto only_front_faces_neg = [adj_front_neg](
FEMethod *fe_method_ptr) {
4194 auto ent = fe_method_ptr->getFEEntityHandle();
4195 if (adj_front_neg.find(ent) == adj_front_neg.end()) {
4201 post_proc_ptr->exeTestHook = only_front_faces_pos;
4203 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4204 post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
4206 post_proc_negative_sense_ptr, 0,
4207 mField.get_comm_size());
4212 CHKERR post_proc_end.writeFile(file.c_str());
4219 std::vector<Tag> tags_to_transfer) {
4224 auto post_proc_mesh = boost::make_shared<moab::Core>();
4225 auto post_proc_ptr =
4226 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
4232 auto hybrid_disp = boost::make_shared<MatrixDouble>();
4233 post_proc_ptr->getOpPtrVector().push_back(
4236 auto hybrid_res = boost::make_shared<MatrixDouble>();
4237 post_proc_ptr->getOpPtrVector().push_back(
4241 post_proc_ptr->getOpPtrVector().push_back(
4245 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4249 {{
"hybrid_disp", hybrid_disp}},
4260 auto hybrid_res = boost::make_shared<MatrixDouble>();
4261 post_proc_ptr->getOpPtrVector().push_back(
4263 hybridSpatialDisp, hybrid_res,
4266 post_proc_ptr->getOpPtrVector().push_back(
4270 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4274 {{
"res_hybrid", hybrid_res}},
4285 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4287 auto post_proc_begin =
4294 CHKERR post_proc_end.writeFile(file.c_str());
4302 constexpr bool debug =
false;
4304 auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
4305 std::vector<Tag> tags;
4306 tags.reserve(names.size());
4307 auto create_and_clean = [&]() {
4309 for (
auto n : names) {
4310 tags.push_back(
Tag());
4311 auto &tag = tags.back();
4313 auto rval = moab.tag_get_handle(
n.first.c_str(), tag);
4314 if (
rval == MB_SUCCESS) {
4315 moab.tag_delete(tag);
4317 double def_val[] = {0., 0., 0.};
4318 CHKERR moab.tag_get_handle(
n.first.c_str(),
n.second, MB_TYPE_DOUBLE,
4319 tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4327 enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
4329 auto tags = get_tags_vec({{
"MaterialForce", 3},
4331 {
"GriffithForce", 1},
4332 {
"FacePressure", 1}});
4334 auto calculate_material_forces = [&]() {
4340 auto get_face_material_force_fe = [&]() {
4342 auto fe_ptr = boost::make_shared<FaceEle>(
mField);
4343 fe_ptr->getRuleHook = [](int, int, int) {
return -1; };
4344 fe_ptr->setRuleHook =
4348 fe_ptr->getOpPtrVector().push_back(
4352 auto op_loop_domain_side =
4355 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4356 boost::make_shared<CGGUserPolynomialBase>();
4358 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4360 op_loop_domain_side->getOpPtrVector().push_back(
4363 op_loop_domain_side->getOpPtrVector().push_back(
4367 op_loop_domain_side->getOpPtrVector().push_back(
4372 op_loop_domain_side->getOpPtrVector().push_back(
4374 fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4380 auto integrate_face_material_force_fe = [&](
auto &&face_energy_fe) {
4388 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
4391 CHKERR VecGetLocalSize(
v.second, &size);
4393 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
4394 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( "
4395 << low <<
" " << high <<
" ) ";
4399 CHKERR print_loc_size(face_exchange,
"material face_exchange",
4403 mField.
get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
4410 "front_skin_faces_material_force_" +
4419 CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
4424 auto calculate_front_material_force = [&]() {
4428 auto get_conn = [&](
auto e) {
4431 "get connectivity");
4435 auto get_conn_range = [&](
auto e) {
4438 "get connectivity");
4442 auto get_adj = [&](
auto e,
auto dim) {
4449 auto get_adj_range = [&](
auto e,
auto dim) {
4452 moab::Interface::UNION),
4457 auto get_material_force = [&](
auto r,
auto th) {
4458 MatrixDouble material_forces(r.size(), 3,
false);
4462 return material_forces;
4467 auto crack_edges = get_adj_range(*
crackFaces, 1);
4468 auto front_nodes = get_conn_range(*
frontEdges);
4475 Range all_skin_faces;
4476 Range all_front_faces;
4479 auto calculate_edge_direction = [&](
auto e) {
4484 "get connectivity");
4485 std::array<double, 6> coords;
4490 &coords[0], &coords[1], &coords[2]};
4492 &coords[3], &coords[4], &coords[5]};
4495 t_dir(
i) = t_p1(
i) - t_p0(
i);
4500 auto calculate_force_through_node = [&]() {
4512 auto body_skin_conn = get_conn_range(body_skin);
4515 for (
auto n : front_nodes) {
4516 auto adj_tets = get_adj(
n, 3);
4518 auto conn = get_conn_range(adj_tets);
4519 adj_tets = get_adj_range(conn, 3);
4522 auto material_forces = get_material_force(skin_faces, tags[0]);
4526 all_skin_faces.merge(skin_faces);
4530 auto calculate_node_material_force = [&]() {
4534 for (
auto face : skin_faces) {
4537 t_face_force_tmp(
I) = t_face_T(
I);
4540 auto face_tets = intersect(get_adj(face, 3), adj_tets);
4542 if (face_tets.empty()) {
4546 if (face_tets.size() != 1) {
4548 "face_tets.size() != 1");
4551 int side_number, sense, offset;
4555 "moab side number");
4556 t_face_force_tmp(
I) *= sense;
4557 t_node_force(
I) += t_face_force_tmp(
I);
4560 return t_node_force;
4563 auto calculate_crack_area_growth_direction = [&](
auto n,
4564 auto &t_node_force) {
4567 auto boundary_node = intersect(
Range(
n,
n), body_skin_conn);
4568 if (boundary_node.size()) {
4569 auto faces = intersect(get_adj(
n, 2), body_skin);
4570 for (
auto f : faces) {
4573 f, &t_normal_face(0));
4574 t_project(
I) += t_normal_face(
I);
4576 t_project.normalize();
4579 auto adj_faces = intersect(get_adj(
n, 2), *
crackFaces);
4580 if (adj_faces.empty()) {
4582 intersect(get_adj(
n, 1), unite(*
frontEdges, body_edges));
4584 for (
auto e : adj_edges) {
4585 auto t_dir = calculate_edge_direction(e);
4590 t_area_dir(
I) = -t_node_force(
I) / t_node_force.
l2();
4591 t_area_dir(
I) *=
l / 2;
4599 if (boundary_node.size()) {
4600 t_Q(
I,
J) -= t_project(
I) * t_project(
J);
4604 auto front_edges = get_adj(
n, 1);
4606 for (
auto f : adj_faces) {
4611 std::array<double, 9> coords;
4617 coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
4618 auto n_it = std::find(conn, conn + num_nodes,
n);
4619 auto n_index = std::distance(conn, n_it);
4622 t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
4623 t_d_normal(0, n_index, 2),
4625 t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
4626 t_d_normal(1, n_index, 2),
4628 t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
4629 t_d_normal(2, n_index, 2)};
4632 t_projected_hessian(
I,
J) =
4633 t_Q(
I,
K) * (t_face_hessian(
K,
L) * t_Q(
L,
J));
4636 t_face_normal(
I) * t_projected_hessian(
I,
K) / 2.;
4642 auto t_node_force = calculate_node_material_force();
4646 &
n, 1, &t_node_force(0)),
4649 auto get_area_dir = [&]() {
4651 calculate_crack_area_growth_direction(
n, t_node_force);
4656 adj_edges.merge(intersect(get_adj_range(seed_n, 1),
4658 seed_n = get_conn_range(adj_edges);
4661 skin_adj_edges.merge(
Range(
n,
n));
4662 seed_n = subtract(seed_n, skin_adj_edges);
4663 for (
auto sn : seed_n) {
4664 auto t_area_dir_sn =
4665 calculate_crack_area_growth_direction(sn, t_node_force);
4666 t_area_dir(
I) += t_area_dir_sn(
I);
4672 auto t_area_dir = get_area_dir();
4678 auto griffith = -t_node_force(
I) * t_area_dir(
I) /
4679 (t_area_dir(
K) * t_area_dir(
K));
4686 auto ave_node_force = [&](
auto th) {
4691 auto conn = get_conn(e);
4692 auto data = get_material_force(conn,
th);
4695 for (
auto n : conn) {
4696 t_edge(
I) += t_node(
I);
4699 t_edge(
I) /= conn.size();
4702 calculate_edge_direction(e);
4719 auto ave_node_griffith_energy = [&](
auto th) {
4724 tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
4727 &e, 1, &t_edge_area_dir(0));
4728 double griffith_energy = -t_edge_force(
I) * t_edge_area_dir(
I) /
4729 (t_edge_area_dir(
K) * t_edge_area_dir(
K));
4735 CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
4736 CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
4737 CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
4742 CHKERR calculate_force_through_node();
4746 auto adj_faces = get_adj(e, 2);
4747 auto crack_face = intersect(get_adj(e, 2), *
crackFaces);
4751 all_front_faces.merge(adj_faces);
4757 &e, 1, &t_edge_force(0));
4759 calculate_edge_direction(e);
4768 for (
auto f : adj_faces) {
4772 int side_number, sense, offset;
4775 auto dot = -sense * t_cross(
I) * t_normal(
I);
4777 tags[ExhangeTags::GRIFFITHFORCE], &
f, 1, &dot),
4785 CHKERR TSGetStepNumber(ts, &ts_step);
4787 "front_edges_material_force_" +
4788 std::to_string(ts_step) +
".vtk",
4791 "front_skin_faces_material_force_" +
4792 std::to_string(ts_step) +
".vtk",
4795 "front_faces_material_force_" +
4796 std::to_string(ts_step) +
".vtk",
4805 mField.
get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
4807 mField.
get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
4814 auto print_results = [&]() {
4817 auto get_conn_range = [&](
auto e) {
4820 "get connectivity");
4824 auto get_tag_data = [&](
auto &ents,
auto tag,
auto dim) {
4825 std::vector<double> data(ents.size() * dim);
4832 auto at_nodes = [&]() {
4835 auto material_force =
4836 get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
4837 auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
4838 auto griffith_force =
4839 get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
4840 std::vector<double> coords(conn.size() * 3);
4844 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Material force at nodes";
4845 for (
size_t i = 0;
i < conn.size(); ++
i) {
4847 <<
"Node " << conn[
i] <<
" coords "
4848 << coords[
i * 3 + 0] <<
" " << coords[
i * 3 + 1] <<
" "
4849 << coords[
i * 3 + 2] <<
" material force "
4850 << material_force[
i * 3 + 0] <<
" "
4851 << material_force[
i * 3 + 1] <<
" "
4852 << material_force[
i * 3 + 2] <<
" area growth "
4853 << area_growth[
i * 3 + 0] <<
" " << area_growth[
i * 3 + 1]
4854 <<
" " << area_growth[
i * 3 + 2] <<
" griffith force "
4855 << griffith_force[
i];
4865 CHKERR calculate_material_forces();
4866 CHKERR calculate_front_material_force();
4873 bool set_orientation) {
4876 constexpr bool debug =
false;
4877 constexpr auto sev = Sev::verbose;
4882 Range body_skin_edges;
4884 moab::Interface::UNION);
4885 Range boundary_skin_verts;
4887 boundary_skin_verts,
true);
4890 Range geometry_edges_verts;
4892 geometry_edges_verts,
true);
4893 Range crack_faces_verts;
4896 Range crack_faces_edges;
4898 *
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
4899 Range crack_faces_tets;
4901 *
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
4907 moab::Interface::UNION);
4908 Range front_verts_edges;
4910 front_verts, 1,
true, front_verts_edges, moab::Interface::UNION);
4912 auto get_tags_vec = [&](
auto tag_name,
int dim) {
4913 std::vector<Tag> tags(1);
4918 auto create_and_clean = [&]() {
4921 auto rval = moab.tag_get_handle(tag_name, tags[0]);
4922 if (
rval == MB_SUCCESS) {
4923 moab.tag_delete(tags[0]);
4925 double def_val[] = {0., 0., 0.};
4926 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
4927 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4936 auto get_adj_front = [&](
bool subtract_crack) {
4939 adj_front, moab::Interface::UNION);
4941 adj_front = subtract(adj_front, *
crackFaces);
4947 auto th_front_position = get_tags_vec(
"FrontPosition", 3);
4948 auto th_max_face_energy = get_tags_vec(
"MaxFaceEnergy", 1);
4952 auto get_crack_adj_tets = [&](
auto r) {
4953 Range crack_faces_conn;
4955 Range crack_faces_conn_tets;
4957 true, crack_faces_conn_tets,
4958 moab::Interface::UNION);
4959 return crack_faces_conn_tets;
4962 auto get_layers_for_sides = [&](
auto &side) {
4963 std::vector<Range> layers;
4967 auto get_adj = [&](
auto &r,
int dim) {
4970 moab::Interface::UNION);
4974 auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
4979 Range front_faces = get_adj(front_nodes, 2);
4980 front_faces = subtract(front_faces, *
crackFaces);
4981 auto front_tets = get_tets(front_nodes);
4982 auto front_side = intersect(side, front_tets);
4983 layers.push_back(front_side);
4986 adj_faces = intersect(adj_faces, front_faces);
4987 auto adj_faces_tets = get_tets(adj_faces);
4988 adj_faces_tets = intersect(adj_faces_tets, front_tets);
4989 layers.push_back(unite(layers.back(), adj_faces_tets));
4990 if (layers.back().size() == layers[layers.size() - 2].size()) {
5001 auto layers_top = get_layers_for_sides(sides_pair.first);
5002 auto layers_bottom = get_layers_for_sides(sides_pair.second);
5014 MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
5016 for (
auto &r : layers_top) {
5017 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
5020 "layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
5025 for (
auto &r : layers_bottom) {
5026 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " << r.size();
5029 "layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
5035 auto get_cross = [&](
auto t_dir,
auto f) {
5047 auto get_sense = [&](
auto f,
auto e) {
5048 int side, sense, offset;
5051 return std::make_tuple(side, sense, offset);
5054 auto calculate_edge_direction = [&](
auto e,
auto normalize =
true) {
5058 std::array<double, 6> coords;
5061 &coords[0], &coords[1], &coords[2]};
5063 &coords[3], &coords[4], &coords[5]};
5066 t_dir(
i) = t_p1(
i) - t_p0(
i);
5072 auto evaluate_face_energy_and_set_orientation = [&](
auto front_edges,
5079 Tag th_material_force;
5089 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5090 "Unknown energy release selector");
5098 auto find_maximal_face_energy = [&](
auto front_edges,
auto front_faces,
5099 auto &edge_face_max_energy_map) {
5108 for (
auto e : front_edges) {
5110 double griffith_force;
5116 faces = subtract(intersect(faces, front_faces), body_skin);
5117 std::vector<double> face_energy(faces.size());
5119 face_energy.data());
5120 auto max_energy_it =
5121 std::max_element(face_energy.begin(), face_energy.end());
5123 max_energy_it != face_energy.end() ? *max_energy_it : 0;
5125 edge_face_max_energy_map[e] =
5126 std::make_tuple(faces[max_energy_it - face_energy.begin()],
5127 griffith_force,
static_cast<double>(0));
5129 <<
"Edge " << e <<
" griffith force " << griffith_force
5130 <<
" max face energy " << max_energy <<
" factor "
5131 << max_energy / griffith_force;
5133 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
5155 auto calculate_face_orientation = [&](
auto &edge_face_max_energy_map) {
5158 auto up_down_face = [&](
5160 auto &face_angle_map_up,
5161 auto &face_angle_map_down
5166 for (
auto &
m : edge_face_max_energy_map) {
5168 auto [max_face, energy, opt_angle] =
m.second;
5172 faces = intersect(faces, front_faces);
5176 moab::Interface::UNION);
5177 if (adj_tets.size()) {
5182 moab::Interface::UNION);
5183 if (adj_tets.size()) {
5185 Range adj_tets_faces;
5188 adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
5189 moab::Interface::UNION);
5190 adj_tets_faces = intersect(adj_tets_faces, faces);
5195 get_cross(calculate_edge_direction(e,
true), max_face);
5196 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
5197 t_cross_max(
i) *= sense_max;
5199 for (
auto t : adj_tets) {
5200 Range adj_tets_faces;
5202 &
t, 1,
SPACE_DIM - 1,
false, adj_tets_faces);
5203 adj_tets_faces = intersect(adj_tets_faces, faces);
5205 subtract(adj_tets_faces,
Range(max_face, max_face));
5207 if (adj_tets_faces.size() == 1) {
5211 auto t_cross = get_cross(calculate_edge_direction(e,
true),
5213 auto [side, sense, offset] =
5214 get_sense(adj_tets_faces[0], e);
5215 t_cross(
i) *= sense;
5216 double dot = t_cross(
i) * t_cross_max(
i);
5217 auto angle = std::acos(dot);
5221 th_face_energy, adj_tets_faces, &face_energy);
5223 auto [side_face, sense_face, offset_face] =
5224 get_sense(
t, max_face);
5226 if (sense_face > 0) {
5227 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
5231 face_angle_map_down[e] = std::make_tuple(
5232 face_energy, -angle, adj_tets_faces[0]);
5243 auto calc_optimal_angle = [&](
5245 auto &face_angle_map_up,
5246 auto &face_angle_map_down
5251 for (
auto &
m : edge_face_max_energy_map) {
5253 auto &[max_face, e0,
a0] =
m.second;
5255 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
5257 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
5258 face_angle_map_down.find(e) == face_angle_map_down.end()) {
5266 Tag th_material_force;
5271 th_material_force, &e, 1, &t_material_force(0));
5272 auto material_force_magnitude = t_material_force.
l2();
5273 if (material_force_magnitude <
5274 std::numeric_limits<double>::epsilon()) {
5279 auto t_edge_dir = calculate_edge_direction(e,
true);
5280 auto t_cross_max = get_cross(t_edge_dir, max_face);
5281 auto [side, sense, offset] = get_sense(max_face, e);
5282 t_cross_max(sense) *= sense;
5289 t_cross_max.normalize();
5292 t_material_force(
J) * t_cross_max(
K);
5293 a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
5296 <<
"Optimal angle " <<
a0 <<
" energy " << e0;
5302 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5303 "Unknown energy release selector");
5313 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5315 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5316 face_angle_map_down;
5317 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
5318 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
5322 auto th_angle = get_tags_vec(
"Angle", 1);
5324 for (
auto &
m : face_angle_map_up) {
5325 auto [e,
a, face] =
m.second;
5330 for (
auto &
m : face_angle_map_down) {
5331 auto [e,
a, face] =
m.second;
5336 Range max_energy_faces;
5337 for (
auto &
m : edge_face_max_energy_map) {
5338 auto [face, e, angle] =
m.second;
5339 max_energy_faces.insert(face);
5355 auto get_conn = [&](
auto e) {
5362 auto get_adj = [&](
auto e,
auto dim) {
5365 e, dim,
false, adj, moab::Interface::UNION),
5370 auto get_coords = [&](
auto v) {
5378 auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
5381 auto t_edge_dir = calculate_edge_direction(e,
true);
5382 auto [side, sense, offset] = get_sense(
f, e);
5383 t_edge_dir(
i) *= sense;
5384 t_edge_dir.normalize();
5385 t_edge_dir(
i) *= angle;
5390 t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
5391 return std::make_tuple(t_normal, t_rotated_normal);
5394 auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
5395 auto &t_move,
auto gamma) {
5396 auto index = adj_vertex_tets_verts.index(
v);
5398 for (
auto ii : {0, 1, 2}) {
5399 coords[3 * index + ii] += gamma * t_move(ii);
5406 auto tets_quality = [&](
auto quality,
auto &adj_vertex_tets_verts,
5407 auto &adj_vertex_tets,
auto &coords) {
5408 for (
auto t : adj_vertex_tets) {
5412 std::array<double, 12> tet_coords;
5413 for (
auto n = 0;
n != 4; ++
n) {
5414 auto index = adj_vertex_tets_verts.index(conn[
n]);
5418 for (
auto ii = 0; ii != 3; ++ii) {
5419 tet_coords[3 *
n + ii] = coords[3 * index + ii];
5423 if (!std::isnormal(q))
5425 quality = std::min(quality, q);
5431 auto calculate_free_face_node_displacement =
5432 [&](
auto &edge_face_max_energy_map) {
5434 auto get_vertex_edges = [&](
auto vertex) {
5441 vertex_edges = subtract(vertex_edges, front_verts_edges);
5443 if (boundary_skin_verts.size() &&
5444 boundary_skin_verts.find(vertex[0]) !=
5445 boundary_skin_verts.end()) {
5446 MOFEM_LOG(
"EP", sev) <<
"Boundary vertex";
5447 vertex_edges = intersect(vertex_edges, body_skin_edges);
5449 if (geometry_edges_verts.size() &&
5450 geometry_edges_verts.find(vertex[0]) !=
5451 geometry_edges_verts.end()) {
5452 MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
5453 vertex_edges = intersect(vertex_edges, geometry_edges);
5455 if (crack_faces_verts.size() &&
5456 crack_faces_verts.find(vertex[0]) !=
5457 crack_faces_verts.end()) {
5458 MOFEM_LOG(
"EP", sev) <<
"Crack face vertex";
5459 vertex_edges = intersect(vertex_edges, crack_faces_edges);
5466 return vertex_edges;
5471 using Bundle = std::vector<
5477 std::map<EntityHandle, Bundle> edge_bundle_map;
5479 for (
auto &
m : edge_face_max_energy_map) {
5481 auto edge =
m.first;
5482 auto &[max_face, energy, opt_angle] =
m.second;
5485 auto [t_normal, t_rotated_normal] =
5486 get_rotated_normal(edge, max_face, opt_angle);
5488 auto front_vertex = get_conn(
Range(
m.first,
m.first));
5489 auto adj_tets = get_adj(
Range(max_face, max_face), 3);
5490 auto adj_tets_faces = get_adj(adj_tets, 2);
5491 auto adj_front_faces = subtract(
5492 intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
5494 if (adj_front_faces.size() > 3)
5496 "adj_front_faces.size()>3");
5500 &t_material_force(0));
5501 std::vector<double> griffith_energy(adj_front_faces.size());
5503 th_face_energy, adj_front_faces, griffith_energy.data());
5506 auto set_edge_bundle = [&](
auto min_gamma) {
5507 for (
auto rotated_f : adj_front_faces) {
5509 double rotated_face_energy =
5510 griffith_energy[adj_front_faces.index(rotated_f)];
5512 auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
5514 if (vertex.size() != 1) {
5516 "Wrong number of vertex to move");
5518 auto front_vertex_edges_vertex = get_conn(
5519 intersect(get_adj(front_vertex, 1), crack_faces_edges));
5521 vertex, front_vertex_edges_vertex);
5522 if (vertex.empty()) {
5526 auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
5529 subtract(body_skin_edges, crack_faces_edges));
5532 for (;
c < 10; ++
c) {
5534 subtract(get_adj(faces, 1), seen_front_edges);
5535 if (front_edges.size() == 0) {
5538 auto front_connected_edges =
5539 intersect(front_edges, whole_front);
5540 if (front_connected_edges.size()) {
5541 seen_front_edges.merge(front_connected_edges);
5544 faces.merge(get_adj(front_edges, 2));
5551 double rotated_face_cardinality = face_cardinality(
5557 rotated_face_cardinality = std::max(rotated_face_cardinality,
5560 auto t_vertex_coords = get_coords(vertex);
5561 auto vertex_edges = get_vertex_edges(vertex);
5570 for (
auto e_used_to_move_detection : vertex_edges) {
5571 auto edge_conn = get_conn(
Range(e_used_to_move_detection,
5572 e_used_to_move_detection));
5573 edge_conn = subtract(edge_conn, vertex);
5583 t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
5587 (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5589 (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
5592 constexpr double eps =
5593 std::numeric_limits<double>::epsilon();
5594 if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
5597 t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
5599 auto check_rotated_face_directoon = [&]() {
5601 t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
5604 (t_material_force(
i) / t_material_force.
l2()) *
5606 return -dot > 0 ? true :
false;
5609 if (check_rotated_face_directoon()) {
5612 <<
"Crack edge " << edge <<
" moved face "
5614 <<
" edge: " << e_used_to_move_detection
5615 <<
" face direction/energy " << rotated_face_energy
5616 <<
" face cardinality " << rotated_face_cardinality
5617 <<
" gamma: " << gamma;
5619 auto &bundle = edge_bundle_map[edge];
5620 bundle.emplace_back(rotated_f, e_used_to_move_detection,
5621 vertex[0], t_move, 1,
5622 rotated_face_cardinality, gamma);
5629 set_edge_bundle(std::numeric_limits<double>::epsilon());
5630 if (edge_bundle_map[edge].empty()) {
5631 set_edge_bundle(-1.);
5635 return edge_bundle_map;
5638 auto get_sort_by_energy = [&](
auto &edge_face_max_energy_map) {
5639 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
5642 for (
auto &
m : edge_face_max_energy_map) {
5644 auto &[max_face, energy, opt_angle] =
m.second;
5645 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
5648 return sort_by_energy;
5651 auto set_tag = [&](
auto &&adj_edges_map,
auto &&sort_by_energy) {
5654 Tag th_face_pressure;
5656 mField.
get_moab().tag_get_handle(
"FacePressure", th_face_pressure),
5658 auto get_face_pressure = [&](
auto face) {
5667 <<
"Number of edges to check " << sort_by_energy.size();
5669 enum face_energy { POSITIVE, NEGATIVE };
5670 constexpr bool skip_negative =
true;
5672 for (
auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
5676 for (
auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
5679 auto energy = it->first;
5680 auto [max_edge, max_face, opt_angle] = it->second;
5682 auto face_pressure = get_face_pressure(max_face);
5683 if (skip_negative) {
5684 if (fe == face_energy::POSITIVE) {
5685 if (face_pressure < 0) {
5687 <<
"Skip negative face " << max_face <<
" with energy "
5688 << energy <<
" and pressure " << face_pressure;
5695 <<
"Check face " << max_face <<
" edge " << max_edge
5696 <<
" energy " << energy <<
" optimal angle " << opt_angle
5697 <<
" face pressure " << face_pressure;
5699 auto jt = adj_edges_map.find(max_edge);
5700 if (jt == adj_edges_map.end()) {
5702 <<
"Edge " << max_edge <<
" not found in adj_edges_map";
5705 auto &bundle = jt->second;
5707 auto find_max_in_bundle_impl = [&](
auto edge,
auto &bundle,
5714 double max_quality = -2;
5715 double max_quality_evaluated = -2;
5716 double min_cardinality = std::numeric_limits<double>::max();
5720 for (
auto &b : bundle) {
5721 auto &[face, move_edge, vertex, t_move, quality, cardinality,
5724 auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
5725 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
5726 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
5728 adj_vertex_tets_verts, coords.data()),
5731 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
5732 quality = tets_quality(quality, adj_vertex_tets_verts,
5733 adj_vertex_tets, coords);
5735 auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
5739 return ((edge_gamma < 0) ? (q / 2) : q) / pow(
c, 2);
5743 if (eval_quality(quality, cardinality, edge_gamma) >=
5744 max_quality_evaluated) {
5745 max_quality = quality;
5746 min_cardinality = cardinality;
5747 vertex_max = vertex;
5749 move_edge_max = move_edge;
5750 t_move_last(
i) = t_move(
i);
5751 max_quality_evaluated =
5752 eval_quality(max_quality, min_cardinality, edge_gamma);
5756 return std::make_tuple(vertex_max, face_max, t_move_last,
5757 max_quality, min_cardinality);
5760 auto find_max_in_bundle = [&](
auto edge,
auto &bundle) {
5761 auto b_org_bundle = bundle;
5762 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
5763 auto &[vertex_max, face_max, t_move_last, max_quality,
5765 if (max_quality < 0) {
5766 for (
double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
5767 bundle = b_org_bundle;
5768 r = find_max_in_bundle_impl(edge, bundle, gamma);
5769 auto &[vertex_max, face_max, t_move_last, max_quality,
5772 <<
"Back tracking: gamma " << gamma <<
" edge " << edge
5773 <<
" quality " << max_quality <<
" cardinality "
5775 if (max_quality > 0.01) {
5777 t_move_last(
I) *= gamma;
5788 auto set_tag_to_vertex_and_face = [&](
auto &&r,
auto &quality) {
5790 auto &[
v,
f, t_move, q, cardinality] = r;
5792 if ((q > 0 && std::isnormal(q)) && energy > 0) {
5795 <<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" "
5796 << max_edge <<
" move " << t_move <<
" energy " << energy
5797 <<
" quality " << q <<
" cardinality " << cardinality;
5808 double quality = -2;
5809 CHKERR set_tag_to_vertex_and_face(
5811 find_max_in_bundle(max_edge, bundle),
5817 if (quality > 0 && std::isnormal(quality) && energy > 0) {
5819 <<
"Crack face set with quality: " << quality;
5832 MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
5833 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
5834 edge_face_max_energy_map;
5835 CHKERR find_maximal_face_energy(front_edges, front_faces,
5836 edge_face_max_energy_map);
5837 CHKERR calculate_face_orientation(edge_face_max_energy_map);
5839 MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
5842 calculate_free_face_node_displacement(edge_face_max_energy_map),
5843 get_sort_by_energy(edge_face_max_energy_map)
5851 CHKERR evaluate_face_energy_and_set_orientation(
5852 *
frontEdges, get_adj_front(
false), sides_pair, th_front_position);
5870 auto get_max_moved_faces = [&]() {
5871 Range max_moved_faces;
5872 auto adj_front = get_adj_front(
false);
5873 std::vector<double> face_energy(adj_front.size());
5875 face_energy.data());
5876 for (
int i = 0;
i != adj_front.size(); ++
i) {
5877 if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
5878 max_moved_faces.insert(adj_front[
i]);
5882 return boost::make_shared<Range>(max_moved_faces);
5893 "max_moved_faces_" +
5908 Tag th_front_position;
5910 mField.
get_moab().tag_get_handle(
"FrontPosition", th_front_position);
5915 std::vector<double> coords(3 * verts.size());
5917 std::vector<double> pos(3 * verts.size());
5919 for (
int i = 0;
i != 3 * verts.size(); ++
i) {
5920 coords[
i] += pos[
i];
5923 double zero[] = {0., 0., 0.};
5928 constexpr bool debug =
false;
5933 "set_coords_faces_" +
5944 constexpr bool potential_crack_debug =
false;
5945 if constexpr (potential_crack_debug) {
5948 Range crack_front_verts;
5953 Range crack_front_faces;
5955 true, crack_front_faces,
5956 moab::Interface::UNION);
5957 crack_front_faces = intersect(crack_front_faces, add_ents);
5964 auto get_crack_faces = [&]() {
5972 auto get_extended_crack_faces = [&]() {
5973 auto get_faces_of_crack_front_verts = [&](
auto crack_faces_org) {
5974 ParallelComm *pcomm =
5979 if (!pcomm->rank()) {
5981 auto get_nodes = [&](
auto &&e) {
5984 "get connectivity");
5988 auto get_adj = [&](
auto &&e,
auto dim,
5989 auto t = moab::Interface::UNION) {
6001 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6004 auto front_block_nodes = get_nodes(front_block_edges);
6008 s = crack_faces.size();
6010 auto crack_face_nodes = get_nodes(crack_faces_org);
6011 auto crack_faces_edges =
6012 get_adj(crack_faces_org, 1, moab::Interface::UNION);
6015 front_block_edges = subtract(front_block_edges, crack_skin);
6016 auto crack_skin_nodes = get_nodes(crack_skin);
6017 crack_skin_nodes.merge(front_block_nodes);
6019 auto crack_skin_faces =
6020 get_adj(crack_skin, 2, moab::Interface::UNION);
6022 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
6024 crack_faces = crack_faces_org;
6025 for (
auto f : crack_skin_faces) {
6026 auto edges = intersect(
6027 get_adj(
Range(
f,
f), 1, moab::Interface::UNION), crack_skin);
6031 if (edges.size() == 2) {
6033 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6037 if (edges.size() == 2) {
6038 auto edge_conn = get_nodes(
Range(edges));
6039 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
6041 if (faces.size() == 2) {
6042 auto edge0_conn = get_nodes(
Range(edges[0], edges[0]));
6043 auto edge1_conn = get_nodes(
Range(edges[1], edges[1]));
6044 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
6050 moab::Interface::INTERSECT),
6055 if (node_edges.size()) {
6060 auto get_t_dir = [&](
auto e_conn) {
6061 auto other_node = subtract(e_conn,
edges_conn);
6065 t_dir(
i) -= t_v0(
i);
6071 get_t_dir(edge0_conn)(
i) + get_t_dir(edge1_conn)(
i);
6074 t_crack_surface_ave_dir(
i) = 0;
6075 for (
auto e : node_edges) {
6076 auto e_conn = get_nodes(
Range(e, e));
6077 auto t_dir = get_t_dir(e_conn);
6078 t_crack_surface_ave_dir(
i) += t_dir(
i);
6081 auto dot = t_ave_dir(
i) * t_crack_surface_ave_dir(
i);
6084 if (dot < -std::numeric_limits<double>::epsilon()) {
6085 crack_faces.insert(
f);
6088 crack_faces.insert(
f);
6092 }
else if (edges.size() == 3) {
6093 crack_faces.insert(
f);
6097 if (edges.size() == 1) {
6099 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6102 intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
6103 front_block_edges));
6104 if (edges.size() == 2) {
6105 crack_faces.insert(
f);
6111 crack_faces_org = crack_faces;
6113 }
while (s != crack_faces.size());
6119 return get_faces_of_crack_front_verts(get_crack_faces());
6124 get_extended_crack_faces());
6127 auto reconstruct_crack_faces = [&](
auto crack_faces) {
6128 ParallelComm *pcomm =
6134 Range new_crack_faces;
6135 if (!pcomm->rank()) {
6137 auto get_nodes = [&](
auto &&e) {
6140 "get connectivity");
6144 auto get_adj = [&](
auto &&e,
auto dim,
6145 auto t = moab::Interface::UNION) {
6153 auto get_test_on_crack_surface = [&]() {
6154 auto crack_faces_nodes =
6155 get_nodes(crack_faces);
6156 auto crack_faces_tets =
6157 get_adj(crack_faces_nodes, 3,
6158 moab::Interface::UNION);
6162 auto crack_faces_tets_nodes =
6163 get_nodes(crack_faces_tets);
6164 crack_faces_tets_nodes =
6165 subtract(crack_faces_tets_nodes, crack_faces_nodes);
6167 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
6168 moab::Interface::UNION));
6170 get_adj(crack_faces_tets, 2,
6171 moab::Interface::UNION);
6173 new_crack_faces.merge(crack_faces);
6175 return std::make_tuple(new_crack_faces, crack_faces_tets);
6178 auto carck_faces_test_edges = [&](
auto faces,
auto tets) {
6179 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
6180 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
6181 moab::Interface::UNION);
6182 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
6185 adj_faces_edges.merge(geometry_edges);
6186 adj_faces_edges.merge(front_block_edges);
6188 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
6189 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
6190 auto boundary_test_nodes_edges =
6191 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
6192 auto boundary_test_nodes_edges_nodes = subtract(
6193 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
6195 boundary_tets_edges =
6196 subtract(boundary_test_nodes_edges,
6197 get_adj(boundary_test_nodes_edges_nodes, 1,
6198 moab::Interface::UNION));
6205 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6206 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
6208 body_skin = intersect(body_skin, adj_tets_faces);
6209 body_skin_edges = subtract(
6210 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
6213 for (
auto e : body_skin_edges) {
6214 auto adj_tet = intersect(
6215 get_adj(
Range(e, e), 3, moab::Interface::INTERSECT), tets);
6216 if (adj_tet.size() == 1) {
6217 boundary_tets_edges.insert(e);
6221 return boundary_tets_edges;
6224 auto p = get_test_on_crack_surface();
6225 auto &[new_crack_faces, crack_faces_tets] = p;
6236 auto boundary_tets_edges =
6237 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
6239 boundary_tets_edges);
6241 auto resolve_surface = [&](
auto boundary_tets_edges,
6242 auto crack_faces_tets) {
6243 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
6244 auto crack_faces_tets_faces =
6245 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6247 Range all_removed_faces;
6248 Range all_removed_tets;
6252 while (size != crack_faces_tets.size()) {
6254 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6258 auto skin_skin_nodes = get_nodes(skin_skin);
6260 size = crack_faces_tets.size();
6262 <<
"Crack faces tets size " << crack_faces_tets.size()
6263 <<
" crack faces size " << crack_faces_tets_faces.size();
6264 auto skin_tets_nodes = subtract(
6265 get_nodes(skin_tets),
6266 boundary_tets_edges_nodes);
6268 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
6270 Range removed_nodes;
6271 Range tets_to_remove;
6272 Range faces_to_remove;
6273 for (
auto n : skin_tets_nodes) {
6275 intersect(get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT),
6277 if (tets.size() == 0) {
6281 auto hole_detetction = [&]() {
6283 get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT);
6288 if (adj_tets.size() == 0) {
6289 return std::make_pair(
6291 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6296 std::vector<Range> tets_groups;
6297 auto test_adj_tets = adj_tets;
6298 while (test_adj_tets.size()) {
6300 Range seed =
Range(test_adj_tets[0], test_adj_tets[0]);
6301 while (seed.size() != seed_size) {
6303 subtract(get_adj(seed, 2, moab::Interface::UNION),
6306 seed_size = seed.size();
6308 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
6311 tets_groups.push_back(seed);
6312 test_adj_tets = subtract(test_adj_tets, seed);
6314 if (tets_groups.size() == 1) {
6316 return std::make_pair(
6318 get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
6324 Range tets_to_remove;
6325 Range faces_to_remove;
6326 for (
auto &r : tets_groups) {
6327 auto f = get_adj(r, 2, moab::Interface::UNION);
6328 auto t = intersect(get_adj(
f, 3, moab::Interface::UNION),
6331 if (
f.size() > faces_to_remove.size() ||
6332 faces_to_remove.size() == 0) {
6333 faces_to_remove =
f;
6338 <<
"Hole detection: faces to remove "
6339 << faces_to_remove.size() <<
" tets to remove "
6340 << tets_to_remove.size();
6341 return std::make_pair(faces_to_remove, tets_to_remove);
6344 if (tets.size() < tets_to_remove.size() ||
6345 tets_to_remove.size() == 0) {
6347 auto [h_faces_to_remove, h_tets_to_remove] =
6349 faces_to_remove = h_faces_to_remove;
6350 tets_to_remove = h_tets_to_remove;
6358 all_removed_faces.merge(faces_to_remove);
6359 all_removed_tets.merge(tets_to_remove);
6362 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
6363 crack_faces_tets_faces =
6364 subtract(crack_faces_tets_faces, faces_to_remove);
6369 boost::lexical_cast<std::string>(counter) +
".vtk",
6372 "faces_to_remove_" +
6373 boost::lexical_cast<std::string>(counter) +
".vtk",
6377 boost::lexical_cast<std::string>(counter) +
".vtk",
6380 "crack_faces_tets_faces_" +
6381 boost::lexical_cast<std::string>(counter) +
".vtk",
6382 crack_faces_tets_faces);
6384 "crack_faces_tets_" +
6385 boost::lexical_cast<std::string>(counter) +
".vtk",
6391 auto cese_internal_faces = [&]() {
6394 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
6396 subtract(adj_faces, skin_tets);
6397 auto adj_tets = get_adj(adj_faces, 3,
6398 moab::Interface::UNION);
6401 subtract(crack_faces_tets,
6404 crack_faces_tets_faces =
6405 subtract(crack_faces_tets_faces, adj_faces);
6407 all_removed_faces.merge(adj_faces);
6408 all_removed_tets.merge(adj_tets);
6412 <<
"Remove internal faces size " << adj_faces.size()
6413 <<
" tets size " << adj_tets.size();
6417 auto case_only_one_free_edge = [&]() {
6420 for (
auto t :
Range(crack_faces_tets)) {
6422 auto adj_faces = get_adj(
6424 moab::Interface::UNION);
6425 auto crack_surface_edges =
6426 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6429 moab::Interface::UNION);
6432 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6433 crack_surface_edges);
6434 adj_edges = subtract(
6436 boundary_tets_edges);
6438 if (adj_edges.size() == 1) {
6440 subtract(crack_faces_tets,
6444 auto faces_to_remove =
6445 get_adj(adj_edges, 2, moab::Interface::UNION);
6448 crack_faces_tets_faces =
6449 subtract(crack_faces_tets_faces, faces_to_remove);
6451 all_removed_faces.merge(faces_to_remove);
6452 all_removed_tets.merge(
Range(
t,
t));
6454 MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Remove free one edges ";
6458 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6459 crack_faces_tets_faces =
6460 subtract(crack_faces_tets_faces, all_removed_faces);
6465 auto cese_flat_tet = [&](
auto max_adj_edges) {
6472 auto body_skin_edges =
6473 get_adj(body_skin, 1, moab::Interface::UNION);
6475 for (
auto t :
Range(crack_faces_tets)) {
6477 auto adj_faces = get_adj(
6479 moab::Interface::UNION);
6480 auto crack_surface_edges =
6481 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6484 moab::Interface::UNION);
6487 subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
6488 crack_surface_edges);
6489 adj_edges = subtract(adj_edges, body_skin_edges);
6491 auto tet_edges = get_adj(
Range(
t,
t), 1,
6492 moab::Interface::UNION);
6494 tet_edges = subtract(tet_edges, adj_edges);
6496 for (
auto e : tet_edges) {
6497 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
6498 auto get_side = [&](
auto e) {
6499 int side, sense, offset;
6502 "get side number failed");
6505 auto get_side_ent = [&](
auto side) {
6512 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
6515 if (adj_edges.size() <= max_adj_edges) {
6518 Range faces_to_remove;
6519 for (
auto e : adj_edges) {
6520 auto edge_adj_faces =
6521 get_adj(
Range(e, e), 2, moab::Interface::UNION);
6522 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
6523 if (edge_adj_faces.size() != 2) {
6525 "Adj faces size is not 2 for edge " +
6526 boost::lexical_cast<std::string>(e));
6529 auto get_normal = [&](
auto f) {
6533 "get tri normal failed");
6536 auto t_n0 = get_normal(edge_adj_faces[0]);
6537 auto t_n1 = get_normal(edge_adj_faces[1]);
6538 auto get_sense = [&](
auto f) {
6539 int side, sense, offset;
6542 "get side number failed");
6545 auto sense0 = get_sense(edge_adj_faces[0]);
6546 auto sense1 = get_sense(edge_adj_faces[1]);
6551 auto dot_e = (sense0 * sense1) * t_n0(
i) * t_n1(
i);
6552 if (dot_e < dot || e == adj_edges[0]) {
6554 faces_to_remove = edge_adj_faces;
6558 all_removed_faces.merge(faces_to_remove);
6559 all_removed_tets.merge(
Range(
t,
t));
6562 <<
"Remove free edges on flat tet, with considered nb. of "
6564 << adj_edges.size();
6568 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6569 crack_faces_tets_faces =
6570 subtract(crack_faces_tets_faces, all_removed_faces);
6576 "Case only one free edge failed");
6577 for (
auto max_adj_edges : {0, 1, 2, 3}) {
6579 "Case only one free edge failed");
6582 "Case internal faces failed");
6586 "crack_faces_tets_faces_" +
6587 boost::lexical_cast<std::string>(counter) +
".vtk",
6588 crack_faces_tets_faces);
6590 "crack_faces_tets_" +
6591 boost::lexical_cast<std::string>(counter) +
".vtk",
6595 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
6596 all_removed_faces, all_removed_tets);
6599 auto [resolved_faces, resolved_tets, all_removed_faces,
6601 resolve_surface(boundary_tets_edges, crack_faces_tets);
6602 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
6610 crack_faces = resolved_faces;
6622 auto resolve_consisten_crack_extension = [&]() {
6624 auto crack_meshset =
6627 auto meshset = crack_meshset->getMeshset();
6631 Range old_crack_faces;
6634 auto extendeded_crack_faces = get_extended_crack_faces();
6635 auto reconstructed_crack_faces =
6636 subtract(reconstruct_crack_faces(extendeded_crack_faces),
6638 if (
nbCrackFaces >= reconstructed_crack_faces.size()) {
6640 <<
"No new crack faces to add, skipping adding to meshset";
6641 extendeded_crack_faces = subtract(
6642 extendeded_crack_faces, subtract(*
crackFaces, old_crack_faces));
6644 <<
"Number crack faces size (extended) "
6645 << extendeded_crack_faces.size();
6651 reconstructed_crack_faces);
6653 <<
"Number crack faces size (reconstructed) "
6654 << reconstructed_crack_faces.size();
6673 CHKERR resolve_consisten_crack_extension();
6686 verts = subtract(verts, conn);
6687 std::vector<double> coords(3 * verts.size());
6689 double def_coords[] = {0., 0., 0.};
6692 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
6693 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
6713 auto post_proc_norm_fe =
6714 boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
6716 auto post_proc_norm_rule_hook = [](int, int,
int p) ->
int {
6719 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
6721 post_proc_norm_fe->getUserPolynomialBase() =
6728 enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
6731 CHKERR VecZeroEntries(norms_vec);
6733 auto u_l2_ptr = boost::make_shared<MatrixDouble>();
6734 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
6735 post_proc_norm_fe->getOpPtrVector().push_back(
6737 post_proc_norm_fe->getOpPtrVector().push_back(
6739 post_proc_norm_fe->getOpPtrVector().push_back(
6741 post_proc_norm_fe->getOpPtrVector().push_back(
6743 post_proc_norm_fe->getOpPtrVector().push_back(
6747 auto piola_ptr = boost::make_shared<MatrixDouble>();
6748 post_proc_norm_fe->getOpPtrVector().push_back(
6750 post_proc_norm_fe->getOpPtrVector().push_back(
6754 post_proc_norm_fe->getOpPtrVector().push_back(
6757 TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
6759 *post_proc_norm_fe);
6760 TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
6762 CHKERR VecAssemblyBegin(norms_vec);
6763 CHKERR VecAssemblyEnd(norms_vec);
6764 const double *norms;
6765 CHKERR VecGetArrayRead(norms_vec, &norms);
6766 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
6767 MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
6769 <<
"norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
6771 <<
"norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
6772 CHKERR VecRestoreArrayRead(norms_vec, &norms);
6786 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6787 if (
auto disp_bc = bc.second->dispBcPtr) {
6792 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6793 MOFEM_LOG(
"EP", Sev::noisy) <<
"Displacement BC: " << *disp_bc;
6795 std::vector<double> block_attributes(6, 0.);
6796 if (disp_bc->data.flag1 == 1) {
6797 block_attributes[0] = disp_bc->data.value1;
6798 block_attributes[3] = 1;
6800 if (disp_bc->data.flag2 == 1) {
6801 block_attributes[1] = disp_bc->data.value2;
6802 block_attributes[4] = 1;
6804 if (disp_bc->data.flag3 == 1) {
6805 block_attributes[2] = disp_bc->data.value3;
6806 block_attributes[5] = 1;
6808 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6816 boost::make_shared<NormalDisplacementBcVec>();
6817 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6818 auto block_name =
"(.*)NORMAL_DISPLACEMENT(.*)";
6819 std::regex reg_name(block_name);
6820 if (std::regex_match(bc.first, reg_name)) {
6824 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6826 block_name, bc.second->bcAttributes,
6827 bc.second->bcEnts.subset_by_dimension(2));
6832 boost::make_shared<AnalyticalDisplacementBcVec>();
6834 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6835 auto block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
6836 std::regex reg_name(block_name);
6837 if (std::regex_match(bc.first, reg_name)) {
6841 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6843 block_name, bc.second->bcAttributes,
6844 bc.second->bcEnts.subset_by_dimension(2));
6848 auto ts_displacement =
6849 boost::make_shared<DynamicRelaxationTimeScale>(
"disp_history.txt");
6852 <<
"Add time scaling displacement BC: " << bc.blockName;
6855 ts_displacement,
"disp_history",
".txt", bc.blockName);
6858 auto ts_normal_displacement =
6859 boost::make_shared<DynamicRelaxationTimeScale>(
"normal_disp_history.txt");
6862 <<
"Add time scaling normal displacement BC: " << bc.blockName;
6865 ts_normal_displacement,
"normal_disp_history",
".txt",
6881 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6882 if (
auto force_bc = bc.second->forceBcPtr) {
6887 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6888 MOFEM_LOG(
"EP", Sev::noisy) <<
"Force BC: " << *force_bc;
6890 std::vector<double> block_attributes(6, 0.);
6891 block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
6892 block_attributes[3] = 1;
6893 block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
6894 block_attributes[4] = 1;
6895 block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
6896 block_attributes[5] = 1;
6897 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6905 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6906 auto block_name =
"(.*)PRESSURE(.*)";
6907 std::regex reg_name(block_name);
6908 if (std::regex_match(bc.first, reg_name)) {
6913 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6915 block_name, bc.second->bcAttributes,
6916 bc.second->bcEnts.subset_by_dimension(2));
6921 boost::make_shared<AnalyticalTractionBcVec>();
6923 for (
auto bc : bc_mng->getBcMapByBlockName()) {
6924 auto block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
6925 std::regex reg_name(block_name);
6926 if (std::regex_match(bc.first, reg_name)) {
6930 <<
"Field name: " <<
field_name <<
" Block name: " << block_name;
6932 block_name, bc.second->bcAttributes,
6933 bc.second->bcEnts.subset_by_dimension(2));
6938 boost::make_shared<DynamicRelaxationTimeScale>(
"traction_history.txt");
6942 ts_traction,
"traction_history",
".txt", bc.blockName);
6946 boost::make_shared<DynamicRelaxationTimeScale>(
"pressure_history.txt");
6950 ts_pressure,
"pressure_history",
".txt", bc.blockName);
6959 auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
6960 const std::string block_name,
6961 const int nb_attributes) {
6964 std::regex((boost::format(
"%s(.*)") % block_name).str()))
6966 std::vector<double> block_attributes;
6967 CHKERR it->getAttributes(block_attributes);
6968 if (block_attributes.size() < nb_attributes) {
6970 "In block %s expected %d attributes, but given %ld",
6971 it->getName().c_str(), nb_attributes, block_attributes.size());
6974 auto get_block_ents = [&]() {
6980 auto Ents = get_block_ents();
6981 ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
6990 auto ts_pre_stretch =
6991 boost::make_shared<DynamicRelaxationTimeScale>(
"externalstrain_history.txt");
6994 <<
"Add time scaling external strain: " << ext_strain_block.blockName;
6997 ts_pre_stretch,
"externalstrain_history",
".txt", ext_strain_block.blockName);
7006 auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
7009 CHKERR VecGetLocalSize(
v.second, &size);
7011 CHKERR VecGetOwnershipRange(
v.second, &low, &high);
7012 MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( " << low
7013 <<
" " << high <<
" ) ";
7045 MOFEM_LOG(
"EP", Sev::inform) <<
"Calculate crack area";
7047 for (
auto f : crack_faces) {
7050 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
7052 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
7054 if (success != MPI_SUCCESS) {
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
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 nme:nme847.
@ 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
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
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 getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
#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
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
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)
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
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)
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
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[]
MoFEMErrorCode MoFEMSNESMonitorEnergy(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
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[]
FTensor::Index< 'm', 3 > m
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
boost::shared_ptr< Range > frontAdjEdges
MoFEMErrorCode createCrackSurfaceMeshset()
const std::string skeletonElement
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
MoFEM::Interface & mField
const std::string spatialL2Disp
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
SmartPetscObj< DM > dM
Coupled problem all fields.
const std::string materialH1Positions
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
static double griffithEnergy
Griffith energy.
MoFEMErrorCode calculateCrackArea(boost::shared_ptr< double > area_ptr)
const std::string elementVolumeName
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
CommInterface::EntitiesPetscVector vertexExchange
boost::shared_ptr< Range > maxMovedFaces
const std::string spatialH1Disp
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
const std::string piolaStress
MoFEMErrorCode gettingNorms()
const std::string bubbleField
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation)
static PetscBool noStretch
MoFEMErrorCode setNewFrontCoordinates()
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
boost::shared_ptr< Range > skeletonFaces
boost::shared_ptr< PhysicalEquations > physicalEquations
static boost::function< double(const double)> f
CommInterface::EntitiesPetscVector edgeExchange
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
static int nbJIntegralLevels
MoFEMErrorCode addCrackSurfaces(const bool debug=false)
MoFEMErrorCode getSpatialDispBc()
CommInterface::EntitiesPetscVector volumeExchange
MoFEMErrorCode saveOrgCoords()
static int addCrackMeshsetId
MoFEMErrorCode createExchangeVectors(Sev sev)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< Range > crackFaces
boost::shared_ptr< Range > frontVertices
static enum EnergyReleaseSelector energyReleaseSelector
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
const std::string hybridSpatialDisp
CommInterface::EntitiesPetscVector faceExchange
boost::shared_ptr< Range > frontEdges
AnalyticalDisplacementBc(std::string name, std::vector< double > attr, Range faces)
AnalyticalTractionBc(std::string name, std::vector< double > attr, Range faces)
BcDisp(std::string name, std::vector< double > attr, Range faces)
BcRot(std::string name, std::vector< double > attr, Range faces)
CGGUserPolynomialBase(boost::shared_ptr< CachePhi > cache_phi=nullptr)
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
boost::shared_ptr< CachePhi > cachePhiPtr
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
static double inv_d_f_linear(const double v)
static PetscBool noStretch
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
static int nbJIntegralLevels
int contactRefinementLevels
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
const std::string contactDisp
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
static PetscBool internalStressVoigt
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
static PetscBool crackingOn
static boost::function< double(const double)> inv_dd_f
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
const std::string piolaStress
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
static double f_log_e(const double v)
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
MoFEMErrorCode setBlockTagsOnSkin()
static double inv_f_log(const double v)
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.
static double crackingStartTime
BitRefLevel bitAdjEntMask
bit ref level for parent parent
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
MoFEMErrorCode setElasticElementOps(const int tag)
static double dd_f_log_e_quadratic(const double v)
static double inv_d_f_log(const double v)
static double d_f_linear(const double v)
boost::shared_ptr< Range > contactFaces
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
static double dd_f_log(const double v)
const std::string bubbleField
static double dynamicTime
static double inv_dd_f_log_e(const double v)
static double inv_d_f_log_e(const double v)
static enum StretchSelector stretchSelector
static std::string internalStressTagName
static PetscBool setSingularity
const std::string skeletonElement
MoFEMErrorCode solveElastic(TS ts, Vec x)
static double inv_dd_f_log(const double v)
boost::shared_ptr< Range > skeletonFaces
static double d_f_log_e_quadratic(const double v)
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
static double d_f_log_e(const double v)
static double griffithEnergy
Griffith energy.
MoFEMErrorCode getOptions()
static enum SymmetrySelector symmetrySelector
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
boost::shared_ptr< Range > crackFaces
MoFEMErrorCode setElasticElementToTs(DM dm)
boost::shared_ptr< Range > frontAdjEdges
std::vector< Tag > listTagsToTransfer
list of tags to transfer to postprocessor
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
static enum RotSelector rotSelector
SmartPetscObj< DM > dM
Coupled problem all fields.
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.
static PetscBool l2UserBaseScale
static double d_f_log(const double v)
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x)
SmartPetscObj< DM > dmElastic
Elastic problem.
const std::string hybridSpatialDisp
const std::string skinElement
static enum EnergyReleaseSelector energyReleaseSelector
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
BitRefLevel bitAdjParent
bit ref level for parent
static double exponentBase
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
boost::shared_ptr< PhysicalEquations > physicalEquations
const std::string spatialL2Disp
const std::string spatialH1Disp
static enum RotSelector gradApproximator
static double inv_dd_f_linear(const double v)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
static double f_log_e_quadratic(const double v)
static double dd_f_linear(const double v)
MoFEMErrorCode gettingNorms()
[Getting norms]
boost::shared_ptr< Range > frontEdges
BitRefLevel bitAdjEnt
bit ref level for parent
static double f_linear(const double v)
static boost::function< double(const double)> f
static double dd_f_log_e(const double v)
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
SmartPetscObj< Vec > solTSStep
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
static double f_log(const double v)
const std::string naturalBcElement
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
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={})
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
static boost::function< double(const double)> dd_f
const std::string stretchTensor
const std::string materialH1Positions
boost::shared_ptr< Range > frontVertices
const std::string rotAxis
static boost::function< double(const double)> d_f
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
static boost::function< double(const double)> inv_d_f
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)
const std::string elementVolumeName
static PetscBool dynamicRelaxation
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
EshelbianCore(MoFEM::Interface &m_field)
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
static double inv_f_log_e(const double v)
static double inv_f_linear(const double v)
static boost::function< double(const double)> inv_f
const std::string contactElement
static int internalStressInterpOrder
BitRefLevel bitAdjParentMask
bit ref level for parent parent
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
MoFEM::Interface & mField
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)
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Apply rotation boundary condition.
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
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode postStepFun(TS ts)
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode preStepFun(TS ts)
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.
Simple interface for fast problem set-up.
static std::pair< std::string, std::string > extractStringFromBlockId(const std::string block_id, const std::string prb_name)
Extract block name and block name form block id.
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 block DOFs.
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.
Class used to pass element data to calculate base functions on tet,triangle,edge.
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Provide data structure for (tensor) field approximation.
Definition of the force bc data structure.
structure to get information form 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.
Get values at integration pts for tensor field rank 1, i.e. vector field.
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
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
Calculate base functions on tetrahedral.
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.
BoundaryEle::UserDataOperator BdyEleOp