Eshelbian plasticity implementation.
#define SINGULARITY
#include <boost/math/constants/constants.hpp>
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#pragma message "With ENABLE_PYTHON_BINDING"
#else
#pragma message "Without ENABLE_PYTHON_BINDING"
#endif
extern "C" {
}
#include <queue>
struct VolUserDataOperatorStabAssembly
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
: public FaceElementForcesAndSourcesCore::UserDataOperator {
using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
}
const EntityType type) {
ParallelComm *pcomm =
auto dim = CN::Dimension(type);
std::vector<int> sendcounts(pcomm->size());
std::vector<int> displs(pcomm->size());
std::vector<int> sendbuf(
r.size());
if (pcomm->rank() == 0) {
for (auto p = 1; p != pcomm->size(); p++) {
->getPartEntities(m_field.
get_moab(), p)
CHKERR m_field.
get_moab().get_adjacencies(part_ents, dim,
true, faces,
moab::Interface::UNION);
faces = intersect(faces, r);
sendcounts[p] = faces.size();
displs[p] = sendbuf.size();
for (auto f : faces) {
sendbuf.push_back(id);
}
}
}
int recv_data;
MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
pcomm->comm());
std::vector<int> recvbuf(recv_data);
MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
if (pcomm->rank() > 0) {
for (auto &f : recvbuf) {
}
}
}
const std::string block_name, int dim) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
faces, true),
"get meshset ents");
}
};
const std::string block_name, int dim) {
std::map<std::string, Range>
r;
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
faces, true),
"get meshset ents");
r[bc->getName()] = faces;
}
}
const unsigned int cubit_bc_type) {
CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
return meshset;
};
static auto save_range(moab::Interface &moab,
const std::string name,
const Range r, std::vector<Tag> tags = {}) {
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1,
tags.data(), tags.size());
} else {
MOFEM_LOG(
"SELF", Sev::warning) <<
"Empty range for " << name;
}
};
ParallelComm *pcomm =
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents),
"filter_pstatus");
return boundary_ents;
};
ParallelComm *pcomm =
CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
&owner_ents),
"filter_pstatus");
return owner_ents;
};
CHK_MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents),
"find_skin");
return skin_ents;
};
ParallelComm *pcomm =
Range crack_skin_without_bdy;
if (pcomm->rank() == 0) {
CHKERR moab.get_adjacencies(crack_faces, 1,
true, crack_edges,
moab::Interface::UNION);
auto crack_skin = get_skin(m_field, crack_faces);
"get_entities_by_dimension");
auto body_skin = get_skin(m_field, body_ents);
CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1,
true, body_skin_edges,
moab::Interface::UNION),
"get_adjacencies");
crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
for (
auto &
m : front_edges_map) {
auto add_front = subtract(
m.second, crack_edges);
auto i = intersect(
m.second, crack_edges);
crack_skin_without_bdy.merge(add_front);
} else {
auto i_skin = get_skin(m_field,
i);
CHKERR moab.get_adjacencies(i_skin, 1,
true, adj_i_skin,
moab::Interface::UNION);
adj_i_skin = subtract(intersect(adj_i_skin,
m.second), crack_edges);
crack_skin_without_bdy.merge(adj_i_skin);
}
}
}
return send_type(m_field, crack_skin_without_bdy, MBEDGE);
}
ParallelComm *pcomm =
MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface";
if (!pcomm->rank()) {
auto impl = [&](auto &saids) {
auto get_adj = [&](auto e, auto dim) {
e, dim, true, adj, moab::Interface::UNION),
"get adj");
return adj;
};
auto get_conn = [&](auto e) {
"get connectivity");
return conn;
};
constexpr bool debug =
false;
body_ents);
auto body_skin = get_skin(m_field, body_ents);
auto body_skin_edges = get_adj(body_skin, 1);
auto crack_skin =
subtract(get_skin(m_field, crack_faces), body_skin_edges);
auto crack_skin_conn = get_conn(crack_skin);
auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
auto crack_edges = get_adj(crack_faces, 1);
crack_edges = subtract(crack_edges, crack_skin);
auto all_tets = get_adj(crack_edges, 3);
crack_edges = subtract(crack_edges, crack_skin_conn_edges);
auto crack_conn = get_conn(crack_edges);
all_tets.merge(get_adj(crack_conn, 3));
crack_edges);
}
if (crack_faces.size()) {
auto grow = [&](
auto r) {
auto crack_faces_conn = get_conn(crack_faces);
auto size_r = 0;
while (size_r !=
r.size() &&
r.size() > 0) {
CHKERR moab.get_connectivity(r,
v,
true);
v = subtract(
v, crack_faces_conn);
moab::Interface::UNION);
r = intersect(r, all_tets);
}
break;
}
}
};
Range all_tets_ord = all_tets;
while (all_tets.size()) {
Range faces = get_adj(unite(saids.first, saids.second), 2);
faces = subtract(crack_faces, faces);
if (faces.size()) {
auto fit = faces.begin();
for (; fit != faces.end(); ++fit) {
tets = intersect(get_adj(
Range(*fit, *fit), 3), all_tets);
if (tets.size() == 2) {
break;
}
}
if (tets.empty()) {
break;
} else {
saids.first.insert(tets[0]);
saids.first = grow(saids.first);
all_tets = subtract(all_tets, saids.first);
if (tets.size() == 2) {
saids.second.insert(tets[1]);
saids.second = grow(saids.second);
all_tets = subtract(all_tets, saids.second);
}
}
} else {
break;
}
}
saids.first = subtract(all_tets_ord, saids.second);
saids.second = subtract(all_tets_ord, saids.first);
}
};
std::pair<Range, Range> saids;
if (crack_faces.size())
return saids;
}
MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface <- done";
return std::pair<Range, Range>();
}
struct SetIntegrationAtFrontVolume {
boost::shared_ptr<Range> front_edges)
int order_col, int order_data) {
constexpr bool debug =
false;
constexpr int numNodes = 4;
constexpr int numEdges = 6;
constexpr int refinementLevels = 4;
auto &m_field = fe_raw_ptr->
mField;
auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
auto fe_handle = fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
int rule = 2 * order_data + 1;
"wrong dimension");
}
}
auto &gauss_pts = fe_ptr->gaussPts;
gauss_pts.resize(4, nb_gauss_pts, false);
&gauss_pts(0, 0), 1);
&gauss_pts(1, 0), 1);
&gauss_pts(2, 0), 1);
&gauss_pts(3, 0), 1);
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
false);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
1);
} else {
}
};
auto get_singular_nodes = [&]() {
int num_nodes;
CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
true);
std::bitset<numNodes> singular_nodes;
for (auto nn = 0; nn != numNodes; ++nn) {
singular_nodes.set(nn);
} else {
singular_nodes.reset(nn);
}
}
return singular_nodes;
};
auto get_singular_edges = [&]() {
std::bitset<numEdges> singular_edges;
for (int ee = 0; ee != numEdges; ee++) {
singular_edges.set(ee);
} else {
singular_edges.reset(ee);
}
}
return singular_edges;
};
auto set_gauss_pts = [&](auto &ref_gauss_pts) {
fe_ptr->gaussPts.swap(ref_gauss_pts);
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
&fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
nb_gauss_pts);
};
auto singular_nodes = get_singular_nodes();
if (singular_nodes.count()) {
auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
CHKERR set_gauss_pts(it_map_ref_coords->second);
} else {
auto refine_quadrature = [&]() {
const int max_level = refinementLevels;
moab::Core moab_ref;
double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
for (int nn = 0; nn != 4; nn++)
CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
{
tets, 1, true, edges, moab::Interface::UNION);
tets, BitRefLevel().set(0),
false,
VERBOSE);
}
for (int nn = 0; nn != numNodes; nn++) {
if (singular_nodes[nn]) {
CHKERR moab_ref.side_element(tet, 0, nn, ent);
nodes_at_front.insert(ent);
}
}
auto singular_edges = get_singular_edges();
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
for (int ee = 0; ee != numEdges; ee++) {
if (singular_edges[ee]) {
CHKERR moab_ref.side_element(tet, 1, ee, ent);
CHKERR moab_ref.add_entities(meshset, &ent, 1);
}
}
for (int ll = 0; ll != max_level; ll++) {
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
BitRefLevel().set(), MBEDGE,
edges);
CHKERR moab_ref.get_adjacencies(
nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
ref_edges = intersect(ref_edges, edges);
CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
ref_edges = intersect(ref_edges, ents);
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
ref_edges, BitRefLevel().set(ll + 1));
CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
->updateMeshsetByEntitiesChildren(meshset,
BitRefLevel().set(ll + 1),
meshset, MBEDGE, true);
}
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
BitRefLevel().set(), MBTET,
tets);
}
MatrixDouble ref_coords(tets.size(), 12, false);
int tt = 0;
for (Range::iterator tit = tets.begin(); tit != tets.end();
tit++, tt++) {
int num_nodes;
CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
}
auto &data = fe_ptr->dataOnElement[
H1];
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
MatrixDouble &shape_n =
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
int gg = 0;
for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
double *tet_coords = &ref_coords(tt, 0);
double det = Tools::tetVolume(tet_coords);
det *= 6;
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (
int dd = 0;
dd != 3;
dd++) {
ref_gauss_pts(dd, gg) =
shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
}
ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
}
}
mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
};
}
}
}
}
private:
using ForcesAndSourcesCore::dataOnElement;
private:
using ForcesAndSourcesCore::ForcesAndSourcesCore;
};
static inline std::map<long int, MatrixDouble>
mapRefCoords;
};
struct SetIntegrationAtFrontFace {
using FunRule = boost::function<int(
int)>;
boost::shared_ptr<Range> front_edges)
boost::shared_ptr<Range> front_edges,
int order_col, int order_data) {
constexpr bool debug =
false;
constexpr int numNodes = 3;
constexpr int numEdges = 3;
constexpr int refinementLevels = 4;
auto &m_field = fe_raw_ptr->
mField;
auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
auto fe_handle = fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
}
}
fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
&fe_ptr->gaussPts(0, 0), 1);
&fe_ptr->gaussPts(1, 0), 1);
&fe_ptr->gaussPts(2, 0), 1);
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
false);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
1);
data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
std::copy(
Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
} else {
}
};
auto get_singular_nodes = [&]() {
int num_nodes;
CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
true);
std::bitset<numNodes> singular_nodes;
for (auto nn = 0; nn != numNodes; ++nn) {
singular_nodes.set(nn);
} else {
singular_nodes.reset(nn);
}
}
return singular_nodes;
};
auto get_singular_edges = [&]() {
std::bitset<numEdges> singular_edges;
for (int ee = 0; ee != numEdges; ee++) {
singular_edges.set(ee);
} else {
singular_edges.reset(ee);
}
}
return singular_edges;
};
auto set_gauss_pts = [&](auto &ref_gauss_pts) {
fe_ptr->gaussPts.swap(ref_gauss_pts);
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
&fe_ptr->gaussPts(1, 0), nb_gauss_pts);
};
auto singular_nodes = get_singular_nodes();
if (singular_nodes.count()) {
auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
CHKERR set_gauss_pts(it_map_ref_coords->second);
} else {
auto refine_quadrature = [&]() {
const int max_level = refinementLevels;
moab::Core moab_ref;
double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
for (int nn = 0; nn != numNodes; nn++)
CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
{
tris, 1, true, edges, moab::Interface::UNION);
tris, BitRefLevel().set(0),
false,
VERBOSE);
}
for (int nn = 0; nn != numNodes; nn++) {
if (singular_nodes[nn]) {
CHKERR moab_ref.side_element(tri, 0, nn, ent);
nodes_at_front.insert(ent);
}
}
auto singular_edges = get_singular_edges();
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
for (int ee = 0; ee != numEdges; ee++) {
if (singular_edges[ee]) {
CHKERR moab_ref.side_element(tri, 1, ee, ent);
CHKERR moab_ref.add_entities(meshset, &ent, 1);
}
}
for (int ll = 0; ll != max_level; ll++) {
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
BitRefLevel().set(), MBEDGE,
edges);
CHKERR moab_ref.get_adjacencies(
nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
ref_edges = intersect(ref_edges, edges);
CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
ref_edges = intersect(ref_edges, ents);
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, tris);
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
ref_edges, BitRefLevel().set(ll + 1));
CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
->updateMeshsetByEntitiesChildren(meshset,
BitRefLevel().set(ll + 1),
meshset, MBEDGE, true);
}
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
BitRefLevel().set(), MBTRI,
tris);
}
MatrixDouble ref_coords(tris.size(), 9, false);
int tt = 0;
for (Range::iterator tit = tris.begin(); tit != tris.end();
tit++, tt++) {
int num_nodes;
CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
}
auto &data = fe_ptr->dataOnElement[
H1];
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
MatrixDouble &shape_n =
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
int gg = 0;
for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
double *tri_coords = &ref_coords(tt, 0);
CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
auto det = t_normal.
l2();
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (
int dd = 0;
dd != 2;
dd++) {
ref_gauss_pts(dd, gg) =
shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
}
ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
}
}
mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
};
}
}
}
}
private:
using ForcesAndSourcesCore::dataOnElement;
private:
using ForcesAndSourcesCore::ForcesAndSourcesCore;
};
static inline std::map<long int, MatrixDouble>
mapRefCoords;
};
return 0;
}
}
}
const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
const char *list_symm[] = {"symm", "not_symm"};
const char *list_release[] = {"griffith_force", "griffith_skelton"};
const char *list_stretches[] = {"linear", "log", "log_quadratic"};
PetscInt choice_stretch = StretchSelector::LOG;
char analytical_expr_file_name[255] = "analytical_expr.py";
PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
CHKERR PetscOptionsScalar(
"-viscosity_alpha_omega",
"rot viscosity",
"",
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
PETSC_NULLPTR);
CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
list_rots, NO_H1_CONFIGURATION + 1,
list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
CHKERR PetscOptionsEList(
"-symm",
"symmetric variant",
"", list_symm, 2,
list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
CHKERR PetscOptionsEList(
"-stretches",
"stretches",
"", list_stretches,
StretchSelector::STRETCH_SELECTOR_LAST,
list_stretches[choice_stretch], &choice_stretch,
PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
PETSC_NULLPTR);
CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
CHKERR PetscOptionsEList(
"-energy_release_variant",
"energy release variant",
"", list_release, 2, list_release[choice_release],
&choice_release, PETSC_NULLPTR);
CHKERR PetscOptionsInt(
"-nb_J_integral_levels",
"Number of J integarl levels",
char tag_name[255] = "";
CHKERR PetscOptionsString(
"-internal_stress_tag_name",
"internal stress tag name", "", "", tag_name, 255,
PETSC_NULLPTR);
"-internal_stress_interp_order", "internal stress interpolation order",
CHKERR PetscOptionsBool(
"-internal_stress_voigt",
"Voigt index notation",
"",
PETSC_NULLPTR);
analytical_expr_file_name, 255, PETSC_NULLPTR);
PetscOptionsEnd();
"Unsupported internal stress interpolation order %d",
}
}
case StretchSelector::LINEAR:
break;
case StretchSelector::LOG:
std::numeric_limits<float>::epsilon()) {
} else {
}
break;
case StretchSelector::LOG_QUADRATIC:
"No logarithmic quadratic stretch for this case");
return 0;
};
break;
default:
break;
};
<<
"alphaOmega: -viscosity_alpha_omega " <<
alphaOmega;
MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
else
MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
<< "Stretch: -stretches " << list_stretches[choice_stretch];
? "yes"
: "no";
? "yes"
: "no";
? "yes"
: "no";
? "yes"
: "no";
? "yes"
: "no";
<< "Energy release variant: -energy_release_variant "
<< "Number of J integral levels: -nb_J_integral_levels "
#ifdef ENABLE_PYTHON_BINDING
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
if (file_exists(analytical_expr_file_name)) {
MOFEM_LOG(
"EP", Sev::inform) << analytical_expr_file_name <<
" file found";
analytical_expr_file_name);
} else {
<< analytical_expr_file_name << " file NOT found";
}
#endif
}
auto get_tets = [&]() {
return tets;
};
auto get_tets_skin = [&]() {
CHKERR skin.find_skin(0, get_tets(),
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
return tets_skin;
};
auto subtract_boundary_conditions = [&](auto &&tets_skin) {
tets_skin = subtract(tets_skin,
v.faces);
}
tets_skin = subtract(tets_skin,
v.faces);
}
tets_skin = subtract(tets_skin,
v.faces);
}
return tets_skin;
};
auto add_blockset = [&](auto block_name, auto &&tets_skin) {
auto crack_faces =
tets_skin.merge(crack_faces);
return tets_skin;
};
auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
auto contact_range =
tets_skin = subtract(tets_skin, contact_range);
return tets_skin;
};
auto get_stress_trace_faces = [&](auto &&tets_skin) {
faces, moab::Interface::UNION);
Range trace_faces = subtract(faces, tets_skin);
return trace_faces;
};
auto tets = get_tets();
auto trace_faces = get_stress_trace_faces(
subtract_blockset("CONTACT",
subtract_boundary_conditions(get_tets_skin()))
);
boost::make_shared<Range>(subtract(trace_faces, *
contactFaces));
#ifndef NDEBUG
"contact_faces_" +
"skeleton_faces_" +
#endif
auto add_broken_hdiv_field = [
this, meshset](
const std::string
field_name,
auto get_side_map_hdiv = [&]() {
return std::vector<
std::pair<EntityType,
>>{
{MBTET,
dofs_side_map);
}}
};
};
get_side_map_hdiv(), MB_TAG_DENSE,
MF_ZERO);
};
auto add_l2_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_h1_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_l2_field_by_range = [
this](
const std::string
field_name,
const int order,
const int dim,
const int field_dim,
Range &&
r) {
};
auto add_bubble_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
auto field_order_table =
const_cast<Field *
>(field_ptr)->getFieldOrderTable();
auto get_cgg_bubble_order_zero = [](int p) { return 0; };
auto get_cgg_bubble_order_tet = [](int p) {
};
field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
field_order_table[MBTRI] = get_cgg_bubble_order_zero;
field_order_table[MBTET] = get_cgg_bubble_order_tet;
};
auto add_user_l2_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
auto field_order_table =
const_cast<Field *
>(field_ptr)->getFieldOrderTable();
auto zero_dofs = [](int p) { return 0; };
field_order_table[MBVERTEX] = zero_dofs;
field_order_table[MBEDGE] = zero_dofs;
field_order_table[MBTRI] = zero_dofs;
field_order_table[MBTET] = dof_l2_tet;
};
auto get_hybridised_disp = [&]() {
auto skin = subtract_boundary_conditions(get_tets_skin());
faces.merge(intersect(bc.faces, skin));
}
return faces;
};
get_hybridised_disp());
}
auto project_ho_geometry = [&](auto field) {
};
auto get_adj_front_edges = [&](auto &front_edges) {
Range crack_front_edges_with_both_nodes_not_at_front;
CHKERR moab.get_connectivity(front_edges, front_crack_nodes,
true);
crack_front_edges, moab::Interface::UNION);
crack_front_edges =
intersect(subtract(crack_front_edges, front_edges), meshset_ents);
Range crack_front_edges_nodes;
CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
true);
crack_front_edges_nodes =
subtract(crack_front_edges_nodes, front_crack_nodes);
crack_front_edges_nodes,
SPACE_DIM - 2,
false,
crack_front_edges_with_both_nodes_not_at_front,
moab::Interface::UNION);
crack_front_edges_with_both_nodes_not_at_front = intersect(
crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
crack_front_edges_with_both_nodes_not_at_front = intersect(
crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
}
crack_front_edges_with_both_nodes_not_at_front =
send_type(
mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
boost::make_shared<Range>(
crack_front_edges_with_both_nodes_not_at_front));
};
auto [front_vertices, front_adj_edges] = get_adj_front_edges(*
frontEdges);
#ifndef NDEBUG
(boost::format("crack_faces_%d.vtk") % rank).str(),
(boost::format("front_edges_%d.vtk") % rank).str(),
}
#endif
auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
double beta = 0;
PETSC_NULLPTR);
MOFEM_LOG(
"EP", Sev::inform) <<
"Singularity eps " << beta;
for (auto edge : front_adj_edges) {
int num_nodes;
CHKERR moab.get_connectivity(edge, conn, num_nodes,
false);
double coords[6];
CHKERR moab.get_coords(conn, num_nodes, coords);
const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
coords[5] - coords[2]};
double dof[3] = {0, 0, 0};
if (front_vertices.find(conn[0]) != front_vertices.end()) {
for (
int dd = 0;
dd != 3;
dd++) {
}
} else if (front_vertices.find(conn[1]) != front_vertices.end()) {
for (
int dd = 0;
dd != 3;
dd++) {
}
}
const int idx = dit->get()->getEntDofIdx();
if (idx > 2) {
dit->get()->getFieldData() = 0;
} else {
dit->get()->getFieldData() = dof[idx];
}
}
}
};
}
auto add_field_to_fe = [this](const std::string fe,
};
}
}
auto set_fe_adjacency = [&](auto fe_name) {
boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
};
auto add_field_to_fe = [this](const std::string fe,
};
Range natural_bc_elements;
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
}
auto get_skin = [&](auto &body_ents) {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
body_ents);
}
}
}
}
}
BitRefLevel().set());
auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
for (int d : {0, 1, 2}) {
std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
->getSideDofsOnBrokenSpaceEntities(
dofs_to_remove);
dofs_to_remove);
}
};
CHKERR remove_dofs_on_broken_skin(
"ESHELBY_PLASTICITY");
}
auto set_zero_block = [&]() {
}
}
};
auto set_section = [&]() {
PetscSection section;
§ion);
CHKERR PetscSectionDestroy(§ion);
};
}
BcDisp::BcDisp(std::string name, std::vector<double> attr,
Range faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
<< "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " << faces.size();
}
BcRot::BcRot(std::string name, std::vector<double> attr,
Range faces)
: blockName(name), faces(faces) {
vals.resize(attr.size(), false);
for (int ii = 0; ii != attr.size(); ++ii) {
vals[ii] = attr[ii];
}
theta = attr[3];
}
TractionBc::TractionBc(std::string name, std::vector<double> attr,
Range faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
<< "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " << faces.size();
}
NormalDisplacementBc::NormalDisplacementBc(std::string name,
std::vector<double> attr,
: blockName(name), faces(faces) {
blockName = name;
if (attr.size() < 1) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Wrong size of normal displacement BC");
}
val = attr[0];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc " << name;
MOFEM_LOG(
"EP", Sev::inform) <<
"Add NormalDisplacementBc val " << val;
<< "Add NormalDisplacementBc nb. of faces " << faces.size();
}
PressureBc::PressureBc(std::string name, std::vector<double> attr,
Range faces)
: blockName(name), faces(faces) {
blockName = name;
if (attr.size() < 1) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Wrong size of normal displacement BC");
}
val = attr[0];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc " << name;
MOFEM_LOG(
"EP", Sev::inform) <<
"Add PressureBc val " << val;
<< "Add PressureBc nb. of faces " << faces.size();
}
ExternalStrain::ExternalStrain(std::string name, std::vector<double> attr,
: blockName(name), ents(ents) {
blockName = name;
if (attr.size() < 2) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Wrong size of external strain attribute");
}
val = attr[0];
bulkModulusK = attr[1];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain " << name;
MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain val " << val;
MOFEM_LOG(
"EP", Sev::inform) <<
"Add ExternalStrain bulk modulus K "
<< bulkModulusK;
<< "Add ExternalStrain nb. of tets " << ents.size();
}
AnalyticalDisplacementBc::AnalyticalDisplacementBc(std::string name,
std::vector<double> attr,
: blockName(name), faces(faces) {
blockName = name;
if (attr.size() < 3) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Wrong size of analytical displacement BC");
}
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
flags[ii] = attr[ii];
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalDisplacementBc " << name;
<< "Add AnalyticalDisplacementBc flags " << flags[0] << " " << flags[1]
<< " " << flags[2];
<< "Add AnalyticalDisplacementBc nb. of faces " << faces.size();
}
AnalyticalTractionBc::AnalyticalTractionBc(std::string name,
std::vector<double> attr,
: blockName(name), faces(faces) {
blockName = name;
if (attr.size() < 3) {
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
"Wrong size of analytical traction BC");
}
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
flags[ii] = attr[ii];
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add AnalyticalTractionBc " << name;
<< "Add AnalyticalTractionBc flags " << flags[0] << " " << flags[1]
<< " " << flags[2];
<< "Add AnalyticalTractionBc nb. of faces " << faces.size();
}
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
bc_ptr->resize(3);
for (
int dd = 0;
dd != 3; ++
dd)
(*bc_ptr)[
dd] = tets_skin;
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
true);
(*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
}
}
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data + 1;
}
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * (p_data + 1);
}
};
CGGUserPolynomialBase::CGGUserPolynomialBase(
boost::shared_ptr<CachePhi> cache_phi_otr)
boost::typeindex::type_index type_index,
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
CGGUserPolynomialBase::getValue(MatrixDouble &pts,
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
const auto base = this->cTx->bAse;
switch (this->cTx->sPace) {
CHKERR getValueHdivForCGGBubble(pts);
break;
data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
false);
&pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
CHKERR getValueL2AinsworthBase(pts);
break;
default:
}
}
CGGUserPolynomialBase::getValueHdivForCGGBubble(MatrixDouble &pts) {
"Wrong base, should be USER_BASE");
}
const int nb_gauss_pts = pts.size2();
auto check_cache = [
this](
int order,
int nb_gauss_pts) ->
bool {
if (cachePhiPtr) {
return cachePhiPtr->get<0>() ==
order &&
cachePhiPtr->get<1>() == nb_gauss_pts;
} else {
return false;
}
};
if (check_cache(
order, nb_gauss_pts)) {
auto &mat = cachePhiPtr->get<2>();
phi.resize(mat.size1(), mat.size2(),
false);
} else {
shapeFun.resize(nb_gauss_pts, 4, false);
&pts(2, 0), nb_gauss_pts);
double diff_shape_fun[12];
phi.resize(nb_gauss_pts, 9 * nb_base_functions,
false);
nb_gauss_pts);
if (cachePhiPtr) {
cachePhiPtr->get<0>() =
order;
cachePhiPtr->get<1>() = nb_gauss_pts;
cachePhiPtr->get<2>().resize(
phi.size1(),
phi.size2(),
false);
noalias(cachePhiPtr->get<2>()) =
phi;
}
}
}
const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
auto bubble_cache =
boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
fe->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
fe->getRuleHook = [](int, int, int) { return -1; };
}
fe->getOpPtrVector().push_back(
} else {
fe->getOpPtrVector().push_back(
}
if (calc_rates) {
} else {
fe->getOpPtrVector().push_back(
fe->getOpPtrVector().push_back(
MBTET));
}
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
}
}
var_vec));
var_vec, MBMAXTYPE));
fe->getOpPtrVector().push_back(
} else {
fe->getOpPtrVector().push_back(
}
}
} else {
}
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto rule_contact = [](int, int, int o) { return -1; };
auto set_rule_contact = [refine](
int order_col, int order_data
) {
auto rule = 2 * order_data;
fe_raw_ptr->
gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
};
fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe_rhs);
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
} else {
case 0:
fe_rhs->getOpPtrVector().push_back(
break;
case 1:
fe_rhs->getOpPtrVector().push_back(
break;
default:
"Unsupported internal stress interpolation order %d",
}
fe_rhs->getOpPtrVector().push_back(
} else {
fe_rhs->getOpPtrVector().push_back(
}
}
fe_rhs->getOpPtrVector().push_back(op);
"OpSpatialPhysicalExternalStrain not implemented for this "
"material");
}
fe_rhs->getOpPtrVector().push_back(
}
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
auto set_hybridisation = [&](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
return -1;
};
op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
flux_mat_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_ptr));
op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact = [&](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
};
auto add_ops_contact_rhs = [&](auto &pip) {
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(new OpTreeSearch(
nullptr));
contact_sfd_map_range_ptr));
pip.push_back(
};
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
pip.push_back(op_loop_skeleton_side);
};
CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
CHKERR set_contact(fe_rhs->getOpPtrVector());
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<
GAUSS>;
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
auto body_time_scale =
boost::make_shared<DynamicRelaxationTimeScale>("body_force.txt");
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
"BODY_FORCE", Sev::inform);
}
fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe_lhs);
if (add_elastic) {
fe_lhs->getOpPtrVector().push_back(
fe_lhs->getOpPtrVector().push_back(
} else {
fe_lhs->getOpPtrVector().push_back(
}
}
}
auto set_hybridisation = [&](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
return -1;
};
op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(
boost::make_shared<double>(1.0), true, false));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact = [&](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
};
auto add_ops_contact_lhs = [&](auto &pip) {
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(new OpTreeSearch(
nullptr));
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
contact_sfd_map_range_ptr));
pip.push_back(
pip.push_back(
};
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
pip.push_back(op_loop_skeleton_side);
};
CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
CHKERR set_contact(fe_lhs->getOpPtrVector());
}
if (add_material) {
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
fe_rhs->getRuleHook = [](int, int, int) { return -1; };
fe_lhs->getRuleHook = [](int, int, int) { return -1; };
if (add_elastic) {
auto get_broken_op_side = [this](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
boost::shared_ptr<double> piola_scale_ptr(new double);
op_loop_domain_side->getOpPtrVector().push_back(
auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
piola_stress_mat_ptr));
pip.push_back(op_loop_domain_side);
return std::make_tuple(broken_data_ptr, piola_scale_ptr,
piola_stress_mat_ptr);
};
auto set_rhs = [&]() {
auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
get_broken_op_side(fe_rhs->getOpPtrVector());
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
fe_rhs->getOpPtrVector().push_back(
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
fe_rhs->getOpPtrVector().push_back(
hybrid_ptr));
auto get_normal_disp_bc_faces = [&]() {
auto faces =
return boost::make_shared<Range>(faces);
};
using BdyEleOp = BoundaryEle::UserDataOperator;
GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
get_normal_disp_bc_faces()));
};
auto set_lhs = [&]() {
auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
get_broken_op_side(fe_lhs->getOpPtrVector());
auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
fe_rhs->getOpPtrVector().push_back(
auto get_normal_disp_bc_faces = [&]() {
auto faces =
return boost::make_shared<Range>(faces);
};
using BdyEleOp = BoundaryEle::UserDataOperator;
GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
fe_lhs->getOpPtrVector().push_back(new OpC(
true, true, get_normal_disp_bc_faces()));
};
}
}
boost::shared_ptr<ContactTree> &fe_contact_tree
) {
auto get_body_range = [this](auto name, int dim, auto sev) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
MOFEM_LOG(
"EPSYNC", sev) <<
"Meshset: " << m_ptr->getMeshsetId() <<
" "
<< ents.size() << " entities";
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto calcs_side_traction = [&](auto &pip) {
using SideEleOp = EleOnSide::UserDataOperator;
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr(),
boost::make_shared<double>(1.0)));
pip.push_back(op_loop_domain_side);
};
auto add_contact_three = [&]() {
auto tree_moab_ptr = boost::make_shared<moab::Core>();
fe_contact_tree = boost::make_shared<ContactTree>(
get_body_range(
"CONTACT",
SPACE_DIM - 1, Sev::inform));
fe_contact_tree->getOpPtrVector().push_back(
contactDisp, contact_common_data_ptr->contactDispPtr()));
CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
fe_contact_tree->getOpPtrVector().push_back(
fe_contact_tree->getOpPtrVector().push_back(
new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
};
}
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_contact_bc = [&]() {
return op_loop_side;
};
}
boost::shared_ptr<FEMethod> null;
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
null);
null);
null);
null);
} else {
null);
null);
null);
null);
}
}
struct solve_elastic_setup {
bool set_ts_monitor) {
#ifdef ENABLE_PYTHON_BINDING
auto setup_sdf = [&]() {
boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
char sdf_file_name[255] = "sdf.py";
sdf_file_name, 255, PETSC_NULLPTR);
if (file_exists(sdf_file_name)) {
MOFEM_LOG(
"EP", Sev::inform) << sdf_file_name <<
" file found";
sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
MOFEM_LOG(
"EP", Sev::inform) <<
"SdfPython initialized";
} else {
MOFEM_LOG(
"EP", Sev::warning) << sdf_file_name <<
" file NOT found";
}
return sdf_python_ptr;
};
auto sdf_python_ptr = setup_sdf();
#endif
auto setup_ts_monitor = [&]() {
boost::shared_ptr<TsCtx>
ts_ctx;
if (set_ts_monitor) {
CHKERR TSMonitorSet(ts, TsMonitorSet,
ts_ctx.get(), PETSC_NULLPTR);
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"TS monitor setup";
return std::make_tuple(
ts_ctx);
};
auto setup_snes_monitor = [&]() {
SNES snes;
void *))MoFEMSNESMonitorEnergy,
(void *)(snes_ctx.get()), PETSC_NULLPTR);
MOFEM_LOG(
"EP", Sev::inform) <<
"SNES monitor setup";
};
auto setup_snes_conergence_test = [&]() {
auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
PetscReal snorm, PetscReal fnorm,
SNESConvergedReason *reason, void *cctx) {
CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
PETSC_NULLPTR);
CHKERR SNESGetSolutionUpdate(snes, &x_update);
CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
};
};
auto setup_section = [&]() {
PetscSection section_raw;
int num_fields;
CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
}
return section_raw;
};
auto set_vector_on_mesh = [&]() {
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG(
"EP", Sev::inform) <<
"Vector set on mesh";
};
auto setup_schur_block_solver = [&]() {
MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver";
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
if constexpr (A == AssemblyType::BLOCK_MAT) {
schur_ptr =
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Setting up Schur block solver done";
return schur_ptr;
};
#ifdef ENABLE_PYTHON_BINDING
return std::make_tuple(setup_sdf(), setup_ts_monitor(),
setup_snes_monitor(), setup_snes_conergence_test(),
setup_section(), set_vector_on_mesh(),
setup_schur_block_solver());
#else
return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
setup_snes_conergence_test(), setup_section(),
set_vector_on_mesh(), setup_schur_block_solver());
#endif
}
};
PetscBool debug_model = PETSC_FALSE;
PETSC_NULLPTR);
if (debug_model == PETSC_TRUE) {
auto post_proc = [&](TS ts, PetscReal
t, Vec u, Vec u_t, Vec u_tt, Vec
F,
void *ctx) {
SNES snes;
int it;
CHKERR SNESGetIterationNumber(snes, &it);
std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
std::string file_skel_name =
"snes_iteration_skel_" + std::to_string(it) + ".h5m";
auto get_material_force_tag = [&]() {
"can't get tag");
return tag;
};
{get_material_force_tag()});
};
ts_ctx_ptr->tsDebugHook = post_proc;
}
auto storage = solve_elastic_setup::setup(this, ts, x, true);
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
Vec xx;
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSSolve(ts, PETSC_NULLPTR);
TetPolynomialBase::switchCacheBaseOff<HDIV>(
SNES snes;
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
PETSC_NULLPTR);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
if (lin_solver_iterations > expected_lin_solver_iterations)
SETERRQ(
"Expected number of iterations is different than expected %d > %d",
lin_solver_iterations, expected_lin_solver_iterations);
}
PetscBool test_sslv116_flag = PETSC_FALSE;
&test_sslv116_flag, PETSC_NULLPTR);
if (test_sslv116_flag) {
double max_val = 0.0;
double min_val = 0.0;
auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
auto ent_type = ent_ptr->getEntType();
if (ent_type == MBVERTEX) {
max_val = std::max(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], max_val);
min_val = std::min(ent_ptr->getEntFieldData()[
SPACE_DIM - 1], min_val);
}
};
double global_max_val = 0.0;
double global_min_val = 0.0;
MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
double ref_max_val = 0.00767;
double ref_min_val = -0.00329;
if (std::abs(global_max_val - ref_max_val) > 1e-5) {
"Incorrect max value of the displacement field: %f != %f",
global_max_val, ref_max_val);
}
if (std::abs(global_min_val - ref_min_val) > 4e-5) {
"Incorrect min value of the displacement field: %f != %f",
global_min_val, ref_min_val);
}
}
}
auto storage = solve_elastic_setup::setup(this, ts, x, false);
double final_time = 1;
double delta_time = 0.1;
int max_it = 10;
PetscBool ts_h1_update = PETSC_FALSE;
PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
"none");
CHKERR PetscOptionsScalar(
"-dynamic_final_time",
"dynamic relaxation final time", "", final_time,
&final_time, PETSC_NULLPTR);
CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
"dynamic relaxation final time", "", delta_time,
&delta_time, PETSC_NULLPTR);
CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
max_it, &max_it, PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
PetscOptionsEnd();
auto setup_ts_monitor = [&]() {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
return monitor_ptr;
};
auto monitor_ptr = setup_ts_monitor();
TetPolynomialBase::switchCacheBaseOn<HDIV>(
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
if (ts_h1_update) {
}
monitor_ptr->ts_u = PETSC_NULLPTR;
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
if (!ts_h1_update) {
}
CHKERR TSSolve(ts, PETSC_NULLPTR);
if (!ts_h1_update) {
}
monitor_ptr->ts_u = PETSC_NULLPTR;
break;
}
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
auto set_block = [&](auto name, int dim) {
std::map<int, Range> map;
auto set_tag_impl = [&](auto name) {
std::regex((boost::format("%s(.*)") % name).str())
);
for (auto bc : bcs) {
true);
map[bc->getMeshsetId()] =
r;
}
};
return std::make_pair(name, map);
};
auto set_skin = [&](auto &&map) {
for (
auto &
m : map.second) {
}
return map;
};
auto set_tag = [&](auto &&map) {
auto name = map.first;
int def_val[] = {-1};
MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
"create tag");
for (
auto &
m : map.second) {
"clear tag");
}
};
set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
}
Vec f_residual, Vec var_vector,
std::vector<Tag> tags_to_transfer) {
if (rval == MB_SUCCESS) {
}
int def_val[] = {0};
"CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
int mark[] = {1};
tags_to_transfer.push_back(th);
auto get_tag = [&](auto name, auto dim) {
double def_val[] = {0., 0., 0.};
MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
"create tag");
return tag;
};
tags_to_transfer.push_back(get_tag("MaterialForce", 3));
}
tags_to_transfer.push_back(
t);
}
}
struct exclude_sdf {
exclude_sdf(
Range &&r) : map(
r) {}
bool operator()(
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (map.find(ent) != map.end()) {
return false;
}
return true;
}
private:
};
auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
auto post_proc_ptr =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
auto domain_ops = [&](auto &fe, int sense) {
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
auto piola_scale_ptr = boost::make_shared<double>(1.0);
fe.getOpPtrVector().push_back(
} else {
fe.getOpPtrVector().push_back(
}
if (var_vector) {
boost::make_shared<double>(1), vec));
boost::make_shared<double>(1), vec, MBMAXTYPE));
fe.getOpPtrVector().push_back(
} else {
fe.getOpPtrVector().push_back(
}
}
fe.getOpPtrVector().push_back(
if (auto op =
fe.getOpPtrVector().push_back(op);
}
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
struct OpSidePPMap :
public OpPPMap {
OpSidePPMap(moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
DataMapVec data_map_scalar, DataMapMat data_map_vec,
DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
int sense)
:
OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
data_map_vec, data_map_mat, data_symm_map_mat),
tagSense(sense) {}
if (tagSense != OpPPMap::getSkeletonSense())
CHKERR OpPPMap::doWork(side, type, data);
}
private:
int tagSense;
};
OpPPMap::DataMapMat vec_fields;
vec_fields[
"SpatialDisplacementL2"] =
dataAtPts->getSmallWL2AtPts();
vec_fields[
"SpatialDisplacementH1"] =
dataAtPts->getSmallWH1AtPts();
vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
vec_fields[
"ContactDisplacement"] =
dataAtPts->getContactL2AtPts();
vec_fields[
"AngularMomentum"] =
dataAtPts->getLeviKirchhoffAtPts();
vec_fields[
"X"] =
dataAtPts->getLargeXH1AtPts();
vec_fields[
"EiegnLogStreach"] =
dataAtPts->getEigenValsAtPts();
}
if (var_vector) {
vec_fields[
"VarOmega"] =
dataAtPts->getVarRotAxisPts();
vec_fields["VarSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], vec, MBTET));
}
if (f_residual) {
vec_fields["ResSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], vec, MBTET));
vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
rotAxis, vec_fields[
"ResOmega"], vec, MBTET));
}
OpPPMap::DataMapMat mat_fields;
mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
if (var_vector) {
mat_fields[
"VarPiolaStress"] =
dataAtPts->getVarPiolaPts();
}
if (f_residual) {
mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
boost::make_shared<double>(1), vec));
boost::make_shared<double>(1), vec, MBMAXTYPE));
}
case 0:
fe.getOpPtrVector().push_back(
break;
case 1:
fe.getOpPtrVector().push_back(
break;
default:
"Unsupported internal stress interpolation order %d",
}
}
OpPPMap::DataMapMat mat_fields_symm;
mat_fields_symm["LogSpatialStretch"] =
mat_fields_symm[
"SpatialStretch"] =
dataAtPts->getStretchTensorAtPts();
if (var_vector) {
mat_fields_symm["VarLogSpatialStretch"] =
}
if (f_residual) {
mat_fields_symm["ResLogSpatialStretch"] =
boost::make_shared<MatrixDouble>();
fe.getOpPtrVector().push_back(
MBTET));
}
}
fe.getOpPtrVector().push_back(
new OpSidePPMap(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
vec_fields,
mat_fields,
mat_fields_symm,
sense
)
);
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
};
post_proc_ptr->getOpPtrVector().push_back(
auto X_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
domain_ops(*(op_loop_side->getSideFEPtr()), sense);
post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
return post_proc_ptr;
};
auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
using SideEleOp = EleOnSide::UserDataOperator;
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr(),
boost::make_shared<double>(1.0)));
pip.push_back(op_loop_domain_side);
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
contactDisp, contact_common_data_ptr->contactDispPtr()));
pip.push_back(new OpTreeSearch(
&post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
if (f_residual) {
pip.push_back(op_this);
auto contact_residual = boost::make_shared<MatrixDouble>();
op_this->getOpPtrVector().push_back(
op_this->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"res_contact", contact_residual}},
{},
{}
)
);
}
};
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
post_proc_ptr->getOpPtrVector());
auto own_tets =
CommInterface::getPartEntities(mField.get_moab(), mField.get_comm_rank())
own_faces, moab::Interface::UNION);
auto get_post_negative = [&](auto &&ents) {
auto crack_faces_pos = ents;
auto crack_faces_neg = crack_faces_pos;
auto skin = get_skin(mField, own_tets);
auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
for (auto f : crack_on_proc_skin) {
tet = intersect(tet, own_tets);
int side_number, sense, offset;
CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
offset);
if (sense == 1) {
crack_faces_neg.erase(f);
} else {
crack_faces_pos.erase(f);
}
}
return std::make_pair(crack_faces_pos, crack_faces_neg);
};
auto get_crack_faces = [&](auto crack_faces) {
auto get_adj = [&](auto e, auto dim) {
CHKERR mField.get_moab().get_adjacencies(e, dim,
true, adj,
moab::Interface::UNION);
return adj;
};
auto tets = get_adj(crack_faces, 3);
auto faces = subtract(get_adj(tets, 2), crack_faces);
tets = subtract(tets, get_adj(faces, 3));
return subtract(crack_faces, get_adj(tets, 2));
};
auto [crack_faces_pos, crack_faces_neg] =
get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
auto only_crack_faces = [&](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
return false;
}
return true;
};
auto only_negative_crack_faces = [&](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
return false;
}
return true;
};
auto get_adj_front = [&]() {
auto skeleton_faces = *skeletonFaces;
CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2,
true, adj_front,
moab::Interface::UNION);
adj_front = intersect(adj_front, skeleton_faces);
adj_front = subtract(adj_front, *crackFaces);
adj_front = intersect(own_faces, adj_front);
return adj_front;
};
post_proc_ptr->setTagsToTransfer(tags_to_transfer);
post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
auto post_proc_begin =
post_proc_ptr->exeTestHook = only_crack_faces;
post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
post_proc_negative_sense_ptr, 0,
mField.get_comm_size());
constexpr bool debug =
false;
auto [adj_front_pos, adj_front_neg] =
auto only_front_faces_pos = [adj_front_pos](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (adj_front_pos.find(ent) == adj_front_pos.end()) {
return false;
}
return true;
};
auto only_front_faces_neg = [adj_front_neg](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (adj_front_neg.find(ent) == adj_front_neg.end()) {
return false;
}
return true;
};
post_proc_ptr->exeTestHook = only_front_faces_pos;
dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
post_proc_negative_sense_ptr, 0,
mField.get_comm_size());
}
CHKERR post_proc_end.writeFile(file.c_str());
}
Vec f_residual,
std::vector<Tag> tags_to_transfer) {
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_ptr =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
auto hybrid_disp = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
auto hybrid_res = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"hybrid_disp", hybrid_disp}},
{},
{}
)
);
if (f_residual) {
auto hybrid_res = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
hybridSpatialDisp, hybrid_res,
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"res_hybrid", hybrid_res}},
{},
{}
)
);
}
post_proc_ptr->setTagsToTransfer(tags_to_transfer);
auto post_proc_begin =
CHKERR post_proc_end.writeFile(file.c_str());
}
constexpr bool debug =
false;
auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
std::vector<Tag> tags;
tags.reserve(names.size());
auto create_and_clean = [&]() {
auto &tag = tags.back();
auto rval = moab.tag_get_handle(
n.first.c_str(), tag);
if (rval == MB_SUCCESS) {
moab.tag_delete(tag);
}
double def_val[] = {0., 0., 0.};
CHKERR moab.tag_get_handle(
n.first.c_str(),
n.second, MB_TYPE_DOUBLE,
tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
}
};
return tags;
};
enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
auto tags = get_tags_vec({{"MaterialForce", 3},
{"AreaGrowth", 3},
{"GriffithForce", 1},
{"FacePressure", 1}});
auto calculate_material_forces = [&]() {
auto get_face_material_force_fe = [&]() {
auto fe_ptr = boost::make_shared<FaceEle>(
mField);
fe_ptr->getRuleHook = [](int, int, int) { return -1; };
fe_ptr->setRuleHook =
fe_ptr->getOpPtrVector().push_back(
auto op_loop_domain_side =
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_domain_side->getOpPtrVector().push_back(
}
op_loop_domain_side->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
return fe_ptr;
};
auto integrate_face_material_force_fe = [&](auto &&face_energy_fe) {
auto face_exchange = CommInterface::createEntitiesPetscVector(
auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
int size;
CHKERR VecGetLocalSize(
v.second, &size);
int low, high;
CHKERR VecGetOwnershipRange(
v.second, &low, &high);
MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( "
<< low << " " << high << " ) ";
};
CHKERR print_loc_size(face_exchange,
"material face_exchange",
Sev::verbose);
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
"front_skin_faces_material_force_" +
}
};
CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
};
auto calculate_front_material_force = [&]() {
auto get_conn = [&](auto e) {
"get connectivity");
return conn;
};
auto get_conn_range = [&](auto e) {
"get connectivity");
return conn;
};
auto get_adj = [&](auto e, auto dim) {
"get adj");
return adj;
};
auto get_adj_range = [&](auto e, auto dim) {
moab::Interface::UNION),
"get adj");
return adj;
};
auto get_material_force = [&](
auto r,
auto th) {
MatrixDouble material_forces(
r.size(), 3,
false);
"get data");
return material_forces;
};
auto calculate_edge_direction = [&](auto e) {
int num_nodes;
"get connectivity");
std::array<double, 6> coords;
"get coords");
&coords[0], &coords[1], &coords[2]};
&coords[3], &coords[4], &coords[5]};
t_dir(
i) = t_p1(
i) - t_p0(
i);
return t_dir;
};
auto calculate_force_through_node = [&]() {
body_ents);
auto body_skin = get_skin(
mField, body_ents);
auto body_skin_conn = get_conn_range(body_skin);
for (
auto n : front_nodes) {
auto adj_tets = get_adj(
n, 3);
auto conn = get_conn_range(adj_tets);
adj_tets = get_adj_range(conn, 3);
}
auto skin_faces = get_skin(
mField, adj_tets);
auto material_forces = get_material_force(skin_faces, tags[0]);
all_skin_faces.merge(skin_faces);
}
auto calculate_node_material_force = [&]() {
auto t_face_T =
getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
for (auto face : skin_faces) {
t_face_force_tmp(
I) = t_face_T(
I);
++t_face_T;
auto face_tets = intersect(get_adj(face, 3), adj_tets);
if (face_tets.empty()) {
continue;
}
if (face_tets.size() != 1) {
"face_tets.size() != 1");
}
int side_number, sense, offset;
side_number, sense,
offset),
"moab side number");
t_face_force_tmp(
I) *= sense;
t_node_force(
I) += t_face_force_tmp(
I);
}
return t_node_force;
};
auto calculate_crack_area_growth_direction = [&](
auto n,
auto &t_node_force) {
auto boundary_node = intersect(
Range(
n,
n), body_skin_conn);
if (boundary_node.size()) {
auto faces = intersect(get_adj(
n, 2), body_skin);
t_project(
I) += t_normal_face(
I);
}
t_project.normalize();
}
auto adj_faces = intersect(get_adj(
n, 2), *
crackFaces);
if (adj_faces.empty()) {
auto adj_edges =
intersect(get_adj(
n, 1), unite(*
frontEdges, body_edges));
for (auto e : adj_edges) {
auto t_dir = calculate_edge_direction(e);
}
t_area_dir(
I) = -t_node_force(
I) / t_node_force.
l2();
return t_area_dir;
}
if (boundary_node.size()) {
t_Q(
I,
J) -= t_project(
I) * t_project(
J);
}
auto front_edges = get_adj(
n, 1);
for (
auto f : adj_faces) {
int num_nodes;
true);
std::array<double, 9> coords;
coords.data());
coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
auto n_it = std::find(conn, conn + num_nodes,
n);
auto n_index = std::distance(conn, n_it);
t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
t_d_normal(0, n_index, 2),
t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
t_d_normal(1, n_index, 2),
t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
t_d_normal(2, n_index, 2)};
t_projected_hessian(
I,
J) =
t_Q(
I, K) * (t_face_hessian(K, L) * t_Q(L,
J));
t_area_dir(K) +=
t_face_normal(
I) * t_projected_hessian(
I, K) / 2.;
}
return t_area_dir;
};
auto t_node_force = calculate_node_material_force();
&
n, 1, &t_node_force(0)),
"set data");
auto get_area_dir = [&]() {
auto t_area_dir =
calculate_crack_area_growth_direction(
n, t_node_force);
adj_edges.merge(intersect(get_adj_range(seed_n, 1),
seed_n = get_conn_range(adj_edges);
}
auto skin_adj_edges = get_skin(
mField, adj_edges);
seed_n = subtract(seed_n, skin_adj_edges);
for (auto sn : seed_n) {
auto t_area_dir_sn =
calculate_crack_area_growth_direction(sn, t_node_force);
t_area_dir(
I) += t_area_dir_sn(
I);
}
}
return t_area_dir;
};
auto t_area_dir = get_area_dir();
1, &t_area_dir(0)),
"set data");
auto griffith = -t_node_force(
I) * t_area_dir(
I) /
(t_area_dir(K) * t_area_dir(K));
"set data");
}
auto ave_node_force = [&](
auto th) {
auto conn = get_conn(e);
auto data = get_material_force(conn, th);
auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
++t_node;
}
t_edge(
I) /= conn.size();
calculate_edge_direction(e);
t_cross(K) =
}
};
auto ave_node_griffith_energy = [&](
auto th) {
tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
&e, 1, &t_edge_area_dir(0));
double griffith_energy = -t_edge_force(
I) * t_edge_area_dir(
I) /
(t_edge_area_dir(K) * t_edge_area_dir(K));
}
};
CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
};
CHKERR calculate_force_through_node();
auto adj_faces = get_adj(e, 2);
auto crack_face = intersect(get_adj(e, 2), *
crackFaces);
all_front_faces.merge(adj_faces);
}
&e, 1, &t_edge_force(0));
calculate_edge_direction(e);
for (
auto f : adj_faces) {
int side_number, sense, offset;
offset);
auto dot = -sense * t_cross(
I) * t_normal(
I);
tags[ExhangeTags::GRIFFITHFORCE], &
f, 1, &dot),
"set data");
}
}
int ts_step;
CHKERR TSGetStepNumber(ts, &ts_step);
"front_edges_material_force_" +
std::to_string(ts_step) + ".vtk",
"front_skin_faces_material_force_" +
std::to_string(ts_step) + ".vtk",
all_skin_faces);
"front_faces_material_force_" +
std::to_string(ts_step) + ".vtk",
all_front_faces);
}
}
auto edge_exchange = CommInterface::createEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
};
auto print_results = [&]() {
auto get_conn_range = [&](auto e) {
"get connectivity");
return conn;
};
auto get_tag_data = [&](auto &ents, auto tag, auto dim) {
std::vector<double> data(ents.size() * dim);
"get data");
return data;
};
auto at_nodes = [&]() {
auto material_force =
get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
auto griffith_force =
get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
std::vector<double> coords(conn.size() * 3);
"get coords");
if (conn.size())
MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Material force at nodes";
for (
size_t i = 0;
i < conn.size(); ++
i) {
<<
"Node " << conn[
i] <<
" coords " << coords[
i * 3 + 0] <<
" " << coords[
i * 3 + 1] <<
" " << coords[
i * 3 + 2] <<
" material force " << material_force[
i * 3 + 0] <<
" " << material_force[
i * 3 + 1] <<
" " << material_force[
i * 3 + 2] <<
" area growth " << area_growth[
i * 3 + 0] <<
" " << area_growth[
i * 3 + 1]
<<
" " << area_growth[
i * 3 + 2] <<
" griffith force " }
};
at_nodes();
}
};
CHKERR calculate_material_forces();
CHKERR calculate_front_material_force();
}
bool set_orientation) {
constexpr bool debug =
false;
constexpr auto sev = Sev::verbose;
auto body_skin = get_skin(
mField, body_ents);
moab::Interface::UNION);
Range boundary_skin_verts;
boundary_skin_verts, true);
Range geometry_edges_verts;
geometry_edges_verts, true);
true);
*
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
*
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
moab::Interface::UNION);
front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
auto get_tags_vec = [&](auto tag_name, int dim) {
std::vector<Tag> tags(1);
if (dim > 3)
auto create_and_clean = [&]() {
auto rval = moab.tag_get_handle(tag_name, tags[0]);
if (rval == MB_SUCCESS) {
moab.tag_delete(tags[0]);
}
double def_val[] = {0., 0., 0.};
CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
};
return tags;
};
auto get_adj_front = [&](bool subtract_crack) {
adj_front, moab::Interface::UNION);
if (subtract_crack)
return adj_front;
};
auto th_front_position = get_tags_vec("FrontPosition", 3);
auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
auto get_crack_adj_tets = [&](
auto r) {
Range crack_faces_conn_tets;
true, crack_faces_conn_tets,
moab::Interface::UNION);
return crack_faces_conn_tets;
};
auto get_layers_for_sides = [&](auto &side) {
std::vector<Range> layers;
auto get = [&]() {
auto get_adj = [&](
auto &
r,
int dim) {
moab::Interface::UNION);
return adj;
};
auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
true);
Range front_faces = get_adj(front_nodes, 2);
auto front_tets = get_tets(front_nodes);
auto front_side = intersect(side, front_tets);
layers.push_back(front_side);
for (;;) {
auto adj_faces = get_skin(
mField, layers.back());
adj_faces = intersect(adj_faces, front_faces);
auto adj_faces_tets = get_tets(adj_faces);
adj_faces_tets = intersect(adj_faces_tets, front_tets);
layers.push_back(unite(layers.back(), adj_faces_tets));
if (layers.back().size() == layers[layers.size() - 2].size()) {
break;
}
}
};
return layers;
};
auto layers_top = get_layers_for_sides(sides_pair.first);
auto layers_bottom = get_layers_for_sides(sides_pair.second);
#ifndef NDEBUG
"crack_tets_" +
sides_pair.second);
MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
for (auto &r : layers_top) {
MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
"layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
}
for (auto &r : layers_bottom) {
MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
"layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
}
}
#endif
auto get_cross = [&](
auto t_dir,
auto f) {
return t_cross;
};
auto get_sense = [&](
auto f,
auto e) {
int side, sense, offset;
"get sense");
return std::make_tuple(side, sense, offset);
};
auto calculate_edge_direction = [&](auto e, auto normalize = true) {
int num_nodes;
std::array<double, 6> coords;
&coords[0], &coords[1], &coords[2]};
&coords[3], &coords[4], &coords[5]};
t_dir(
i) = t_p1(
i) - t_p0(
i);
if (normalize)
return t_dir;
};
auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
auto front_faces,
auto &sides_pair,
auto th_position) {
th_face_energy);
th_material_force);
break;
default:
SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
"Unknown energy release selector");
};
auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
auto &edge_face_max_energy_map) {
auto body_skin = get_skin(
mField, body_ents);
for (auto e : front_edges) {
double griffith_force;
&griffith_force);
faces = subtract(intersect(faces, front_faces), body_skin);
std::vector<double> face_energy(faces.size());
face_energy.data());
auto max_energy_it =
std::max_element(face_energy.begin(), face_energy.end());
double max_energy =
max_energy_it != face_energy.end() ? *max_energy_it : 0;
edge_face_max_energy_map[e] =
std::make_tuple(faces[max_energy_it - face_energy.begin()],
griffith_force, static_cast<double>(0));
<< "Edge " << e << " griffith force " << griffith_force
<< " max face energy " << max_energy << " factor "
<< max_energy / griffith_force;
max_faces.insert(faces[max_energy_it - face_energy.begin()]);
}
#ifndef NDEBUG
"max_faces_" +
".vtk",
max_faces);
}
#endif
};
auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
auto up_down_face = [&](
auto &face_angle_map_up,
auto &face_angle_map_down
) {
for (
auto &
m : edge_face_max_energy_map) {
auto [max_face, energy, opt_angle] =
m.second;
faces = intersect(faces, front_faces);
false, adj_tets,
moab::Interface::UNION);
if (adj_tets.size()) {
false, adj_tets,
moab::Interface::UNION);
if (adj_tets.size()) {
adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
moab::Interface::UNION);
adj_tets_faces = intersect(adj_tets_faces, faces);
auto t_cross_max =
get_cross(calculate_edge_direction(e, true), max_face);
auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
t_cross_max(
i) *= sense_max;
for (
auto t : adj_tets) {
adj_tets_faces = intersect(adj_tets_faces, faces);
adj_tets_faces =
subtract(adj_tets_faces,
Range(max_face, max_face));
if (adj_tets_faces.size() == 1) {
auto t_cross = get_cross(calculate_edge_direction(e, true),
adj_tets_faces[0]);
auto [side, sense, offset] =
get_sense(adj_tets_faces[0], e);
double dot = t_cross(
i) * t_cross_max(
i);
auto angle = std::acos(dot);
double face_energy;
th_face_energy, adj_tets_faces, &face_energy);
auto [side_face, sense_face, offset_face] =
if (sense_face > 0) {
face_angle_map_up[e] = std::make_tuple(face_energy, angle,
adj_tets_faces[0]);
} else {
face_angle_map_down[e] = std::make_tuple(
face_energy, -angle, adj_tets_faces[0]);
}
}
}
}
}
}
};
auto calc_optimal_angle = [&](
auto &face_angle_map_up,
auto &face_angle_map_down
) {
for (
auto &
m : edge_face_max_energy_map) {
auto &[max_face, e0,
a0] =
m.second;
if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
face_angle_map_down.find(e) == face_angle_map_down.end()) {
} else {
th_material_force);
th_material_force, &e, 1, &t_material_force(0));
auto material_force_magnitude = t_material_force.
l2();
if (material_force_magnitude <
std::numeric_limits<double>::epsilon()) {
} else {
auto t_edge_dir = calculate_edge_direction(e, true);
auto t_cross_max = get_cross(t_edge_dir, max_face);
auto [side, sense, offset] = get_sense(max_face, e);
t_cross_max(sense) *= sense;
t_cross_max.normalize();
t_material_force(
J) * t_cross_max(K);
a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
<<
"Optimal angle " <<
a0 <<
" energy " << e0;
}
break;
}
default: {
SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
"Unknown energy release selector");
}
}
}
}
}
};
std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
face_angle_map_up;
std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
face_angle_map_down;
CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
#ifndef NDEBUG
auto th_angle = get_tags_vec("Angle", 1);
for (
auto &
m : face_angle_map_up) {
auto [e,
a, face] =
m.second;
up.insert(face);
}
for (
auto &
m : face_angle_map_down) {
auto [e,
a, face] =
m.second;
down.insert(face);
}
for (
auto &
m : edge_face_max_energy_map) {
auto [face, e, angle] =
m.second;
max_energy_faces.insert(face);
&angle);
}
max_energy_faces);
}
}
#endif
};
auto get_conn = [&](auto e) {
"get conn");
return conn;
};
auto get_adj = [&](auto e, auto dim) {
e, dim, false, adj, moab::Interface::UNION),
"get adj");
return adj;
};
auto get_coords = [&](
auto v) {
"get coords");
return t_coords;
};
auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
auto t_edge_dir = calculate_edge_direction(e, true);
auto [side, sense, offset] = get_sense(
f, e);
t_edge_dir.normalize();
t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
return std::make_tuple(t_normal, t_rotated_normal);
};
auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
auto &t_move, auto gamma) {
auto index = adj_vertex_tets_verts.index(
v);
if (index >= 0) {
for (auto ii : {0, 1, 2}) {
coords[3 * index + ii] += gamma * t_move(ii);
}
return true;
}
return false;
};
auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
auto &adj_vertex_tets, auto &coords) {
for (
auto t : adj_vertex_tets) {
int num_nodes;
std::array<double, 12> tet_coords;
for (
auto n = 0;
n != 4; ++
n) {
auto index = adj_vertex_tets_verts.index(conn[
n]);
if (index < 0) {
}
for (auto ii = 0; ii != 3; ++ii) {
tet_coords[3 *
n + ii] = coords[3 * index + ii];
}
}
double q = Tools::volumeLengthQuality(tet_coords.data());
if (!std::isnormal(q))
q = -2;
quality = std::min(quality, q);
};
return quality;
};
auto calculate_free_face_node_displacement =
[&](auto &edge_face_max_energy_map) {
auto get_vertex_edges = [&](auto vertex) {
auto impl = [&]() {
vertex_edges);
vertex_edges = subtract(vertex_edges, front_verts_edges);
if (boundary_skin_verts.size() &&
boundary_skin_verts.find(vertex[0]) !=
boundary_skin_verts.end()) {
vertex_edges = intersect(vertex_edges, body_skin_edges);
}
if (geometry_edges_verts.size() &&
geometry_edges_verts.find(vertex[0]) !=
geometry_edges_verts.end()) {
MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
vertex_edges = intersect(vertex_edges, geometry_edges);
}
if (crack_faces_verts.size() &&
crack_faces_verts.find(vertex[0]) !=
crack_faces_verts.end()) {
vertex_edges = intersect(vertex_edges, crack_faces_edges);
}
};
return vertex_edges;
};
using Bundle = std::vector<
>;
std::map<EntityHandle, Bundle> edge_bundle_map;
for (
auto &
m : edge_face_max_energy_map) {
auto &[max_face, energy, opt_angle] =
m.second;
auto [t_normal, t_rotated_normal] =
get_rotated_normal(edge, max_face, opt_angle);
auto front_vertex = get_conn(
Range(
m.first,
m.first));
auto adj_tets = get_adj(
Range(max_face, max_face), 3);
auto adj_tets_faces = get_adj(adj_tets, 2);
auto adj_front_faces = subtract(
intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
if (adj_front_faces.size() > 3)
"adj_front_faces.size()>3");
&t_material_force(0));
std::vector<double> griffith_energy(adj_front_faces.size());
th_face_energy, adj_front_faces, griffith_energy.data());
auto set_edge_bundle = [&](auto min_gamma) {
for (auto rotated_f : adj_front_faces) {
double rotated_face_energy =
griffith_energy[adj_front_faces.index(rotated_f)];
auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
front_vertex);
if (vertex.size() != 1) {
"Wrong number of vertex to move");
}
auto front_vertex_edges_vertex = get_conn(
intersect(get_adj(front_vertex, 1), crack_faces_edges));
vertex = subtract(
vertex, front_vertex_edges_vertex);
if (vertex.empty()) {
continue;
}
auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
auto whole_front =
subtract(body_skin_edges, crack_faces_edges));
auto front_edges =
subtract(get_adj(faces, 1), seen_front_edges);
if (front_edges.size() == 0) {
return 0;
}
auto front_connected_edges =
intersect(front_edges, whole_front);
if (front_connected_edges.size()) {
seen_front_edges.merge(front_connected_edges);
}
faces.merge(get_adj(front_edges, 2));
}
};
double rotated_face_cardinality = face_cardinality(
rotated_f,
seen_edges);
rotated_face_cardinality = std::max(rotated_face_cardinality,
1.);
auto t_vertex_coords = get_coords(vertex);
auto vertex_edges = get_vertex_edges(vertex);
for (auto e_used_to_move_detection : vertex_edges) {
auto edge_conn = get_conn(
Range(e_used_to_move_detection,
e_used_to_move_detection));
edge_conn = subtract(edge_conn, vertex);
t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
(t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
auto b =
(t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
std::numeric_limits<double>::epsilon();
if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
gamma > -0.1) {
t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
auto check_rotated_face_directoon = [&]() {
t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
auto dot =
(t_material_force(
i) / t_material_force.
l2()) *
return -dot > 0 ? true : false;
};
if (check_rotated_face_directoon()) {
<< "Crack edge " << edge << " moved face "
<< rotated_f
<< " edge: " << e_used_to_move_detection
<< " face direction/energy " << rotated_face_energy
<< " face cardinality " << rotated_face_cardinality
<< " gamma: " << gamma;
auto &bundle = edge_bundle_map[edge];
bundle.emplace_back(rotated_f, e_used_to_move_detection,
vertex[0], t_move, 1,
rotated_face_cardinality, gamma);
}
}
}
}
};
set_edge_bundle(std::numeric_limits<double>::epsilon());
if (edge_bundle_map[edge].empty()) {
set_edge_bundle(-1.);
}
}
return edge_bundle_map;
};
auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
sort_by_energy;
for (
auto &
m : edge_face_max_energy_map) {
auto &[max_face, energy, opt_angle] =
m.second;
sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
}
return sort_by_energy;
};
auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
"get tag");
auto get_face_pressure = [&](auto face) {
double pressure;
1, &pressure),
"get rag data");
return pressure;
};
<< "Number of edges to check " << sort_by_energy.size();
enum face_energy { POSITIVE, NEGATIVE };
constexpr bool skip_negative = true;
for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
++it) {
auto energy = it->first;
auto [max_edge, max_face, opt_angle] = it->second;
auto face_pressure = get_face_pressure(max_face);
if (skip_negative) {
if (fe == face_energy::POSITIVE) {
if (face_pressure < 0) {
<< "Skip negative face " << max_face << " with energy "
<< energy << " and pressure " << face_pressure;
continue;
}
}
}
<< "Check face " << max_face << " edge " << max_edge
<< " energy " << energy << " optimal angle " << opt_angle
<< " face pressure " << face_pressure;
auto jt = adj_edges_map.find(max_edge);
if (jt == adj_edges_map.end()) {
<< "Edge " << max_edge << " not found in adj_edges_map";
continue;
}
auto &bundle = jt->second;
auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
auto gamma) {
double max_quality = -2;
double max_quality_evaluated = -2;
double min_cardinality = std::numeric_limits<double>::max();
for (auto &b : bundle) {
auto &[face, move_edge, vertex, t_move, quality, cardinality,
edge_gamma] = b;
auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
std::vector<double> coords(3 * adj_vertex_tets_verts.size());
adj_vertex_tets_verts, coords.data()),
"get coords");
set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
quality = tets_quality(quality, adj_vertex_tets_verts,
adj_vertex_tets, coords);
auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
if (q < 0) {
return q;
} else {
return ((edge_gamma < 0) ? (q / 2) : q) / pow(
c, 2);
}
};
if (eval_quality(quality, cardinality, edge_gamma) >=
max_quality_evaluated) {
max_quality = quality;
min_cardinality = cardinality;
vertex_max = vertex;
face_max = face;
move_edge_max = move_edge;
t_move_last(
i) = t_move(
i);
max_quality_evaluated =
eval_quality(max_quality, min_cardinality, edge_gamma);
}
}
return std::make_tuple(vertex_max, face_max, t_move_last,
max_quality, min_cardinality);
};
auto find_max_in_bundle = [&](auto edge, auto &bundle) {
auto b_org_bundle = bundle;
auto r = find_max_in_bundle_impl(edge, bundle, 1.);
auto &[vertex_max, face_max, t_move_last, max_quality,
if (max_quality < 0) {
for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
bundle = b_org_bundle;
r = find_max_in_bundle_impl(edge, bundle, gamma);
auto &[vertex_max, face_max, t_move_last, max_quality,
<< "Back tracking: gamma " << gamma << " edge " << edge
<< " quality " << max_quality << " cardinality "
<< cardinality;
if (max_quality > 0.01) {
}
}
}
};
auto set_tag_to_vertex_and_face = [&](
auto &&
r,
auto &quality) {
auto &[
v,
f, t_move, q, cardinality] =
r;
if ((q > 0 && std::isnormal(q)) && energy > 0) {
<<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" " << max_edge << " move " << t_move << " energy " << energy
<< " quality " << q << " cardinality " << cardinality;
&t_move(0));
1, &energy);
}
quality = q;
};
double quality = -2;
CHKERR set_tag_to_vertex_and_face(
find_max_in_bundle(max_edge, bundle),
quality
);
if (quality > 0 && std::isnormal(quality) && energy > 0) {
<< "Crack face set with quality: " << quality;
}
}
if (!skip_negative)
break;
}
};
MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
edge_face_max_energy_map;
CHKERR find_maximal_face_energy(front_edges, front_faces,
edge_face_max_energy_map);
CHKERR calculate_face_orientation(edge_face_max_energy_map);
MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
calculate_free_face_node_displacement(edge_face_max_energy_map),
get_sort_by_energy(edge_face_max_energy_map)
);
};
CHKERR evaluate_face_energy_and_set_orientation(
*
frontEdges, get_adj_front(
false), sides_pair, th_front_position);
}
SCATTER_FORWARD);
SCATTER_FORWARD);
SCATTER_FORWARD);
auto get_max_moved_faces = [&]() {
auto adj_front = get_adj_front(false);
std::vector<double> face_energy(adj_front.size());
face_energy.data());
for (
int i = 0;
i != adj_front.size(); ++
i) {
if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
max_moved_faces.insert(adj_front[
i]);
}
}
return boost::make_shared<Range>(max_moved_faces);
};
#ifndef NDEBUG
"max_moved_faces_" +
}
#endif
}
auto rval =
std::vector<double> coords(3 * verts.size());
std::vector<double> pos(3 * verts.size());
for (
int i = 0;
i != 3 * verts.size(); ++
i) {
}
double zero[] = {0., 0., 0.};
}
#ifndef NDEBUG
constexpr bool debug =
false;
"set_coords_faces_" +
}
#endif
}
constexpr bool potential_crack_debug = false;
if constexpr (potential_crack_debug) {
true);
crack_front_verts);
true, crack_front_faces,
moab::Interface::UNION);
crack_front_faces = intersect(crack_front_faces, add_ents);
crack_front_faces);
}
auto get_crack_faces = [&]() {
} else {
}
};
auto get_extended_crack_faces = [&]() {
auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
ParallelComm *pcomm =
if (!pcomm->rank()) {
auto get_nodes = [&](auto &&e) {
"get connectivity");
return nodes;
};
auto get_adj = [&](auto &&e, auto dim,
auto t = moab::Interface::UNION) {
"get adj");
return adj;
};
body_ents);
auto body_skin = get_skin(
mField, body_ents);
auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
auto front_block_nodes = get_nodes(front_block_edges);
size_t s;
do {
s = crack_faces.size();
auto crack_face_nodes = get_nodes(crack_faces_org);
auto crack_faces_edges =
get_adj(crack_faces_org, 1, moab::Interface::UNION);
auto crack_skin = get_skin(
mField, crack_faces_org);
front_block_edges = subtract(front_block_edges, crack_skin);
auto crack_skin_nodes = get_nodes(crack_skin);
crack_skin_nodes.merge(front_block_nodes);
auto crack_skin_faces =
get_adj(crack_skin, 2, moab::Interface::UNION);
crack_skin_faces =
subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
crack_faces = crack_faces_org;
for (
auto f : crack_skin_faces) {
auto edges = intersect(
get_adj(
Range(
f,
f), 1, moab::Interface::UNION), crack_skin);
if (edges.size() == 2) {
edges.merge(
intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
body_skin_edges));
}
if (edges.size() == 2) {
auto edge_conn = get_nodes(
Range(edges));
auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
crack_faces_org);
if (faces.size() == 2) {
auto edge0_conn = get_nodes(
Range(edges[0], edges[0]));
auto edge1_conn = get_nodes(
Range(edges[1], edges[1]));
auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
crack_skin_nodes);
if (edges_conn.size() == 1) {
auto node_edges =
subtract(intersect(get_adj(edges_conn, 1,
moab::Interface::INTERSECT),
crack_faces_edges),
crack_skin);
if (node_edges.size()) {
auto get_t_dir = [&](auto e_conn) {
auto other_node = subtract(e_conn, edges_conn);
&t_dir(0));
return t_dir;
};
get_t_dir(edge0_conn)(
i) + get_t_dir(edge1_conn)(
i);
t_crack_surface_ave_dir(
i) = 0;
for (auto e : node_edges) {
auto e_conn = get_nodes(
Range(e, e));
auto t_dir = get_t_dir(e_conn);
t_crack_surface_ave_dir(
i) += t_dir(
i);
}
auto dot = t_ave_dir(
i) * t_crack_surface_ave_dir(
i);
if (dot < -std::numeric_limits<double>::epsilon()) {
}
} else {
}
}
}
} else if (edges.size() == 3) {
}
if (edges.size() == 1) {
edges.merge(
intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
geometry_edges));
edges.merge(
intersect(get_adj(
Range(
f,
f), 1, moab::Interface::UNION),
front_block_edges));
if (edges.size() == 2) {
continue;
}
}
}
crack_faces_org = crack_faces;
} while (s != crack_faces.size());
};
return crack_faces;
};
return get_faces_of_crack_front_verts(get_crack_faces());
};
get_extended_crack_faces());
}
auto reconstruct_crack_faces = [&](auto crack_faces) {
ParallelComm *pcomm =
auto impl = [&]() {
if (!pcomm->rank()) {
auto get_nodes = [&](auto &&e) {
"get connectivity");
return nodes;
};
auto get_adj = [&](auto &&e, auto dim,
auto t = moab::Interface::UNION) {
"get adj");
return adj;
};
auto get_test_on_crack_surface = [&]() {
auto crack_faces_nodes =
get_nodes(crack_faces);
auto crack_faces_tets =
get_adj(crack_faces_nodes, 3,
moab::Interface::UNION);
auto crack_faces_tets_nodes =
get_nodes(crack_faces_tets);
crack_faces_tets_nodes =
subtract(crack_faces_tets_nodes, crack_faces_nodes);
crack_faces_tets =
subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
moab::Interface::UNION));
new_crack_faces =
get_adj(crack_faces_tets, 2,
moab::Interface::UNION);
new_crack_faces.merge(crack_faces);
return std::make_tuple(new_crack_faces, crack_faces_tets);
};
auto carck_faces_test_edges = [&](auto faces, auto tets) {
auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
moab::Interface::UNION);
auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
adj_faces_edges.merge(geometry_edges);
adj_faces_edges.merge(front_block_edges);
auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
auto boundary_test_nodes = get_nodes(boundary_tets_edges);
auto boundary_test_nodes_edges =
get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
auto boundary_test_nodes_edges_nodes = subtract(
get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
boundary_tets_edges =
subtract(boundary_test_nodes_edges,
get_adj(boundary_test_nodes_edges_nodes, 1,
moab::Interface::UNION));
body_ents);
auto body_skin = get_skin(
mField, body_ents);
auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
body_skin_edges);
body_skin = intersect(body_skin, adj_tets_faces);
body_skin_edges = subtract(
body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
for (auto e : body_skin_edges) {
auto adj_tet = intersect(
get_adj(
Range(e, e), 3, moab::Interface::INTERSECT), tets);
if (adj_tet.size() == 1) {
boundary_tets_edges.insert(e);
}
}
return boundary_tets_edges;
};
auto p = get_test_on_crack_surface();
auto &[new_crack_faces, crack_faces_tets] = p;
crack_faces);
new_crack_faces);
crack_faces_tets);
}
auto boundary_tets_edges =
carck_faces_test_edges(new_crack_faces, crack_faces_tets);
boundary_tets_edges);
auto resolve_surface = [&](auto boundary_tets_edges,
auto crack_faces_tets) {
auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
auto crack_faces_tets_faces =
get_adj(crack_faces_tets, 2, moab::Interface::UNION);
int counter = 0;
int size = 0;
while (size != crack_faces_tets.size()) {
auto tets_faces =
get_adj(crack_faces_tets, 2, moab::Interface::UNION);
auto skin_tets = get_skin(
mField, crack_faces_tets);
auto skin_skin =
get_skin(
mField, subtract(crack_faces_tets_faces, tets_faces));
auto skin_skin_nodes = get_nodes(skin_skin);
size = crack_faces_tets.size();
<< "Crack faces tets size " << crack_faces_tets.size()
<< " crack faces size " << crack_faces_tets_faces.size();
auto skin_tets_nodes = subtract(
get_nodes(skin_tets),
boundary_tets_edges_nodes);
skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
for (
auto n : skin_tets_nodes) {
auto tets =
intersect(get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT),
crack_faces_tets);
if (tets.size() == 0) {
continue;
}
auto hole_detetction = [&]() {
auto adj_tets =
get_adj(
Range(
n,
n), 3, moab::Interface::INTERSECT);
adj_tets =
subtract(adj_tets,
crack_faces_tets);
if (adj_tets.size() == 0) {
return std::make_pair(
intersect(
get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
tets_faces),
tets);
}
std::vector<Range> tets_groups;
auto test_adj_tets = adj_tets;
while (test_adj_tets.size()) {
auto seed_size = 0;
Range seed =
Range(test_adj_tets[0], test_adj_tets[0]);
while (seed.size() != seed_size) {
auto adj_faces =
subtract(get_adj(seed, 2, moab::Interface::UNION),
tets_faces);
seed_size = seed.size();
seed.merge(
intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
test_adj_tets ));
}
tets_groups.push_back(seed);
test_adj_tets = subtract(test_adj_tets, seed);
}
if (tets_groups.size() == 1) {
return std::make_pair(
intersect(
get_adj(
Range(
n,
n), 2, moab::Interface::INTERSECT),
tets_faces),
tets);
}
for (auto &r : tets_groups) {
auto f = get_adj(r, 2, moab::Interface::UNION);
auto t = intersect(get_adj(
f, 3, moab::Interface::UNION),
crack_faces_tets);
if (
f.size() > faces_to_remove.size() ||
faces_to_remove.size() == 0) {
}
}
<< "Hole detection: faces to remove "
<< faces_to_remove.size() << " tets to remove "
<< tets_to_remove.size();
return std::make_pair(faces_to_remove, tets_to_remove);
};
if (tets.size() < tets_to_remove.size() ||
tets_to_remove.size() == 0) {
auto [h_faces_to_remove, h_tets_to_remove] =
hole_detetction();
faces_to_remove = h_faces_to_remove;
tets_to_remove = h_tets_to_remove;
}
all_removed_faces.merge(faces_to_remove);
all_removed_tets.merge(tets_to_remove);
}
crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
crack_faces_tets_faces =
subtract(crack_faces_tets_faces, faces_to_remove);
"removed_nodes_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
removed_nodes);
"faces_to_remove_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
faces_to_remove);
"tets_to_remove_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
tets_to_remove);
"crack_faces_tets_faces_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
crack_faces_tets_faces);
"crack_faces_tets_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
crack_faces_tets);
}
counter++;
}
auto cese_internal_faces = [&]() {
auto skin_tets = get_skin(
mField, crack_faces_tets);
auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
adj_faces =
subtract(adj_faces, skin_tets);
auto adj_tets = get_adj(adj_faces, 3,
moab::Interface::UNION);
crack_faces_tets =
subtract(crack_faces_tets,
adj_tets);
crack_faces_tets_faces =
subtract(crack_faces_tets_faces, adj_faces);
all_removed_faces.merge(adj_faces);
all_removed_tets.merge(adj_tets);
<< "Remove internal faces size " << adj_faces.size()
<< " tets size " << adj_tets.size();
};
auto case_only_one_free_edge = [&]() {
for (
auto t :
Range(crack_faces_tets)) {
auto adj_faces = get_adj(
moab::Interface::UNION);
auto crack_surface_edges =
get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
adj_faces),
1,
moab::Interface::UNION);
auto adj_edges =
subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
crack_surface_edges);
adj_edges = subtract(
adj_edges,
boundary_tets_edges);
if (adj_edges.size() == 1) {
crack_faces_tets =
subtract(crack_faces_tets,
auto faces_to_remove =
get_adj(adj_edges, 2, moab::Interface::UNION);
crack_faces_tets_faces =
subtract(crack_faces_tets_faces, faces_to_remove);
all_removed_faces.merge(faces_to_remove);
all_removed_tets.merge(
Range(
t,
t));
MOFEM_LOG(
"EPSELF", Sev::inform) <<
"Remove free one edges ";
}
}
crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
crack_faces_tets_faces =
subtract(crack_faces_tets_faces, all_removed_faces);
};
auto cese_flat_tet = [&](auto max_adj_edges) {
body_ents);
auto body_skin = get_skin(
mField, body_ents);
auto body_skin_edges =
get_adj(body_skin, 1, moab::Interface::UNION);
for (
auto t :
Range(crack_faces_tets)) {
auto adj_faces = get_adj(
moab::Interface::UNION);
auto crack_surface_edges =
get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
adj_faces),
1,
moab::Interface::UNION);
auto adj_edges =
subtract(get_adj(
Range(
t,
t), 1, moab::Interface::INTERSECT),
crack_surface_edges);
adj_edges = subtract(adj_edges, body_skin_edges);
auto tet_edges = get_adj(
Range(
t,
t), 1,
moab::Interface::UNION);
tet_edges = subtract(tet_edges, adj_edges);
for (auto e : tet_edges) {
constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
auto get_side = [&](auto e) {
int side, sense, offset;
"get side number failed");
return side;
};
auto get_side_ent = [&](auto side) {
"get side failed");
return side_edge;
};
adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
}
if (adj_edges.size() <= max_adj_edges) {
double dot = 1;
for (auto e : adj_edges) {
auto edge_adj_faces =
get_adj(
Range(e, e), 2, moab::Interface::UNION);
edge_adj_faces = intersect(edge_adj_faces, adj_faces);
if (edge_adj_faces.size() != 2) {
"Adj faces size is not 2 for edge " +
boost::lexical_cast<std::string>(e));
}
auto get_normal = [&](
auto f) {
"get tri normal failed");
return t_n;
};
auto t_n0 = get_normal(edge_adj_faces[0]);
auto t_n1 = get_normal(edge_adj_faces[1]);
auto get_sense = [&](
auto f) {
int side, sense, offset;
sense, offset),
"get side number failed");
return sense;
};
auto sense0 = get_sense(edge_adj_faces[0]);
auto sense1 = get_sense(edge_adj_faces[1]);
t_n0.normalize();
t_n1.normalize();
auto dot_e = (sense0 * sense1) * t_n0(
i) * t_n1(
i);
if (dot_e < dot || e == adj_edges[0]) {
dot = dot_e;
faces_to_remove = edge_adj_faces;
}
}
all_removed_faces.merge(faces_to_remove);
all_removed_tets.merge(
Range(
t,
t));
<< "Remove free edges on flat tet, with considered nb. of "
"edges "
<< adj_edges.size();
}
}
crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
crack_faces_tets_faces =
subtract(crack_faces_tets_faces, all_removed_faces);
};
"Case only one free edge failed");
for (auto max_adj_edges : {0, 1, 2, 3}) {
"Case only one free edge failed");
}
"Case internal faces failed");
"crack_faces_tets_faces_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
crack_faces_tets_faces);
"crack_faces_tets_" +
boost::lexical_cast<std::string>(counter) + ".vtk",
crack_faces_tets);
}
return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
all_removed_faces, all_removed_tets);
};
auto [resolved_faces, resolved_tets, all_removed_faces,
all_removed_tets] =
resolve_surface(boundary_tets_edges, crack_faces_tets);
resolved_faces.merge(subtract(crack_faces, all_removed_faces));
resolved_faces);
resolved_tets);
}
crack_faces = resolved_faces;
}
};
return crack_faces;
};
auto resolve_consisten_crack_extension = [&]() {
auto crack_meshset =
auto meshset = crack_meshset->getMeshset();
old_crack_faces);
auto extendeded_crack_faces = get_extended_crack_faces();
auto reconstructed_crack_faces =
subtract(reconstruct_crack_faces(extendeded_crack_faces),
<< "No new crack faces to add, skipping adding to meshset";
extendeded_crack_faces = subtract(
extendeded_crack_faces, subtract(*
crackFaces, old_crack_faces));
<< "Number crack faces size (extended) "
<< extendeded_crack_faces.size();
} else {
reconstructed_crack_faces);
<< "Number crack faces size (reconstructed) "
<< reconstructed_crack_faces.size();
}
}
crack_faces);
}
}
};
CHKERR resolve_consisten_crack_extension();
};
auto crack_faces =
verts = subtract(verts, conn);
std::vector<double> coords(3 * verts.size());
double def_coords[] = {0., 0., 0.};
"ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
}
};
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
return 2 * (p);
};
post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
post_proc_norm_fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
auto norms_vec =
CHKERR VecZeroEntries(norms_vec);
auto u_l2_ptr = boost::make_shared<MatrixDouble>();
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
u_h1_ptr));
auto piola_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
MBMAXTYPE));
post_proc_norm_fe->getOpPtrVector().push_back(
TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
*post_proc_norm_fe);
TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
CHKERR VecAssemblyBegin(norms_vec);
CHKERR VecAssemblyEnd(norms_vec);
const double *norms;
CHKERR VecGetArrayRead(norms_vec, &norms);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
<< "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
<< "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
}
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto disp_bc = bc.second->dispBcPtr) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
MOFEM_LOG(
"EP", Sev::noisy) <<
"Displacement BC: " << *disp_bc;
std::vector<double> block_attributes(6, 0.);
if (disp_bc->data.flag1 == 1) {
block_attributes[0] = disp_bc->data.value1;
block_attributes[3] = 1;
}
if (disp_bc->data.flag2 == 1) {
block_attributes[1] = disp_bc->data.value2;
block_attributes[4] = 1;
}
if (disp_bc->data.flag3 == 1) {
block_attributes[2] = disp_bc->data.value3;
block_attributes[5] = 1;
}
auto faces = bc.second->bcEnts.subset_by_dimension(2);
}
}
boost::make_shared<NormalDisplacementBcVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
auto block_name = "(.*)NORMAL_DISPLACEMENT(.*)";
std::regex reg_name(block_name);
if (std::regex_match(bc.first, reg_name)) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
block_name, bc.second->bcAttributes,
bc.second->bcEnts.subset_by_dimension(2));
}
}
boost::make_shared<AnalyticalDisplacementBcVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
auto block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
std::regex reg_name(block_name);
if (std::regex_match(bc.first, reg_name)) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
block_name, bc.second->bcAttributes,
bc.second->bcEnts.subset_by_dimension(2));
}
}
auto ts_displacement =
boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt");
<< "Add time scaling displacement BC: " << bc.blockName;
ts_displacement, "disp_history", ".txt", bc.blockName);
}
auto ts_normal_displacement =
boost::make_shared<DynamicRelaxationTimeScale>("normal_disp_history.txt");
<< "Add time scaling normal displacement BC: " << bc.blockName;
ts_normal_displacement, "normal_disp_history", ".txt",
bc.blockName);
}
}
false, false);
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto force_bc = bc.second->forceBcPtr) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
MOFEM_LOG(
"EP", Sev::noisy) <<
"Force BC: " << *force_bc;
std::vector<double> block_attributes(6, 0.);
block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
block_attributes[3] = 1;
block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
block_attributes[4] = 1;
block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
block_attributes[5] = 1;
auto faces = bc.second->bcEnts.subset_by_dimension(2);
faces);
}
}
for (auto bc : bc_mng->getBcMapByBlockName()) {
auto block_name = "(.*)PRESSURE(.*)";
std::regex reg_name(block_name);
if (std::regex_match(bc.first, reg_name)) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
block_name, bc.second->bcAttributes,
bc.second->bcEnts.subset_by_dimension(2));
}
}
boost::make_shared<AnalyticalTractionBcVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
auto block_name = "(.*)ANALYTICAL_TRACTION(.*)";
std::regex reg_name(block_name);
if (std::regex_match(bc.first, reg_name)) {
BcManager::extractStringFromBlockId(bc.first, "");
<<
"Field name: " <<
field_name <<
" Block name: " << block_name;
block_name, bc.second->bcAttributes,
bc.second->bcEnts.subset_by_dimension(2));
}
}
auto ts_traction =
boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
ts_traction, "traction_history", ".txt", bc.blockName);
}
auto ts_pressure =
boost::make_shared<DynamicRelaxationTimeScale>("pressure_history.txt");
ts_pressure, "pressure_history", ".txt", bc.blockName);
}
}
auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
const std::string block_name,
const int nb_attributes) {
std::regex((boost::format("%s(.*)") % block_name).str()))
) {
std::vector<double> block_attributes;
CHKERR it->getAttributes(block_attributes);
if (block_attributes.size() < nb_attributes) {
"In block %s expected %d attributes, but given %ld",
it->getName().c_str(), nb_attributes, block_attributes.size());
}
auto get_block_ents = [&]() {
true);
return ents;
};
auto Ents = get_block_ents();
ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
get_block_ents());
}
};
auto ts_pre_stretch =
boost::make_shared<DynamicRelaxationTimeScale>("externalstrain_history.txt");
<< "Add time scaling external strain: " << ext_strain_block.blockName;
ts_pre_stretch, "externalstrain_history", ".txt", ext_strain_block.blockName);
}
}
auto print_loc_size = [
this](
auto v,
auto str,
auto sev) {
int size;
CHKERR VecGetLocalSize(
v.second, &size);
int low, high;
CHKERR VecGetOwnershipRange(
v.second, &low, &high);
MOFEM_LOG(
"EPSYNC", sev) << str <<
" local size " << size <<
" ( " << low
<< " " << high << " ) ";
};
}
if (!area_ptr) {
}
int success;
*area_ptr = 0;
MOFEM_LOG(
"EP", Sev::inform) <<
"Calculate crack area";
for (
auto f : crack_faces) {
}
success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
} else {
success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0,
mField.
get_comm());
}
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_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_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.
#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
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
auto id_from_handle(const EntityHandle h)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
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.
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr double t
plate stiffness
constexpr auto field_name
#define QUAD_2D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_3D_TABLE[]
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
FTensor::Index< 'm', 3 > m
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static enum StretchSelector stretchSelector
boost::shared_ptr< Range > frontAdjEdges
MoFEMErrorCode createCrackSurfaceMeshset()
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
const std::string skeletonElement
static double inv_f_linear(const double v)
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
boost::shared_ptr< Range > contactFaces
static double dd_f_log_e_quadratic(const double v)
static double dynamicTime
static double dd_f_log(const double v)
BitRefLevel bitAdjEnt
bit ref level for parent
MoFEM::Interface & mField
const std::string spatialL2Disp
static double inv_d_f_log(const double v)
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
static PetscBool l2UserBaseScale
SmartPetscObj< DM > dM
Coupled problem all fields.
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
static int internalStressInterpOrder
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
const std::string materialH1Positions
MoFEMErrorCode setBlockTagsOnSkin()
std::vector< Tag > listTagsToTransfer
list of tags to transfer to postprocessor
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
static PetscBool crackingOn
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
Remove all, but entities where kinematic constrains are applied.
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates, SmartPetscObj< Vec > ver_vec, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe)
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
static double griffithEnergy
Griffith energy.
MoFEMErrorCode calculateCrackArea(boost::shared_ptr< double > area_ptr)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
const std::string elementVolumeName
static double dd_f_log_e(const double v)
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
CommInterface::EntitiesPetscVector vertexExchange
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
boost::shared_ptr< Range > maxMovedFaces
static PetscBool dynamicRelaxation
const std::string spatialH1Disp
MoFEMErrorCode solveElastic(TS ts, Vec x)
static double d_f_log(const double v)
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
static double crackingStartTime
MoFEMErrorCode getOptions()
const std::string piolaStress
MoFEMErrorCode setElasticElementToTs(DM dm)
int contactRefinementLevels
MoFEMErrorCode gettingNorms()
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)
const std::string bubbleField
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x)
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation)
static double inv_f_log(const double v)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
static PetscBool noStretch
MoFEMErrorCode setNewFrontCoordinates()
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
static double exponentBase
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
static double dd_f_linear(const double v)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
static boost::function< double(const double)> d_f
boost::shared_ptr< Range > skeletonFaces
boost::shared_ptr< PhysicalEquations > physicalEquations
static boost::function< double(const double)> f
const std::string rotAxis
static double inv_d_f_linear(const double v)
BitRefLevel bitAdjParentMask
bit ref level for parent parent
static double inv_dd_f_log(const double v)
const std::string contactDisp
static std::string internalStressTagName
static enum SymmetrySelector symmetrySelector
CommInterface::EntitiesPetscVector edgeExchange
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
static boost::function< double(const double)> inv_d_f
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
static boost::function< double(const double)> inv_f
const std::string skinElement
static PetscBool internalStressVoigt
static double inv_dd_f_linear(const double v)
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
static PetscBool setSingularity
static double d_f_log_e(const double v)
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
static double f_log_e_quadratic(const double v)
static int nbJIntegralLevels
MoFEMErrorCode addCrackSurfaces(const bool debug=false)
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
MoFEMErrorCode getSpatialDispBc()
BitRefLevel bitAdjParent
bit ref level for parent
MoFEMErrorCode postProcessResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, Vec var_vec=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
static double d_f_log_e_quadratic(const double v)
CommInterface::EntitiesPetscVector volumeExchange
MoFEMErrorCode saveOrgCoords()
const std::string naturalBcElement
static double f_log_e(const double v)
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
static double d_f_linear(const double v)
const std::string hybridSpatialDisp
SmartPetscObj< Vec > solTSStep
static double f_log(const double v)
CommInterface::EntitiesPetscVector faceExchange
SmartPetscObj< DM > dmElastic
Elastic problem.
static boost::function< double(const double)> inv_dd_f
boost::shared_ptr< Range > frontEdges
static boost::function< double(const double)> dd_f
const std::string stretchTensor
BitRefLevel bitAdjEntMask
bit ref level for parent parent
static double f_linear(const double v)
const std::string contactElement
static boost::function< double(const double)> inv_dd_f
static double f_log_e(const double v)
static double inv_dd_f_log_e(const double v)
static double inv_d_f_log_e(const double v)
static PetscBool setSingularity
static double d_f_log_e(const double v)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
static double exponentBase
static boost::function< double(const double)> f
static double dd_f_log_e(const double v)
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
static boost::function< double(const double)> inv_d_f
EshelbianCore(MoFEM::Interface &m_field)
static double inv_f_log_e(const double v)
static boost::function< double(const double)> inv_f
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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
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
Set integration rule to boundary elements.
int operator()(int, int, int) const
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.
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.
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
Provide data structure for (tensor) field approximation.
structure to get information form mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
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.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Apply rotation boundary condition.
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)
Set integration rule to volume elements.
int operator()(int, int, int) const
BoundaryEle::UserDataOperator BdyEleOp
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
