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;
#endif
extern "C" {
}
#include <queue>
struct VolUserDataOperatorStabAssembly
using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
}
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) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
auto meshset = bc->getMeshset();
"get meshset ents");
}
};
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_release[] = {"griffith_force", "griffith_skeleton"};
const char *list_stretches[] = {"linear", "log", "log_quadratic"};
const char *list_broken_hdiv_bases[] = {"demkowicz", "ainsworth"};
PetscInt choice_stretch = StretchSelector::LOG;
PetscInt choice_broken_hdiv_base = 0;
PetscBool l2_user_base_scale_set = PETSC_FALSE;
choice_broken_hdiv_base = 0;
break;
choice_broken_hdiv_base = 1;
break;
default:
"Unsupported broken HDIV base %s",
}
char analytical_expr_file_name[255] = "analytical_expr.py";
PetscBool no_stretch =
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(
"-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",
"",
no_stretch, &no_stretch, PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
&l2_user_base_scale_set);
"-broken_hdiv_base", "broken HDIV stress approximation base", "",
list_broken_hdiv_bases, 2,
list_broken_hdiv_bases[choice_broken_hdiv_base],
&choice_broken_hdiv_base, PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
CHKERR PetscOptionsScalar(
"-physical_final_time",
"physical final time",
"",
PETSC_NULLPTR);
CHKERR PetscOptionsScalar(
"-physical_delta_time",
"physical delta time",
"",
PETSC_NULLPTR);
CHKERR PetscOptionsInt(
"-physical_max_steps",
"physical max iterations",
"",
PETSC_NULLPTR);
"-physical_h1_update", "update each physicalsolver step", "",
}
CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
PETSC_NULLPTR);
CHKERR PetscOptionsBool(
"-cohesive_interface_on",
"cohesive interface ON",
"",
"-cohesive_interface_remove_level", "cohesive interface remove level", "",
CHKERR PetscOptionsScalar(
"-cracking_add_time",
"cracking add time",
"",
CHKERR PetscOptionsScalar(
"-cracking_start_time",
"cracking start time",
"",
PETSC_NULLPTR);
CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
CHKERR PetscOptionsScalar(
"-cracking_rtol",
"Cracking relative tolerance",
"",
CHKERR PetscOptionsScalar(
"-cracking_atol",
"Cracking absolute tolerance",
"",
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);
CHKERR PetscOptionsBool(
"-internal_stress_voigt",
"Voigt index notation",
"",
PETSC_NULLPTR);
char tag_heterogeneous_youngs_modulus_name[255] = "";
"-heterogeneous_youngs_modulus", "heterogeneous Young's modulus tag name",
"", "", tag_heterogeneous_youngs_modulus_name, 255, PETSC_NULLPTR);
"-analytical_expr_file",
analytical_expr_file_name, 255, PETSC_NULLPTR);
PetscOptionsEnd();
PetscOptionsBegin(PETSC_COMM_WORLD, "mesh_transfer_", "mesh data transfer",
"none");
char tag_mesh_transfer_source_file_name[255] = "";
CHKERR PetscOptionsString(
"-source_file",
"source mesh file name",
"",
"source.h5m", tag_mesh_transfer_source_file_name,
CHKERR PetscOptionsInt(
"-interp_order",
"interpolation order",
"", 0,
CHKERR PetscOptionsBool(
"-hybrid_interp",
"use hybrid interpolation",
"",
PETSC_NULLPTR);
PetscOptionsEnd();
"Unsupported mesh transfer interpolation order %d",
}
}
switch (choice_broken_hdiv_base) {
case 0:
break;
case 1:
break;
default:
"Unknown broken HDIV base option");
}
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;
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];
<< "No stretch: -no_stretch "
? "yes"
: "no";
MOFEM_LOG(
"EP", Sev::inform) <<
"Solver type: -solver_type "
? "yes"
: "no";
? "yes"
: "no";
<< "Broken HDIV base: -broken_hdiv_base "
<< list_broken_hdiv_bases[choice_broken_hdiv_base];
? "yes"
: "no";
<<
"Cracking relative tolerance: -cracking_rtol " <<
crackingRtol;
<<
"Cracking absolute tolerance: -cracking_atol " <<
crackingAtol;
<< "Energy release variant: -energy_release_variant "
<< "Number of J integral contours: -nb_J_integral_contours "
<< "Cohesive interface on: -cohesive_interface_on "
<< "Cohesive interface remove level: -cohesive_interface_remove_level "
#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
}
const bool add_bubble) {
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, broken_hdiv_base](
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;
};
if (add_bubble) {
}
: -1,
6);
get_hybridised_disp());
};
}
double time) {
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));
};
} else {
}
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);
};
auto get_interface_from_block = [&](auto block_name) {
faces, moab::Interface::UNION);
faces = subtract(faces, skin);
<< "Number of vol interface elements: " << vol_eles.size()
<< " and faces: " << faces.size();
return faces;
};
auto remove_interface_from_block = [&](auto block_name, auto level) {
for (
auto l = 0;
l < level; ++
l) {
ents,
SPACE_DIM,
true, adj_tets, moab::Interface::UNION);
adj_tets_faces,
moab::Interface::UNION);
ents.merge(adj_tets_faces);
}
auto faces = ents.subset_by_dimension(
SPACE_DIM - 1);
if (faces.size()) {
<< "Removed ents " << faces.size()
}
<< "Interface faces after remove " << intreface_faces;
}
};
}
#ifdef INCLUDE_MBCOUPLER
double toler = 5.e-10;
<< "No source mesh specified. Skipping projection";
}
MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Interpolation Young's modulus tag name: "
MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Using hybrid interpolation: "
auto rval_check_tag = moab.tag_get_handle(tag_name.c_str(), old_interp_tag);
if (rval_check_tag == MB_SUCCESS) {
<< "Deleting existing tag on target mesh: " << 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;
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);
std::vector<Tag> interp_tags;
std::vector<int> tag_length;
std::vector<DataType> dtype;
std::vector<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(
if (rval_source_mesh != MB_SUCCESS) {
MOFEM_LOG(
"WORLD", Sev::warning) <<
"Error loading source mesh 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);
std::string tag_to_use = iterp_tag_name;
CHKERR moab.tag_get_handle(tag_to_use.c_str(), 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.push_back(interp_tag_len);
dtype.push_back(DataType());
storage.push_back(TagType());
interp_tags.push_back(interp_tag);
CHKERR moab.tag_get_data_type(interp_tag, dtype.back());
CHKERR moab.tag_get_type(interp_tag, storage.back());
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);
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(tag_to_use) + "_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]);
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]);
}
};
<< "Performing interpolation for tag: " << tag_to_use;
<< "Number of target points to interpolate: " << num_pts;
<< "Interpolation method: "
<< (method == Coupler::CONSTANT ? "constant" : "linear FE");
<< "Number of components in tag: " << interp_tag_len;
<< "Source tag data range: ["
<< *std::min_element(source_data.begin(), source_data.end()) << ", "
<< *std::max_element(source_data.begin(), source_data.end()) << "]";
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]);
<< "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);
}
int tag_size = tag_length.size();
MPI_Bcast(&tag_size, 1, MPI_INT, 0, PETSC_COMM_WORLD);
if (rank != 0) {
interp_tags.resize(tag_size);
tag_length.resize(tag_size);
dtype.resize(tag_size);
storage.resize(tag_size);
}
MPI_Bcast(interp_tags.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
MPI_Bcast(tag_length.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
MPI_Bcast(dtype.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
MPI_Bcast(storage.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
for (size_t index = 0; index < interp_tags.size(); index++) {
if (world_rank) {
auto rval_check_tag =
if (rval_check_tag == MB_SUCCESS) {
<< "Deleting existing tag on target mesh (post-projection): "
CHKERR moab.tag_delete(old_interp_tag);
}
}
unsigned flags =
MB_TAG_CREAT | storage[index];
std::vector<double> def_val(tag_length[index], 0.);
tag_length[index], dtype[index],
interp_tag_all, flags, def_val.data());
if (rval != MB_SUCCESS && world_rank) {
"Unable to create projection tag %s",
}
MPI_Barrier(PETSC_COMM_WORLD);
auto vertex_exchange = CommInterface::createEntitiesPetscVector(
Sev::inform);
auto volume_exchange = CommInterface::createEntitiesPetscVector(
Sev::inform);
CHKERR CommInterface::updateEntitiesPetscVector(
CHKERR CommInterface::updateEntitiesPetscVector(
}
CHKERR moab.delete_entities(&target_root, 1);
#endif
}
const bool add_bubble) {
auto add_field_to_fe = [this](const std::string fe,
};
if (add_bubble) {
}
}
}
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 zero_kinetic_constraints_block = [&]() {
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"REMOVE_Y",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"REMOVE_Z",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"REMOVE_ALL",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"FIX_X",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"FIX_Y",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"FIX_Z",
CHKERR bc_mng->removeBlockDOFsOnEntities(
"MATERIAL_PROBLEM",
"FIX_ALL",
};
};
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 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,
const bool add_bubble) {
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; };
vol_rule, bubble_cache);
}
if (add_bubble) {
}
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()) {
}
}
} else {
}
}
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
bool has_nonhomogeneous_mat_block =
}
}
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
fe_lhs->getOpPtrVector().push_back(
}
}
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
fe_lhs->getOpPtrVector().push_back(
}
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
fe_lhs->getOpPtrVector().push_back(
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto local_tau_sacale = boost::make_shared<double>(1.0);
using BdyEleOp = BoundaryEle::UserDataOperator;
struct OpSetTauScale :
public BdyEleOp {
OpSetTauScale(boost::shared_ptr<double> local_tau_sacale, double alphaTau)
localTauSacale(local_tau_sacale), alphaTau(alphaTau) {}
auto &coords = BdyEleOp::getCoords();
auto [centre, barycenter,
h] =
Tools::getTricircumcenter3d(coords.data().data());
*localTauSacale = (alphaTau /
h);
}
private:
boost::shared_ptr<double> localTauSacale;
double alphaTau;
};
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);
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 mesh transfer interpolation order %d, for "
"internal stress",
}
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>(nullptr, true);
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 =
op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
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>(nullptr, true);
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>();
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(
new OpSetTauScale(local_tau_sacale,
alphaTau));
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
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_tau_stabilsation_disp_bc_rhs = [&](auto &pip, auto side_fe_name) {
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 =
op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
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>(nullptr, true);
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>();
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(
new OpSetTauScale(local_tau_sacale,
alphaTau));
op_loop_skeleton_side->getOpPtrVector().push_back(
[local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
return (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
[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) {
};
CHKERR set_hybridisation_rhs(fe_rhs->getOpPtrVector());
CHKERR set_contact_rhs(fe_rhs->getOpPtrVector());
CHKERR set_tau_stabilsation_disp_bc_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);
if (add_elastic) {
} else {
}
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>(nullptr, true);
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 =
op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
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>(nullptr, true);
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(
new OpSetTauScale(local_tau_sacale,
alphaTau));
hybrid_field, hybrid_field,
[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);
},
false, false));
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));
pip.push_back(op_loop_skeleton_side);
};
auto set_tau_stabilsation_disp_bc_lhs = [&](auto &pip, auto side_fe_name) {
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 =
op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
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>(nullptr, true);
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(
new OpSetTauScale(local_tau_sacale,
alphaTau));
op_loop_skeleton_side->getOpPtrVector().push_back(
[local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
[local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
op_loop_skeleton_side->getOpPtrVector().push_back(
[local_tau_sacale](double, double, double) {
return (*local_tau_sacale);
}));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact_lhs = [&](auto &pip) {
};
auto set_cohesive_lhs = [&](auto &pip) {
};
CHKERR set_hybridisation_lhs(fe_lhs->getOpPtrVector());
CHKERR set_contact_lhs(fe_lhs->getOpPtrVector());
CHKERR set_tau_stabilsation_disp_bc_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>(nullptr, true);
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(
pip.push_back(op_loop_domain_side);
return broken_data_ptr;
};
auto set_rhs = [&]() {
auto broken_data_ptr = get_broken_op_side(fe_rhs->getOpPtrVector());
fe_rhs->getOpPtrVector().push_back(
auto piola_scale_ptr = boost::make_shared<double>(1.0);
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 = 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);
auto testing_monitor_ptr =
boost::make_shared<EshelbianTestingMonitor>(*ep_ptr, monitor_ptr);
PetscBool test_cook_flg = PETSC_FALSE;
PetscBool test_cook_pts_flg = PETSC_FALSE;
&test_cook_flg, PETSC_NULLPTR);
&test_cook_pts_flg, PETSC_NULLPTR);
PETSC_NULLPTR);
if (
atom_test || test_cook_flg || test_cook_pts_flg) {
}
}
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);
PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
CHKERR PetscOptionsScalar(
"-dynamic_final_time",
"dynamic relaxation final time", "",
CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
PetscOptionsEnd();
<< "Following options are deprecated, use -physical prefix options "
"instead";
<< "Dynamic relaxation final time -dynamic_final_time = "
<< "Dynamic relaxation delta time -dynamic_delta_time = "
<< "Dynamic relaxation H1 update each step -dynamic_h1_update = "
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);
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TSElasticPostStep::preStepFun(ts);
CHKERR TSElasticPostStep::postStepFun(ts);
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = PETSC_NULLPTR;
"physicalDt must be positive, got %g",
physicalDt);
}
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
CHKERR TSElasticPostStep::preStepFun(ts);
}
CHKERR TSSolve(ts, PETSC_NULLPTR);
CHKERR TSElasticPostStep::postStepFun(ts);
}
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = x;
break;
const double remainingPhysicalTime =
} else {
}
}
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))));
}
std::vector<Tag> tags_to_transfer) {
ParallelComm *pcomm =
std::vector<boost::shared_ptr<TempMeshset>> meshsets_tmp_list;
std::vector<Tag> tags_list;
auto meshset =
m.getMeshset();
std::vector<Tag> tmp_tags_list;
for (
auto t : tmp_tags_list) {
void *tag_vals[1];
int tag_size[1];
t, &meshset, 1, (
const void **)tag_vals, tag_size);
tag_size);
}
std::vector<std::string> remove_tags;
remove_tags.push_back("AKDTree_coord_norm");
remove_tags.push_back("__PARALLEL_");
remove_tags.push_back("_RefBitLevel");
for (
auto t : tmp_tags_list) {
std::string tag_name;
bool add = true;
for (auto &p : remove_tags) {
if (tag_name.compare(0, p.size(), p) == 0) {
add = false;
break;
}
}
if (add)
}
}
for (auto &m_ptr : meshsets_tmp_list) {
}
std::sort(tags_list.begin(), tags_list.end());
auto new_end = std::unique(tags_list.begin(), tags_list.end());
tags_list.resize(std::distance(tags_list.begin(), new_end));
1, &tags_list[0], tags_list.size());
}
}
Vec f_residual, Vec var_vector, Vec gradient,
std::vector<Tag> tags_to_transfer, TS ts) {
if (f_residual != PETSC_NULLPTR || var_vector != PETSC_NULLPTR) {
auto xin = f_residual != PETSC_NULLPTR ? f_residual : var_vector;
if (f_residual) {
CHKERR VecScatterBegin(scatter, f_residual, f_r, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, f_residual, f_r, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(f_r, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(f_r, INSERT_VALUES, SCATTER_FORWARD);
}
if (var_vector) {
CHKERR VecScatterBegin(scatter, var_vector, v_v, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, var_vector, v_v, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(v_v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(v_v, INSERT_VALUES, SCATTER_FORWARD);
}
}
if (gradient) {
CHKERR VecScatterBegin(scatter, gradient,
g, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, gradient,
g, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
}
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};
}
std::string name;
<< "Adding tag " << name << " to transfer list for post-processing";
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(
if (ts != PETSC_NULLPTR) {
post_proc_ptr->data_ctx |= PetscData::CTX_SET_TIME;
CHKERR TSGetTime(ts, &(post_proc_ptr->ts_t));
CHKERR TSGetTimeStep(ts, &(post_proc_ptr->ts_dt));
}
auto domain_ops = [&](auto &fe, int sense) {
auto bubble_cache = boost::make_shared<CGGUserPolynomialBase::CachePhi>(
fe.getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
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), v_v));
boost::make_shared<double>(1), v_v, MBMAXTYPE));
fe.getOpPtrVector().push_back(
} else {
fe.getOpPtrVector().push_back(
}
}
if (gradient) {
}
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())
}
}
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->getEigenVals();
}
if (var_vector) {
vec_fields[
"VarOmega"] =
dataAtPts->getVarRotAxisPts();
vec_fields["VarSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], v_v, MBTET));
}
if (f_residual) {
vec_fields["ResSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], f_r, MBTET));
vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
rotAxis, vec_fields[
"ResOmega"], f_r, MBTET));
}
if (gradient) {
vec_fields[
"Gradient"] =
dataAtPts->getGradientAtPts();
}
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), f_r));
boost::make_shared<double>(1), f_r, MBMAXTYPE));
}
case 0:
fe.getOpPtrVector().push_back(
break;
case 1:
fe.getOpPtrVector().push_back(
break;
default:
"Unsupported mesh transfer interpolation order %d, for "
"internal stress",
}
}
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() =
boost::shared_ptr<BaseFunction>(
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()));
pip.push_back(op_this);
op_this->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"ContactDisplacement", dataAtPts->getContactL2AtPts()}},
{},
{}
)
);
if (f_residual) {
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 get_interface_from_block = [&](auto block_name) {
faces, moab::Interface::UNION);
faces = subtract(faces, skin);
return faces;
};
auto crack_faces = unite(get_crack_faces(*crackFaces), *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());
}
}
Vec f_residual,
std::vector<Tag> tags_to_transfer) {
if (f_residual != PETSC_NULLPTR) {
CHKERR VecScatterBegin(scatter, f_residual, f_r, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter, f_residual, f_r, INSERT_VALUES,
SCATTER_FORWARD);
}
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>(nullptr, true);
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(
} 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(
f_r));
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 =
}
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
auto bubble_cache =
boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0,
MatrixDouble());
post_proc_norm_fe->getUserPolynomialBase() =
post_proc_norm_fe->getRuleHook = [](
int,
int,
int) {
return -1; };
post_proc_norm_fe->setRuleHook = SetIntegrationAtFrontVolume(
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);
PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
CHKERR PetscOptionsScalar(
"-dynamic_final_time",
"dynamic relaxation final time", "",
CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
PetscOptionsEnd();
<< "Dynamic relaxation final time -dynamic_final_time = "
<< "Dynamic relaxation delta time -dynamic_delta_time = "
<< "Dynamic relaxation H1 update each step -dynamic_h1_update = "
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);
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TaoSetType(tao, TAOLMVM);
auto g = cohesive_tao_ctx->duplicateGradientVec();
cohesiveEvaluateObjectiveAndGradient,
(void *)cohesive_tao_ctx.get());
monitor_ptr->ts = PETSC_NULLPTR;
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);
"physicalDt must be positive, got %g",
physicalDt);
}
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 TSElasticPostStep::preStepFun(ts);
}
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 = PETSC_NULLPTR;
monitor_ptr->ts_u = x;
CHKERR TSElasticPostStep::postStepFun(ts);
}
break;
const double remainingPhysicalTime =
} else {
}
}
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
double start_time) {
auto storage = solve_elastic_setup::setup(this, ts, x, false);
auto setup_ts_monitor = [&]() {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
return monitor_ptr;
};
auto monitor_ptr = setup_ts_monitor();
auto test_monitor_ptr =
boost::make_shared<EshelbianTestingMonitor>(*this, monitor_ptr);
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSElasticPostStep::postStepInitialise(
this);
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TSElasticPostStep::preStepFun(ts);
CHKERR TSElasticPostStep::postStepFun(ts);
double load_factor_change_clip = 0.1;
PetscOptionsBegin(PETSC_COMM_WORLD, "", "Load Factor Options", "none");
CHKERR PetscOptionsScalar(
"-initial_load_factor",
"Initial load factor",
"",
CHKERR PetscOptionsScalar(
"-max_crack_ext_area",
"Maximum crack extension area", "",
PETSC_NULLPTR);
CHKERR PetscOptionsScalar(
"-clip_load_factor_percent",
"Upper bound for load factor change", "",
load_factor_change_clip,
&load_factor_change_clip, PETSC_NULLPTR);
PetscOptionsEnd();
monitor_ptr->ts = ts;
monitor_ptr->ts_u = PETSC_NULLPTR;
PetscBool test_cook_flg = PETSC_FALSE;
PETSC_NULLPTR);
PETSC_NULLPTR);
test_monitor_ptr->ts = ts;
test_monitor_ptr->ts_u = PETSC_NULLPTR;
test_monitor_ptr);
}
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
CHKERR TSElasticPostStep::preStepFun(ts);
}
CHKERR TSSolve(ts, PETSC_NULLPTR);
CHKERR TSElasticPostStep::postStepFun(ts);
}
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts = ts;
monitor_ptr->ts_u = x;
test_monitor_ptr->ts = ts;
test_monitor_ptr->ts_u = x;
test_monitor_ptr);
}
const bool has_crack_extension = delta_area > 0.0;
if (has_crack_extension) {
if (denom > 0.0) {
const double updated_load_factor =
loadFactor = std::max(updated_load_factor, 1.0e-6);
}
}
const double initial_step_range = 5;
const double min_load_factor = 1.0e-6;
const double max_load_factor =
<< "Allowable range for load factor [" << min_load_factor << ", "
<< max_load_factor << "]";
}
const double previous_load_factor = is_first_step ? 0. :
oldLoadFactor;
}
CHKERR MPI_Bcast(load_control_data, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
break;
const double remainingPhysicalTime =
} else {
}
}
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
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);
CHKERR TaoSetType(tao, TAOLMVM);
topologicalEvaluateObjectiveAndGradient,
(void *)topological_tao_ctx.get());
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = PETSC_NULLPTR;
CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
int tao_sol_size, tao_sol_loc_size;
CHKERR VecGetSize(tao_sol0, &tao_sol_size);
CHKERR VecGetLocalSize(tao_sol0, &tao_sol_loc_size);
<< "Toplogical data vector size " << tao_sol_size << " local size "
<< tao_sol_loc_size << " number of interface faces "
CHKERR TaoSetFromOptions(tao);
if (delta_time <= 0.) {
"delta_time must be positive, got %g", delta_time);
}
<< " delta time " << delta_time;
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);
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = x;
break;
if (delta_time >= remainingPhysicalTime) {
} else {
}
}
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
int start_step,
double start_time) {
auto storage = solve_elastic_setup::setup(this, ts, x, false);
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
TetPolynomialBase::switchCacheBaseOn<HDIV>(
CHKERR TSElasticPostStep::postStepInitialise(
this);
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
}
CHKERR TSElasticPostStep::preStepFun(ts);
CHKERR TSElasticPostStep::postStepFun(ts);
const bool restart_run =
start_step != 0 ||
std::abs(start_time) > std::numeric_limits<double>::epsilon();
if (!restart_run &&
std::abs(test_time) < std::numeric_limits<double>::epsilon()) {
"Set non-zero -physical_final_time for test_topological_derivative");
}
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = PETSC_NULLPTR;
<< "Solving load step before topological derivative test: "
<< " TS delta time " << ts_delta_time;
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
CHKERR TSElasticPostStep::preStepFun(ts);
}
CHKERR TSSolve(ts, PETSC_NULLPTR);
CHKERR TSElasticPostStep::postStepFun(ts);
}
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
monitor_ptr->ts = PETSC_NULLPTR;
monitor_ptr->ts_u = x;
CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
int tao_sol_size, tao_sol_loc_size;
CHKERR VecGetSize(tao_sol0, &tao_sol_size);
CHKERR VecGetLocalSize(tao_sol0, &tao_sol_loc_size);
<< "Topological data vector size " << tao_sol_size << " local size "
<< tao_sol_loc_size << " number of interface faces "
const char *list_objective_models[ObjectiveModelType::LAST_MODEL] = {
"python_model", "hencky_model"};
#ifdef ENABLE_PYTHON_BINDING
PetscInt choice_objective_model = ObjectiveModelType::PYTHON_MODEL;
#else
PetscInt choice_objective_model = ObjectiveModelType::HENCKY_MODEL;
#endif
PETSC_NULLPTR, PETSC_NULLPTR, "-objective_model_type",
list_objective_models, ObjectiveModelType::LAST_MODEL,
&choice_objective_model, PETSC_NULLPTR);
const auto objective_model_type =
<< "Objective model type: -objective_model_type "
<< list_objective_models[objective_model_type];
PetscReal obj_value;
&obj_value,
g, objective_model_type);
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< VolUserDataOperator >::Assembly< A >::BiLinearForm< GAUSS >::OpMass< 9, 9 > OpStressGram_dBubble_dBubble
FormsIntegrators< FaceElementForcesAndSourcesCore::UserDataOperator >::Assembly< A >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, SPACE_DIM > OpBaseTimesVectorFace
FormsIntegrators< VolUserDataOperator >::Assembly< A >::BiLinearForm< GAUSS >::OpMass< 3, 9 > OpStressGram_dP_dP
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_entities_by_handle(MoFEM::Interface &m_field, const std::string block_name)
static auto get_crack_front_edges(MoFEM::Interface &m_field, Range crack_faces)
Eshelbian plasticity interface.
Contains definition of EshelbianTestingMonitor class.
#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
static const char *const ApproximationBaseNames[]
#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 DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from 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
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
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 removeBlockDOFsOnEntities(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, bool is_distributed_mesh=true)
Remove DOFs from problem based on block entities.
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< 'l', 3 > l
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
boost::shared_ptr< ContactSDFPython > setupContactSdf()
Read SDF file and setup contact SDF.
boost::shared_ptr< TopologicalTAOCtx > createTopologicalTAOCtx(EshelbianCore *ep, ForcesAndSourcesCore::GaussHookFun set_integration_at_interior, ForcesAndSourcesCore::GaussHookFun set_integration_at_face, SmartPetscObj< TS > time_solver)
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)
MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr, Vec sol, PetscReal *f, Vec g, ObjectiveModelType eval_energy_model)
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)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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 PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
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
static PetscBool physicalH1Update
static enum StretchSelector stretchSelector
boost::shared_ptr< Range > frontAdjEdges
static int interfaceRemoveLevel
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
const std::string skeletonElement
static double inv_f_linear(const double v)
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
boost::shared_ptr< Range > contactFaces
static double dd_f_log_e_quadratic(const double v)
static double dd_f_log(const double v)
BitRefLevel bitAdjEnt
bit ref level for parent
static boost::function< double(const double)> inv_dd_f
MoFEM::Interface & mField
const std::string spatialL2Disp
static double inv_d_f_log(const double v)
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
static enum SolverType solverType
static PetscBool l2UserBaseScale
SmartPetscObj< DM > dM
Coupled problem all fields.
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
MoFEMErrorCode solveSchapeOptimisation(TS ts, Vec x, int start_step, double start_time)
Solve shape optimisation problem.
static enum StretchHandling stretchHandling
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
static const char * listSolvers[]
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.
static double griffithEnergy
Griffith energy.
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
MoFEMErrorCode postProcessRestartMesh(const int tag, const std::string file, std::vector< Tag > tags_to_transfer={})
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
PetscBool loadFactorTSSolveExecuted
MoFEMErrorCode postProcessResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, Vec var_vec=PETSC_NULLPTR, Vec gradient=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={}, TS ts=PETSC_NULLPTR)
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
CommInterface::EntitiesPetscVector vertexExchange
static std::vector< std::string > listTagsToProject
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
static std::string heterogeneousYoungModTagName
const std::string spatialH1Disp
static FieldApproximationBase brokenHdivBase
static double maxCrackExtension
static int physicalMaxSteps
MoFEMErrorCode solveElastic(TS ts, Vec x)
@ TestTopologicalDerivative
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
static int physicalStepNumber
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)
static PetscBool physicalTimeFlg
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
SmartPetscObj< DM > dmMaterial
Material problem.
static double inv_f_log(const double v)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
static double crackingAddTime
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 projectGeometry(const EntityHandle meshset=0, double time=0)
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
static double currentPhysicalTime
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
static constexpr enum SymmetrySelector symmetrySelector
static double crackingAtol
Cracking absolute tolerance.
MoFEMErrorCode projectMaterialTags(const EntityHandle meshset=0)
boost::shared_ptr< Range > skeletonFaces
static double crackingRtol
Cracking relative tolerance.
boost::shared_ptr< PhysicalEquations > physicalEquations
const std::string rotAxis
static double inv_d_f_linear(const double v)
static PetscBool meshTransferHybridInterp
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
CommInterface::EntitiesPetscVector edgeExchange
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
static boost::function< double(const double)> f
MoFEMErrorCode solveTestTopologicalDerivative(TS ts, Vec x, int start_step, double start_time)
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
static double finalPhysicalTime
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)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0, const bool add_bubble=true)
static double inv_dd_f_log_e(const double v)
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
MoFEMErrorCode pushNoStretchVolumeA00Ops(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_lhs)
static PetscBool setSingularity
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe, const bool add_bubble=true)
static double d_f_log_e(const double v)
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
boost::shared_ptr< double > currentCrackAreaPtr
static PetscBool meshTransferSourceMeshFileSpecified
static double f_log_e_quadratic(const double v)
double avgGriffithsEnergy
static bool hasNonHomogeneousMaterialBlock
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)
static PetscBool interfaceCrack
MoFEMErrorCode solveLoadFactor(TS ts, Vec x, int start_step, double start_time)
Solve load factor crack growth problem.
static double d_f_log_e_quadratic(const double v)
CommInterface::EntitiesPetscVector volumeExchange
const std::string naturalBcElement
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts, SmartPetscObj< Vec > *adjoint_gradient_vector=nullptr)
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)
MoFEMErrorCode pushStretchVolumeA00Ops(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_lhs)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< Range > crackFaces
static boost::function< double(const double)> d_f
boost::shared_ptr< Range > frontVertices
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_d_f
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
static double d_f_linear(const double v)
static int meshTransferInterpOrder
const std::string hybridSpatialDisp
SmartPetscObj< Vec > solTSStep
static double f_log(const double v)
CommInterface::EntitiesPetscVector faceExchange
SmartPetscObj< DM > dmElastic
Elastic problem.
static std::string meshTransferSourceMeshFileName
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)
MoFEMErrorCode addFields(const EntityHandle meshset=0, const bool add_bubble=true)
MoFEMErrorCode withFieldOrders(Op &&op) const
MoFEMErrorCode pushStressGramOps(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_lhs)
const std::string contactElement
MoFEMErrorCode pushPiolaStressGramOps(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_lhs)
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
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.
Section manager is used to create indexes and sections.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
CubitMeshSet_multiIndex & getMeshsetsMultindex()
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 MatrixDouble 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.
Vector manager is used to create vectors \mofem_vectors.
default operator for TET element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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
