Eshelbian plasticity implementation.
Eshelbian plasticity implementationmofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp
#define SINGULARITY
#ifdef INCLUDE_MBCOUPLER
#include <mbcoupler/Coupler.hpp>
#endif
#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
};
struct FaceUserDataOperatorStabAssembly
};
}
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_nodes,
boost::shared_ptr<Range> front_edges,
boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi = nullptr)
boost::shared_ptr<Range> front_nodes,
boost::shared_ptr<Range> front_edges,
FunRule fun_rule,
boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi = nullptr)
int order_col, int order_data) {
constexpr bool debug =
false;
constexpr int numNodes = 4;
constexpr int numEdges = 6;
constexpr int refinementLevels = 6;
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 = [&]() {
"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);
}
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++) {
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(
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
->updateMeshsetByEntitiesChildren(meshset,
meshset, MBEDGE, true);
}
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
tets);
}
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());
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);
TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
};
}
}
}
}
private:
using ForcesAndSourcesCore::dataOnElement;
private:
using ForcesAndSourcesCore::ForcesAndSourcesCore;
};
boost::shared_ptr<CGGUserPolynomialBase::CachePhi>
cachePhi;
static inline std::map<long int, MatrixDouble>
mapRefCoords;
};
struct SetIntegrationAtFrontFace {
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 = 6;
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);
}
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++) {
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(
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
->updateMeshsetByEntitiesChildren(meshset,
meshset, MBEDGE, true);
}
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
tris);
}
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());
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;
}
if (evalRhs)
if (evalLhs)
}
}
const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
const char *list_symm[] = {"symm", "not_symm"};
const char *list_release[] = {"griffith_force", "griffith_skeleton"};
const char *list_stretches[] = {"linear", "log", "log_quadratic"};
const char *list_solvers[] = {"time_solver", "dynamic_relaxation",
"cohesive"};
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,
PETSC_NULLPTR);
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 PetscOptionsEList(
"-solver_type",
"solver type",
"", list_solvers,
&choice_solver, PETSC_NULLPTR);
CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-cohesive_interface_on",
"cohesive interface ON",
"",
CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
PETSC_NULLPTR);
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",
PETSC_NULLPTR);
"-nb_J_integral_contours", "Number of J integral contours", "",
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",
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;
};
<< "-dynamic_relaxation option is deprecated, use -solver_type "
"dynamic_relaxation instead.";
}
switch (choice_solver) {
break;
break;
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";
<< "Solver type: -solver_type " << list_solvers[choice_solver];
? "yes"
: "no";
? "yes"
: "no";
? "yes"
: "no";
<< "Energy release variant: -energy_release_variant "
<< "Number of J integral contours: -nb_J_integral_contours "
<< "Cohesive interface on: -cohesive_interface_on "
? "yes"
: "no";
#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;
moab.get_connectivity(front_edges, front_crack_nodes, true),
"get_connectivity failed");
false, crack_front_edges,
moab::Interface::UNION),
"get_adjacencies failed");
Range crack_front_edges_nodes;
crack_front_edges_nodes, true),
"get_connectivity failed");
crack_front_edges_nodes =
subtract(crack_front_edges_nodes, front_crack_nodes);
Range crack_front_edges_with_both_nodes_not_at_front;
moab.get_adjacencies(crack_front_edges_nodes, 1, false,
crack_front_edges_with_both_nodes_not_at_front,
moab::Interface::UNION),
"get_adjacencies failed");
crack_front_edges_with_both_nodes_not_at_front = intersect(
crack_front_edges, crack_front_edges_with_both_nodes_not_at_front);
}
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);
<<
"Number of crack faces: " <<
crackFaces->size();
<<
"Number of front edges: " <<
frontEdges->size();
#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;
[&](boost::shared_ptr<FieldEntity> field_entity_ptr) ->
MoFEMErrorCode {
auto nb_dofs = field_entity_ptr->getEntFieldData().size();
if (nb_dofs == 0) {
}
#ifndef NDEBUG
if (field_entity_ptr->getNbOfCoeffs() != 3)
"Expected 3 coefficients per edge");
if (nb_dofs % 3 != 0)
"Expected multiple of 3 coefficients per edge");
#endif
auto get_conn = [&]() {
int num_nodes;
CHKERR moab.get_connectivity(field_entity_ptr->getEnt(), conn,
num_nodes, false);
return std::make_pair(conn, num_nodes);
};
auto get_dir = [&](auto &&conn_p) {
auto [conn, num_nodes] = conn_p;
double coords[6];
CHKERR moab.get_coords(conn, num_nodes, coords);
coords[4] - coords[1],
coords[5] - coords[2]};
return t_edge_dir;
};
auto get_singularity_dof = [&](auto &&conn_p, auto &&t_edge_dir) {
auto [conn, num_nodes] = conn_p;
if (front_vertices.find(conn[0]) != front_vertices.end()) {
t_singularity_dof(
i) = t_edge_dir(
i) * (-
eps);
} else if (front_vertices.find(conn[1]) != front_vertices.end()) {
t_singularity_dof(
i) = t_edge_dir(
i) *
eps;
}
return t_singularity_dof;
};
auto t_singularity_dof =
get_singularity_dof(get_conn(), get_dir(get_conn()));
auto field_data = field_entity_ptr->getEntFieldData();
&field_data[0], &field_data[1], &field_data[2]};
t_dof(
i) = t_singularity_dof(
i);
++t_dof;
for (
auto n = 1;
n < field_data.size() / 3; ++
n) {
++t_dof;
}
};
&front_adj_edges);
};
}
#ifdef INCLUDE_MBCOUPLER
char mesh_source_file[255] = "source.h5m";
char iterp_tag_name[255] = "INTERNAL_STRESS";
int interp_order = 1;
PetscBool hybrid_interp = PETSC_TRUE;
PetscBool project_internal_stress = PETSC_FALSE;
double toler = 5.e-10;
PetscOptionsBegin(PETSC_COMM_WORLD, "mesh_transfer_", "mesh data transfer",
"none");
CHKERR PetscOptionsString(
"-source_file",
"source mesh file name",
"",
"source.h5m", mesh_source_file, 255,
&project_internal_stress);
CHKERR PetscOptionsString(
"-interp_tag",
"interpolation tag name",
"",
"INTERNAL_STRESS", iterp_tag_name, 255,
PETSC_NULLPTR);
CHKERR PetscOptionsInt(
"-interp_order",
"interpolation order",
"", 0,
&interp_order, PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-hybrid_interp",
"use hybrid interpolation",
"",
hybrid_interp, &hybrid_interp, PETSC_NULLPTR);
PetscOptionsEnd();
if (!project_internal_stress) {
<< "No internal stress source mesh specified. Skipping projection of "
"internal stress";
}
<< "Projecting internal stress from source mesh: " << mesh_source_file;
auto rval_check_tag = moab.tag_get_handle(iterp_tag_name, old_interp_tag);
if (rval_check_tag == MB_SUCCESS) {
<< "Deleting existing tag on target mesh: " << iterp_tag_name;
CHKERR moab.tag_delete(old_interp_tag);
}
int world_rank = -1, world_size = -1;
MPI_Comm_rank(PETSC_COMM_WORLD, &world_rank);
MPI_Comm_size(PETSC_COMM_WORLD, &world_size);
Range original_meshset_ents;
CHKERR moab.get_entities_by_handle(0, original_meshset_ents);
MPI_Comm comm_coupler;
if (world_rank == 0) {
MPI_Comm_split(PETSC_COMM_WORLD, 0, 0, &comm_coupler);
} else {
MPI_Comm_split(PETSC_COMM_WORLD, MPI_UNDEFINED, world_rank, &comm_coupler);
}
ParallelComm *pcomm0 = nullptr;
int pcomm0_id = -1;
if (world_rank == 0) {
pcomm0 = new ParallelComm(&moab, comm_coupler, &pcomm0_id);
}
Coupler::Method method;
switch (interp_order) {
case 0:
method = Coupler::CONSTANT;
break;
case 1:
method = Coupler::LINEAR_FE;
break;
default:
"Unsupported interpolation order");
}
int nprocs, rank;
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
CHKERR moab.create_meshset(MESHSET_SET, target_root);
<< "Creating target mesh from existing meshset";
Range target_meshset_ents;
CHKERR moab.get_entities_by_handle(0, target_meshset_ents);
CHKERR moab.add_entities(target_root, target_meshset_ents);
int tag_length;
DataType dtype;
TagType storage;
Range targ_verts, targ_elems;
if (world_rank == 0) {
CHKERR moab.create_meshset(MESHSET_SET, source_root);
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Loading source mesh on rank 0";
auto rval_source_mesh = moab.load_file(mesh_source_file, &source_root, "");
if (rval_source_mesh != MB_SUCCESS) {
<< "Error loading source mesh file: " << mesh_source_file;
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Source mesh loaded.";
CHKERR moab.get_entities_by_dimension(source_root, 3, src_elems);
CHKERR pcomm0->create_part(part_set);
CHKERR moab.add_entities(part_set, src_elems);
CHKERR pcomm0->get_part_entities(src_elems_part, 3);
CHKERR moab.tag_get_handle(iterp_tag_name, interp_tag);
int interp_tag_len;
CHKERR moab.tag_get_length(interp_tag, interp_tag_len);
if (interp_tag_len != 1 && interp_tag_len != 3 && interp_tag_len != 9) {
"Unsupported interpolation tag length: %d", interp_tag_len);
}
tag_length = interp_tag_len;
CHKERR moab.tag_get_data_type(interp_tag, dtype);
CHKERR moab.tag_get_type(interp_tag, storage);
Coupler mbc(&moab, pcomm0, src_elems_part, 0, true);
std::vector<double> vpos;
int num_pts = 0;
CHKERR moab.get_entities_by_dimension(target_root, 3, targ_elems);
if (interp_order == 0) {
targ_verts = targ_elems;
} else {
CHKERR moab.get_adjacencies(targ_elems, 0,
false, targ_verts,
moab::Interface::UNION);
}
CHKERR pcomm0->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
targ_verts = subtract(targ_verts, tmp_verts);
num_pts = (
int)targ_verts.size();
vpos.resize(3 * targ_verts.size());
CHKERR moab.get_coords(targ_verts, &vpos[0]);
boost::shared_ptr<TupleList> tl_ptr;
tl_ptr = boost::make_shared<TupleList>();
CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
false);
auto find_missing_points = [&](
Range &targ_verts,
int &num_pts,
std::vector<double> &vpos,
int missing_pts_num = 0;
auto vit = targ_verts.begin();
for (; vit != targ_verts.end();
i++) {
if (tl_ptr->vi_rd[3 *
i + 1] == -1) {
missing_verts.insert(*vit);
vit = targ_verts.erase(vit);
missing_pts_num++;
} else {
vit++;
}
}
int missing_pts_num_global = 0;
if (missing_pts_num_global) {
<< missing_pts_num_global
<< " points in target mesh were not located in source mesh. ";
}
if (missing_pts_num) {
num_pts = (
int)targ_verts.size();
vpos.resize(3 * targ_verts.size());
CHKERR moab.get_coords(targ_verts, &vpos[0]);
tl_ptr->reset();
CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
false);
}
};
CHKERR find_missing_points(targ_verts, num_pts, vpos, missing_verts);
std::vector<double> source_data(interp_tag_len * src_elems.size(), 0.0);
std::vector<double> target_data(interp_tag_len * num_pts, 0.0);
CHKERR moab.tag_get_data(interp_tag, src_elems, &source_data[0]);
Tag scalar_tag, adj_count_tag;
double def_scl = 0;
string scalar_tag_name = string(iterp_tag_name) + "_COMP";
CHKERR moab.tag_get_handle(scalar_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
scalar_tag, MB_TAG_CREAT | MB_TAG_DENSE,
&def_scl);
string adj_count_tag_name = "ADJ_COUNT";
double def_adj = 0;
CHKERR moab.tag_get_handle(adj_count_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
adj_count_tag, MB_TAG_CREAT | MB_TAG_DENSE,
&def_adj);
auto create_scalar_tags = [&](
const Range &src_elems,
const std::vector<double> &source_data,
int itag) {
std::vector<double> source_data_scalar(src_elems.size());
for (int ielem = 0; ielem < src_elems.size(); ielem++) {
source_data_scalar[ielem] = source_data[itag + ielem * interp_tag_len];
}
CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
if (interp_order == 1) {
CHKERR moab.get_connectivity(src_elems, src_verts,
true);
CHKERR moab.tag_clear_data(scalar_tag, src_verts, &def_scl);
CHKERR moab.tag_clear_data(adj_count_tag, src_verts, &def_adj);
for (auto &tet : src_elems) {
double tet_data = 0;
CHKERR moab.tag_get_data(scalar_tag, &tet, 1, &tet_data);
CHKERR moab.get_connectivity(&tet, 1, adj_verts,
true);
std::vector<double> adj_vert_data(adj_verts.size(), 0.0);
std::vector<double> adj_vert_count(adj_verts.size(), 0.0);
CHKERR moab.tag_get_data(scalar_tag, adj_verts, &adj_vert_data[0]);
CHKERR moab.tag_get_data(adj_count_tag, adj_verts,
&adj_vert_count[0]);
for (int ivert = 0; ivert < adj_verts.size(); ivert++) {
adj_vert_data[ivert] += tet_data;
adj_vert_count[ivert] += 1;
}
CHKERR moab.tag_set_data(scalar_tag, adj_verts, &adj_vert_data[0]);
CHKERR moab.tag_set_data(adj_count_tag, adj_verts,
&adj_vert_count[0]);
}
std::vector<Tag> tags = {scalar_tag, adj_count_tag};
pcomm0->reduce_tags(tags, tags, MPI_SUM, src_verts);
std::vector<double> src_vert_data(src_verts.size(), 0.0);
std::vector<double> src_vert_adj_count(src_verts.size(), 0.0);
CHKERR moab.tag_get_data(scalar_tag, src_verts, &src_vert_data[0]);
CHKERR moab.tag_get_data(adj_count_tag, src_verts,
&src_vert_adj_count[0]);
for (int ivert = 0; ivert < src_verts.size(); ivert++) {
src_vert_data[ivert] /= src_vert_adj_count[ivert];
}
CHKERR moab.tag_set_data(scalar_tag, src_verts, &src_vert_data[0]);
}
};
for (int itag = 0; itag < interp_tag_len; itag++) {
CHKERR create_scalar_tags(src_elems, source_data, itag);
std::vector<double> target_data_scalar(num_pts, 0.0);
CHKERR mbc.interpolate(method, scalar_tag_name, &target_data_scalar[0],
tl_ptr.get());
for (int ielem = 0; ielem < num_pts; ielem++) {
target_data[itag + ielem * interp_tag_len] = target_data_scalar[ielem];
}
}
CHKERR moab.tag_set_data(interp_tag, targ_verts, &target_data[0]);
if (missing_verts.size() && (interp_order == 1) && hybrid_interp) {
MOFEM_LOG(
"WORLD", Sev::warning) <<
"Using hybrid interpolation for "
"missing points in the target mesh.";
CHKERR moab.get_adjacencies(missing_verts, 3,
false, missing_adj_elems,
moab::Interface::UNION);
int num_adj_elems = (
int)missing_adj_elems.size();
std::vector<double> vpos_adj_elems;
vpos_adj_elems.resize(3 * missing_adj_elems.size());
CHKERR moab.get_coords(missing_adj_elems, &vpos_adj_elems[0]);
tl_ptr->reset();
CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
tl_ptr.get(), false);
CHKERR find_missing_points(missing_adj_elems, num_adj_elems,
vpos_adj_elems, missing_tets);
if (missing_tets.size()) {
<< missing_tets.size()
<< "points in target mesh were not located in source mesh. ";
}
std::vector<double> target_data_adj_elems(interp_tag_len * num_adj_elems,
0.0);
for (int itag = 0; itag < interp_tag_len; itag++) {
CHKERR create_scalar_tags(src_elems, source_data, itag);
std::vector<double> target_data_adj_elems_scalar(num_adj_elems, 0.0);
CHKERR mbc.interpolate(method, scalar_tag_name,
&target_data_adj_elems_scalar[0], tl_ptr.get());
for (int ielem = 0; ielem < num_adj_elems; ielem++) {
target_data_adj_elems[itag + ielem * interp_tag_len] =
target_data_adj_elems_scalar[ielem];
}
}
CHKERR moab.tag_set_data(interp_tag, missing_adj_elems,
&target_data_adj_elems[0]);
for (auto &vert : missing_verts) {
CHKERR moab.get_adjacencies(&vert, 1, 3,
false, adj_elems,
moab::Interface::UNION);
std::vector<double> adj_elems_data(adj_elems.size() * interp_tag_len,
0.0);
CHKERR moab.tag_get_data(interp_tag, adj_elems, &adj_elems_data[0]);
std::vector<double> vert_data(interp_tag_len, 0.0);
for (int itag = 0; itag < interp_tag_len; itag++) {
for (
int i = 0;
i < adj_elems.size();
i++) {
vert_data[itag] += adj_elems_data[
i * interp_tag_len + itag];
}
vert_data[itag] /= adj_elems.size();
}
CHKERR moab.tag_set_data(interp_tag, &vert, 1, &vert_data[0]);
}
}
CHKERR moab.tag_delete(scalar_tag);
CHKERR moab.tag_delete(adj_count_tag);
CHKERR moab.get_entities_by_handle(source_root, src_mesh_ents);
CHKERR moab.delete_entities(&source_root, 1);
CHKERR moab.delete_entities(src_mesh_ents);
CHKERR moab.delete_entities(&part_set, 1);
}
MPI_Bcast(&tag_length, 1, MPI_INT, 0, PETSC_COMM_WORLD);
MPI_Bcast(&dtype, 1, MPI_INT, 0, PETSC_COMM_WORLD);
MPI_Bcast(&storage, 1, MPI_INT, 0, PETSC_COMM_WORLD);
unsigned flags =
MB_TAG_CREAT | storage;
std::vector<double> def_val(tag_length, 0.);
CHKERR moab.tag_get_handle(iterp_tag_name, tag_length, dtype, interp_tag_all,
flags, def_val.data());
MPI_Barrier(PETSC_COMM_WORLD);
auto vertex_exchange = CommInterface::createEntitiesPetscVector(
auto edge_exchange = CommInterface::createEntitiesPetscVector(
auto face_exchange = CommInterface::createEntitiesPetscVector(
auto volume_exchange = CommInterface::createEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR moab.delete_entities(&target_root, 1);
#endif
}
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);
}
}
}
}
}
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) {
"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) {
"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) {
"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;
<< "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) {
"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) {
"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;
MOFEM_LOG(
"EP", Sev::inform) <<
"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);
}
};
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);
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe->getRuleHook = [](
int,
int,
int) {
return -1; };
auto vol_rule_lin = [](int o) { return 2 * o + 1; };
auto vol_rule_no_lin = [](int o) { return 2 * o + (o - 1) + 1; };
auto vol_rule = (
SMALL_ROT > 0) ? vol_rule_lin : vol_rule_no_lin;
vol_rule, bubble_cache);
}
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;
};
constexpr bool stab_tau_dot_variant = false;
auto local_tau_sacale = boost::make_shared<double>(1.0);
using BdyEleOp = BoundaryEle::UserDataOperator;
set_scale_op->doWorkRhsHook =
[this,
auto op_ptr =
static_cast<BdyEleOp *
>(raw_op_ptr);
auto h = std::sqrt(op_ptr->getFTensor1Normal().l2());
};
auto not_interface_face = [
this](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (
) {
return false;
};
return true;
};
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",
}
auto ts_internal_stress =
boost::make_shared<DynamicRelaxationTimeScale>(
"internal_stress_history.txt", false, def_scaling_fun);
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_rhs = [&](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>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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_tau_stabilsation_rhs = [&](auto &pip, auto side_fe_name,
auto hybrid_field) {
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},
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
auto broken_disp_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getOpPtrVector().push_back(
broken_disp_data_ptr));
auto disp_mat_ptr = boost::make_shared<MatrixDouble>();
if (stab_tau_dot_variant) {
op_loop_domain_side->getOpPtrVector().push_back(
disp_mat_ptr));
} else {
op_loop_domain_side->getOpPtrVector().push_back(
disp_mat_ptr));
}
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
if (stab_tau_dot_variant) {
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_ptr));
} else {
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_ptr));
}
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_field, hybrid_ptr,
[local_tau_sacale, broken_disp_data_ptr](double, double, double) {
return broken_disp_data_ptr->size() * (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
broken_disp_data_ptr, [local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_field, broken_disp_data_ptr,
[local_tau_sacale](double, double, double) {
return -(*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
broken_disp_data_ptr, hybrid_ptr,
[local_tau_sacale](double, double, double) {
return -(*local_tau_sacale);
}));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact_rhs = [&](auto &pip) {
};
auto set_cohesive_rhs = [&](auto &pip) {
*this,
[](int p) { return 2 * (p + 1) + 1; }),
};
CHKERR set_hybridisation_rhs(fe_rhs->getOpPtrVector());
CHKERR set_contact_rhs(fe_rhs->getOpPtrVector());
}
CHKERR set_cohesive_rhs(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_lhs = [&](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>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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_tau_stabilsation_lhs = [&](auto &pip, auto side_fe_name,
auto hybrid_field) {
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},
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
auto broken_disp_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getOpPtrVector().push_back(
broken_disp_data_ptr));
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
auto time_shift = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a; };
hybrid_field, hybrid_field,
[local_tau_sacale, broken_disp_data_ptr](double, double, double) {
return broken_disp_data_ptr->size() * (*local_tau_sacale);
}));
if (stab_tau_dot_variant) {
op_loop_skeleton_side->getOpPtrVector().back())
.feScalingFun = time_shift;
}
op_loop_skeleton_side->getOpPtrVector().push_back(
broken_disp_data_ptr, [local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
if (stab_tau_dot_variant) {
op_loop_skeleton_side->getOpPtrVector().back())
.feScalingFun = time_shift;
}
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_field, broken_disp_data_ptr,
[local_tau_sacale](double, double, double) {
return -(*local_tau_sacale);
},
false, false));
if (stab_tau_dot_variant) {
op_loop_skeleton_side->getOpPtrVector().back())
.feScalingFun = time_shift;
}
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_field, broken_disp_data_ptr,
[local_tau_sacale](double, double, double) {
return -(*local_tau_sacale);
},
true, true));
if (stab_tau_dot_variant) {
op_loop_skeleton_side->getOpPtrVector().back())
.feScalingFun = time_shift;
}
pip.push_back(op_loop_skeleton_side);
};
auto set_contact_lhs = [&](auto &pip) {
};
auto set_cohesive_lhs = [&](auto &pip) {
*this,
[](int p) { return 2 * (p + 1) + 1; }),
};
CHKERR set_hybridisation_lhs(fe_lhs->getOpPtrVector());
CHKERR set_contact_lhs(fe_lhs->getOpPtrVector());
}
CHKERR set_cohesive_lhs(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; };
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
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>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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 flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
flux_mat_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
pip.push_back(op_loop_domain_side);
return std::make_tuple(broken_data_ptr, piola_scale_ptr);
};
auto set_rhs = [&]() {
auto [broken_data_ptr, piola_scale_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] =
get_broken_op_side(fe_lhs->getOpPtrVector());
auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
fe_lhs->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()));
};
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
}
boost::shared_ptr<ForcesAndSourcesCore> &fe_contact_tree
) {
}
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
#endif
auto setup_ts_monitor = [&]() {
boost::shared_ptr<TsCtx>
ts_ctx;
"get TS ctx");
if (set_ts_monitor) {
TSMonitorSet(ts, TsMonitorSet,
ts_ctx.get(), PETSC_NULLPTR),
"TS monitor set");
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;
"get DM section");
int num_fields;
"get num fields");
for (int ff = 0; ff != num_fields; ff++) {
PetscSectionGetFieldName(section_raw, ff, &
field_name),
"get field name");
}
};
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";
"append options prefix");
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);
<< "Debug model flag is " << (debug_model ? "ON" : "OFF");
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()) {
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
CHKERR TSElasticPostStep::postStepInitialise(
this);
CHKERR TSSolve(ts, PETSC_NULLPTR);
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
#ifndef NDEBUG
CHKERR PipelineGraph::writeTSGraphGraphviz(ts_ctx_ptr.get(),
"solve_elastic_graph.dot");
}
#endif
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);
}
}
}
int start_step,
double start_time) {
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();
<< "Dynamic relaxation final time -dynamic_final_time = " << final_time;
<< "Dynamic relaxation delta time -dynamic_delta_time = " << delta_time;
<< "Dynamic relaxation max iterations -dynamic_max_it = " << max_it;
<< "Dynamic relaxation H1 update each step -dynamic_h1_update = "
<< (ts_h1_update ? "TRUE" : "FALSE");
auto setup_ts_monitor = [&]() {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
return monitor_ptr;
};
auto monitor_ptr = setup_ts_monitor();
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSElasticPostStep::postStepInitialise(
this);
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
if (ts_h1_update) {
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TSElasticPostStep::preStepFun(ts);
CHKERR TSElasticPostStep::postStepFun(ts);
monitor_ptr->ts_u = PETSC_NULLPTR;
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
if (!ts_h1_update) {
CHKERR TSElasticPostStep::preStepFun(ts);
}
CHKERR TSSolve(ts, PETSC_NULLPTR);
if (!ts_h1_update) {
CHKERR TSElasticPostStep::postStepFun(ts);
}
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts_u = x;
break;
}
CHKERR TSElasticPostStep::postStepDestroy();
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;
<< "Block " << name << " id " << bc->getMeshsetId() << " has "
<<
r.size() <<
" entities";
}
};
return std::make_pair(name, map);
};
auto set_skin = [&](auto &&map) {
for (
auto &
m : map.second) {
<<
"Skin for block " << map.first <<
" id " <<
m.first <<
" has " <<
m.second.size() <<
" entities";
}
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) {
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));
}
{
auto get_crack_tag = [&]() {
if (rval == MB_SUCCESS) {
}
int def_val[] = {0};
"CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT,
def_val));
};
Tag th = get_crack_tag();
tags_to_transfer.push_back(th);
int mark[] = {1};
}
tags_to_transfer.push_back(
t);
}
}
auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
auto post_proc_ptr =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
auto domain_ops = [&](auto &fe, int sense) {
fe.getUserPolynomialBase() =
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
}
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 != 0) {
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[
"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(),
};
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) {
using SideEleOp = EleOnSide::UserDataOperator;
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
auto traction_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, traction_ptr, boost::make_shared<double>(1.0)));
pip.push_back(op_loop_domain_side);
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
&post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
op_this->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"ContactDisplacement", dataAtPts->getContactL2AtPts()}},
{},
{}
)
);
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);
auto skin_post_proc_ptr = get_post_proc(post_proc_mesh, 1);
CHKERR calcs_side_traction_and_displacements(
skin_post_proc_ptr, skin_post_proc_ptr->getOpPtrVector());
auto own_tets =
CommInterface::getPartEntities(mField.get_moab(), mField.get_comm_rank())
own_faces, moab::Interface::UNION);
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 side_one_faces = [&](auto &faces) {
std::pair<Range, Range> sides;
for (auto f : faces) {
MOAB_THROW(mField.get_moab().get_adjacencies(&f, 1, 3,
false, adj));
adj = intersect(own_tets, adj);
int side, sense, offset;
MOAB_THROW(mField.get_moab().side_number(
t, f, side, sense, offset));
if (sense == 1) {
sides.first.insert(f);
} else {
sides.second.insert(f);
}
}
}
return sides;
};
auto crack_faces =
unite(get_crack_faces(*crackFaces), get_crack_faces(*interfaceFaces));
auto crack_side_faces = side_one_faces(crack_faces);
auto side_one_crack_faces = [crack_side_faces](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_side_faces.first.find(ent) == crack_side_faces.first.end()) {
return false;
}
return true;
};
auto side_minus_crack_faces = [crack_side_faces](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_side_faces.second.find(ent) == crack_side_faces.second.end()) {
return false;
}
return true;
};
skin_post_proc_ptr->setTagsToTransfer(tags_to_transfer);
post_proc_ptr->setTagsToTransfer(tags_to_transfer);
post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
auto post_proc_begin =
post_proc_ptr->exeTestHook = side_one_crack_faces;
dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
post_proc_negative_sense_ptr->exeTestHook = side_minus_crack_faces;
post_proc_negative_sense_ptr, 0,
mField.get_comm_size());
constexpr bool debug =
false;
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;
};
auto only_front_faces = [adj_front](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (adj_front.find(ent) == adj_front.end()) {
return false;
}
return true;
};
post_proc_ptr->exeTestHook = only_front_faces;
dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
post_proc_negative_sense_ptr->exeTestHook = only_front_faces;
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>>(
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
auto hybrid_disp = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getOpPtrVector().push_back(
auto op_loop_domain_side =
post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
op_loop_domain_side->getOpPtrVector().push_back(
} else {
op_loop_domain_side->getOpPtrVector().push_back(
}
OpPPMap::DataMapMat vec_fields;
vec_fields["HybridDisplacement"] = hybrid_disp;
vec_fields[
"spatialL2Disp"] =
dataAtPts->getSmallWL2AtPts();
vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
OpPPMap::DataMapMat mat_fields;
mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
mat_fields["HybridDisplacementGradient"] =
OpPPMap::DataMapMat mat_fields_symm;
mat_fields_symm[
"LogSpatialStretch"] =
dataAtPts->getLogStretchTensorAtPts();
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
vec_fields,
mat_fields,
mat_fields_symm
)
);
if (f_residual) {
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(),
{},
{{"res_hybrid", hybrid_res}},
{},
{}
)
);
}
post_proc_ptr->setTagsToTransfer(tags_to_transfer);
auto post_proc_begin =
CHKERR post_proc_end.writeFile(file.c_str());
}
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() =
CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
}
}
&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 << " ) ";
};
}
int start_step,
double start_time) {
auto storage = solve_elastic_setup::setup(this, ts, x, false);
this,
[](int p) { return 2 * (p + 1) + 1; }),
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();
<< "Dynamic relaxation final time -dynamic_final_time = " << final_time;
<< "Dynamic relaxation delta time -dynamic_delta_time = " << delta_time;
<< "Dynamic relaxation max iterations -dynamic_max_it = " << max_it;
<< "Dynamic relaxation H1 update each step -dynamic_h1_update = "
<< (ts_h1_update ? "TRUE" : "FALSE");
auto setup_ts_monitor = [&]() {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
return monitor_ptr;
};
auto monitor_ptr = setup_ts_monitor();
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSElasticPostStep::postStepInitialise(
this);
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
if (ts_h1_update) {
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TSElasticPostStep::preStepFun(ts);
CHKERR TSElasticPostStep::postStepFun(ts);
CHKERR TaoSetType(tao, TAOLMVM);
auto g = cohesive_tao_ctx->duplicateGradientVec();
cohesiveEvaluateObjectiveAndGradient,
(void *)cohesive_tao_ctx.get());
monitor_ptr->ts_u = PETSC_NULLPTR;
auto tao_sol0 = cohesive_tao_ctx->duplicateKappaVec();
int tao_sol_size, tao_sol_loc_size;
CHKERR VecGetSize(tao_sol0, &tao_sol_size);
CHKERR VecGetLocalSize(tao_sol0, &tao_sol_loc_size);
<< "Cohesive crack growth initial kappa vector size " << tao_sol_size
<< " local size " << tao_sol_loc_size << " number of interface faces "
CHKERR TaoSetFromOptions(tao);
CHKERR VecSet(xu, PETSC_INFINITY);
CHKERR TaoSetVariableBounds(tao, xl, xu);
CHKERR VecZeroEntries(tao_sol0);
CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR TaoSetSolution(tao, tao_sol0);
CHKERR TaoGetSolution(tao, &tao_sol);
auto &kappa_vec = cohesive_tao_ctx->getKappaVec();
CHKERR VecAXPY(kappa_vec.second, 1.0, tao_sol);
CHKERR VecGhostUpdateBegin(kappa_vec.second, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(kappa_vec.second, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts_u = x;
break;
}
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
}
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Implementation of CGGUserPolynomialBase class.
Auxilary functions for Eshelbian plasticity.
Contains definition of EshelbianMonitor class.
FormsIntegrators< FaceElementForcesAndSourcesCore::UserDataOperator >::Assembly< A >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassVectorFace
FormsIntegrators< FaceElementForcesAndSourcesCore::UserDataOperator >::Assembly< A >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBaseTimesVectorFace
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.
Add cohesive elements to Eshelbian plasticity module.
#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)
static PetscErrorCode ierr
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#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_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)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#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
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark DOFs on block entities for boundary conditions.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
boost::shared_ptr< ContactSDFPython > setupContactSdf()
Read SDF file and setup contact SDF.
ForcesAndSourcesCore::UserDataOperator * getOpContactDetection(EshelbianCore &ep, boost::shared_ptr< ForcesAndSourcesCore > contact_tree_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr, boost::shared_ptr< MatrixDouble > contact_traction_ptr, Range r, moab::Interface *post_proc_mesh_ptr, std::vector< EntityHandle > *map_gauss_pts_ptr)
Push operator for contact detection.
boost::shared_ptr< ForcesAndSourcesCore > createContactDetectionFiniteElement(EshelbianCore &ep)
Create a Contact Tree finite element.
MoFEMErrorCode pushContactOpsRhs(EshelbianCore &ep, boost::shared_ptr< ForcesAndSourcesCore > contact_tree_ptr, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip)
Push contact operations to the right-hand side.
boost::shared_ptr< CohesiveTAOCtx > createCohesiveTAOCtx(EshelbianCore *ep, ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face, SmartPetscObj< TS > time_solver)
MoFEMErrorCode pushContactOpsLhs(EshelbianCore &ep, boost::shared_ptr< ForcesAndSourcesCore > contact_tree_ptr, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip)
Push contact operations to the left-hand side.
MoFEMErrorCode pushCohesiveOpsLhs(EshelbianCore &ep, ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face, boost::shared_ptr< Range > interface_range_ptr, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip)
Tag get_kappa_tag(moab::Interface &moab)
MoFEMErrorCode pushCohesiveOpsRhs(EshelbianCore &ep, ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face, boost::shared_ptr< Range > interface_range_ptr, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip)
MoFEMErrorCode initializeCohesiveKappaField(EshelbianCore &ep)
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.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
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.
PetscErrorCode TaoSetObjectiveAndGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
Sets the objective function value and gradient for a TAO optimization solver.
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 createTao(MPI_Comm comm)
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
constexpr auto field_name
#define QUAD_2D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_3D_TABLE[]
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
FTensor::Index< 'm', 3 > m
CGG User Polynomial Base.
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static enum StretchSelector stretchSelector
boost::shared_ptr< Range > frontAdjEdges
MoFEMErrorCode 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
static boost::function< double(const double)> inv_dd_f
MoFEM::Interface & mField
const std::string spatialL2Disp
static double inv_d_f_log(const double v)
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
static enum SolverType solverType
static PetscBool l2UserBaseScale
SmartPetscObj< DM > dM
Coupled problem all fields.
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
MoFEMErrorCode projectInternalStress(const EntityHandle meshset=0)
static int internalStressInterpOrder
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
const std::string materialH1Positions
static int nbJIntegralContours
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.
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
const std::string elementVolumeName
static double dd_f_log_e(const double v)
static enum RotSelector rotSelector
MoFEMErrorCode addDebugModel(TS ts)
Add debug to model.
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
static PetscBool dynamicRelaxation
const std::string spatialH1Disp
MoFEMErrorCode solveElastic(TS ts, Vec x)
static double d_f_log(const double v)
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
static double crackingStartTime
MoFEMErrorCode getOptions()
const std::string piolaStress
MoFEMErrorCode setElasticElementToTs(DM dm)
static double inv_d_f_log_e(const double v)
MoFEMErrorCode setFaceInterfaceOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
int contactRefinementLevels
MoFEMErrorCode gettingNorms()
[Getting norms]
boost::shared_ptr< Range > interfaceFaces
MoFEMErrorCode setVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
const std::string bubbleField
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
static double inv_f_log(const double v)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
static PetscBool noStretch
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
static double exponentBase
static double dd_f_linear(const double v)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
boost::shared_ptr< Range > skeletonFaces
boost::shared_ptr< PhysicalEquations > physicalEquations
const std::string rotAxis
static double inv_d_f_linear(const double v)
BitRefLevel bitAdjParentMask
bit ref level for parent parent
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x, int start_step, double start_time)
Solve problem using dynamic relaxation method.
static double inv_dd_f_log(const double v)
const std::string contactDisp
static std::string internalStressTagName
static enum SymmetrySelector symmetrySelector
CommInterface::EntitiesPetscVector edgeExchange
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
static boost::function< double(const double)> f
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
boost::shared_ptr< ForcesAndSourcesCore > contactTreeRhs
Make a contact tree.
const std::string skinElement
static PetscBool internalStressVoigt
static double inv_dd_f_linear(const double v)
static double inv_dd_f_log_e(const double v)
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
static PetscBool setSingularity
static double d_f_log_e(const double v)
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
static double f_log_e_quadratic(const double v)
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
MoFEMErrorCode solveCohesiveCrackGrowth(TS ts, Vec x, int start_step, double start_time)
Solve cohesive crack growth problem.
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
BitRefLevel bitAdjParent
bit ref level for parent
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ForcesAndSourcesCore > &fe_contact_tree)
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
const std::string naturalBcElement
static boost::function< double(const double)> dd_f
static double f_log_e(const double v)
static double inv_f_log_e(const double v)
MoFEMErrorCode createExchangeVectors(Sev sev)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< Range > crackFaces
static boost::function< double(const double)> d_f
static PetscBool intefaceCrack
boost::shared_ptr< Range > frontVertices
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_d_f
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
static double d_f_linear(const double v)
const std::string hybridSpatialDisp
SmartPetscObj< Vec > solTSStep
static double f_log(const double v)
CommInterface::EntitiesPetscVector faceExchange
SmartPetscObj< DM > dmElastic
Elastic problem.
EshelbianCore(MoFEM::Interface &m_field)
boost::shared_ptr< Range > frontEdges
static boost::function< double(const double)> inv_f
const std::string stretchTensor
BitRefLevel bitAdjEntMask
bit ref level for parent parent
static double f_linear(const double v)
const std::string contactElement
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
boost::function< int(int)> FunRule
boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cachePhi
SetIntegrationAtFrontVolume(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges, boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cache_phi=nullptr)
boost::shared_ptr< Range > frontEdges
Set integration rule to boundary elements.
int operator()(int, int, int) const
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Boundary condition manager for finite element problem setup.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)> > > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
Field data structure for finite element approximation.
Definition of the force bc data structure.
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Constructor for operators working on finite element spaces.
Section manager is used to create indexes and sections.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Operator for broken loop side.
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients time derivative at integration pts for scalar field rank 0, i....
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Execute "this" element in the operator.
Post post-proc data at points from hash maps.
Template struct for dimension-specific finite element types.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
std::function< double(double)> ScalingFun
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.
IFACE getInterface() const
Get interface pointer to pointer of interface.
default operator for TET element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Apply rotation boundary condition.
Set integration rule to volume elements.
int operator()(int, int, int) const
BoundaryEle::UserDataOperator BdyEleOp
ElementsAndOps< SPACE_DIM >::SideEle SideEle
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
