#define SINGULARITY
#include <boost/math/constants/constants.hpp>
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#pragma message "With PYTHON_SDF"
#else
#pragma message "Without PYTHON_SDF"
#endif
extern "C" {
}
struct VolUserDataOperatorStabAssembly
};
struct FaceUserDataOperatorStabAssembly
};
}
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();
sendbuf.push_back(id);
}
}
}
int recv_data;
MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
pcomm->comm());
std::vector<int> recvbuf(recv_data);
MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
if (pcomm->rank() > 0) {
for (
auto &
f : recvbuf) {
}
}
}
const std::string block_name, int dim) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
faces, true),
"get meshset ents");
}
};
const unsigned int cubit_bc_type) {
CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
return meshset;
};
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 =
if (pcomm->rank() == 0) {
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 = subtract(crack_skin, body_skin_edges);
}
return send_type(m_field, crack_skin, MBEDGE);
}
ParallelComm *pcomm =
MOFEM_LOG(
"EP", Sev::inform) <<
"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::inform) <<
"get_two_sides_of_crack_surface <- done";
return std::pair<Range, Range>();
}
struct DynamicRelaxationTimeScale :
public TimeScale {
double getScale(const double time) override {
} else {
}
}
};
struct SetIntegrationAtFrontVolume {
SetIntegrationAtFrontVolume(boost::shared_ptr<Range> front_nodes,
boost::shared_ptr<Range> front_edges)
: frontNodes(front_nodes), frontEdges(front_edges) {};
int order_col, int order_data) {
constexpr
bool debug =
false;
constexpr int numNodes = 4;
constexpr int numEdges = 6;
constexpr int refinementLevels = 4;
auto &m_field = fe_raw_ptr->
mField;
auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
auto fe_handle = fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
int rule = 2 * order_data + 1;
"wrong dimension");
}
}
auto &gauss_pts = fe_ptr->gaussPts;
gauss_pts.resize(4, nb_gauss_pts, false);
&gauss_pts(0, 0), 1);
&gauss_pts(1, 0), 1);
&gauss_pts(2, 0), 1);
&gauss_pts(3, 0), 1);
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
false);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
1);
} else {
}
};
auto get_singular_nodes = [&]() {
int num_nodes;
CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
true);
std::bitset<numNodes> singular_nodes;
for (auto nn = 0; nn != numNodes; ++nn) {
if (frontNodes->find(conn[nn]) != frontNodes->end()) {
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++) {
if (frontEdges->find(edge) != frontEdges->end()) {
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;
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);
det *= 6;
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (
int dd = 0;
dd != 3;
dd++) {
shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
}
ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
}
}
mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
};
}
}
}
}
private:
private:
};
boost::shared_ptr<Range> frontNodes;
boost::shared_ptr<Range> frontEdges;
static inline std::map<long int, MatrixDouble>
mapRefCoords;
};
struct SetIntegrationAtFrontFace {
SetIntegrationAtFrontFace(boost::shared_ptr<Range> front_nodes,
boost::shared_ptr<Range> front_edges)
: frontNodes(front_nodes), frontEdges(front_edges) {};
int order_col, int order_data) {
constexpr
bool debug =
false;
constexpr int numNodes = 3;
constexpr int numEdges = 3;
constexpr int refinementLevels = 4;
auto &m_field = fe_raw_ptr->
mField;
auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
auto fe_handle = fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
int rule = 2 * (order_data + 1);
}
}
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(
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) {
if (frontNodes->find(conn[nn]) != frontNodes->end()) {
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++) {
if (frontEdges->find(edge) != frontEdges->end()) {
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;
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);
auto det = t_normal.
l2();
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (
int dd = 0;
dd != 2;
dd++) {
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:
private:
};
boost::shared_ptr<Range> frontNodes;
boost::shared_ptr<Range> frontEdges;
static inline std::map<long int, MatrixDouble>
mapRefCoords;
};
return 0;
}
if (evalRhs)
if (evalLhs)
}
}
const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
const char *list_symm[] = {"symm", "not_symm"};
const char *list_stretches[] = {"linear", "log"};
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
CHKERR PetscOptionsScalar(
"-viscosity_alpha_omega",
"rot viscosity",
"",
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
PETSC_NULL);
CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
list_rots[choice_grad], &choice_grad, PETSC_NULL);
CHKERR PetscOptionsEList(
"-symm",
"symmetric variant",
"", list_symm, 2,
list_symm[choice_symm], &choice_symm, PETSC_NULL);
list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
CHKERR PetscOptionsBool(
"-set_singularity",
"set singularity",
"",
CHKERR PetscOptionsBool(
"-l2_user_base_scale",
"streach scale",
"",
CHKERR PetscOptionsBool(
"-dynamic_relaxation",
"dynamic time relaxation",
"",
CHKERR PetscOptionsInt(
"-contact_max_post_proc_ref_level",
"refinement level",
PETSC_NULL);
CHKERR PetscOptionsScalar(
"-griffith_energy",
"Griffith energy",
"",
ierr = PetscOptionsEnd();
}
break;
} else {
}
break;
default:
break;
};
MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
else
MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
MOFEM_LOG(
"EP", Sev::inform) <<
"Stretch " << list_stretches[choice_stretch];
? "yes"
: "no";
: "no";
? "yes"
: "no";
}
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);
}
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));
boost::make_shared<Range>(get_stress_trace_faces(
Range()));
#ifndef NDEBUG
#endif
auto add_broken_hdiv_field = [
this, meshset](
const std::string
field_name,
auto get_side_map_hdiv = [&]() {
return std::vector<
std::pair<EntityType,
>>{
{MBTET,
dofs_side_map);
}}
};
};
get_side_map_hdiv(), MB_TAG_DENSE,
MF_ZERO);
};
auto add_l2_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_h1_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_l2_field_by_range = [
this](
const std::string
field_name,
const int order,
const int dim,
const int field_dim,
Range &&
r) {
};
auto add_bubble_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
auto field_order_table =
const_cast<Field *
>(field_ptr)->getFieldOrderTable();
auto get_cgg_bubble_order_zero = [](int p) { return 0; };
auto get_cgg_bubble_order_tet = [](int p) {
};
field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
field_order_table[MBTRI] = get_cgg_bubble_order_zero;
field_order_table[MBTET] = get_cgg_bubble_order_tet;
};
auto add_user_l2_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
auto field_order_table =
const_cast<Field *
>(field_ptr)->getFieldOrderTable();
auto zero_dofs = [](int p) { return 0; };
field_order_table[MBVERTEX] = zero_dofs;
field_order_table[MBEDGE] = zero_dofs;
field_order_table[MBTRI] = zero_dofs;
field_order_table[MBTET] = dof_l2_tet;
};
}
auto project_ho_geometry = [&](auto field) {
};
auto get_adj_front_edges = [&](auto &front_edges) {
CHKERR moab.get_connectivity(front_edges, front_crack_nodes,
true);
front_crack_nodes);
crack_front_edges, moab::Interface::UNION);
crack_front_edges =
intersect(subtract(crack_front_edges, front_edges), meshset_ents);
crack_front_edges);
Range crack_front_edges_nodes;
CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
true);
crack_front_edges_nodes);
crack_front_edges_nodes =
subtract(crack_front_edges_nodes, front_crack_nodes);
Range crack_front_edges_with_both_nodes_not_at_front;
crack_front_edges_with_both_nodes_not_at_front,
moab::Interface::UNION);
crack_front_edges_with_both_nodes_not_at_front =
intersect(crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
crack_front_edges_with_both_nodes_not_at_front);
crack_front_edges_with_both_nodes_not_at_front = intersect(
crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
boost::make_shared<Range>(
crack_front_edges_with_both_nodes_not_at_front));
};
auto [front_vertices, front_adj_edges] = get_adj_front_edges(*
frontEdges);
(boost::format("crack_faces_%d.vtk") % rank).str(),
}
auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
double beta = 0;
PETSC_NULL);
MOFEM_LOG(
"EP", Sev::inform) <<
"Singularity eps " << beta;
for (auto edge : front_adj_edges) {
int num_nodes;
CHKERR moab.get_connectivity(edge, conn, num_nodes,
false);
double coords[6];
CHKERR moab.get_coords(conn, num_nodes, coords);
const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
coords[5] - coords[2]};
double dof[3] = {0, 0, 0};
if (front_vertices.find(conn[0]) != front_vertices.end()) {
for (
int dd = 0;
dd != 3;
dd++) {
}
} else if (front_vertices.find(conn[1]) != front_vertices.end()) {
for (
int dd = 0;
dd != 3;
dd++) {
}
}
const int idx = dit->get()->getEntDofIdx();
if (idx > 2) {
dit->get()->getFieldData() = 0;
} else {
dit->get()->getFieldData() = dof[idx];
}
}
}
};
}
auto add_field_to_fe = [this](const std::string fe,
};
}
Range front_elements_layer;
front_elements_layer,
moab::Interface::UNION);
front_elements.merge(front_elements_layer);
front_edges.clear();
moab::Interface::UNION);
}
}
}
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 = intersect(natural_bc_elements, meshset_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);
}
}
}
}
Range front_elements_layer;
front_elements_layer,
moab::Interface::UNION);
front_elements.merge(front_elements_layer);
front_edges.clear();
moab::Interface::UNION);
}
Range front_elements_faces;
true, front_elements_faces,
moab::Interface::UNION);
Range material_skeleton_faces =
subtract(front_elements_faces, front_elements_skin);
material_skeleton_faces.merge(intersect(front_elements_skin, body_skin));
#ifndef NDEBUG
front_elements_skin);
material_skeleton_faces);
#endif
}
}
auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
for (
int d : {0, 1, 2}) {
std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
->getSideDofsOnBrokenSpaceEntities(
dofs_to_remove);
dofs_to_remove);
}
};
CHKERR remove_dofs_on_broken_skin(
"ESHELBY_PLASTICITY");
}
auto set_zero_block = [&]() {
}
}
};
auto set_section = [&]() {
PetscSection section;
§ion);
CHKERR PetscSectionDestroy(§ion);
};
->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
}
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(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
}
theta = attr[3];
}
TractionBc::TractionBc(std::string name, std::vector<double> &attr,
: 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();
}
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);
}
std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
true);
(*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
}
}
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data + 1;
}
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * (p_data + 1);
}
};
CGGUserPolynomialBase::CGGUserPolynomialBase(
boost::shared_ptr<CachePhi> cache_phi_otr)
boost::typeindex::type_index type_index,
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
const auto base = this->cTx->bAse;
switch (this->cTx->sPace) {
CHKERR getValueHdivForCGGBubble(pts);
break;
data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
false);
&pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
CHKERR getValueL2AinsworthBase(pts);
break;
default:
}
}
CGGUserPolynomialBase::getValueHdivForCGGBubble(
MatrixDouble &pts) {
"Wrong base, should be USER_BASE");
}
const int nb_gauss_pts = pts.size2();
auto check_cache = [
this](
int order,
int nb_gauss_pts) ->
bool {
if (cachePhiPtr) {
return cachePhiPtr->get<0>() ==
order &&
cachePhiPtr->get<1>() == nb_gauss_pts;
} else {
return false;
}
};
if (check_cache(
order, nb_gauss_pts)) {
auto &mat = cachePhiPtr->get<2>();
phi.resize(mat.size1(), mat.size2(),
false);
} else {
shapeFun.resize(nb_gauss_pts, 4, false);
&pts(2, 0), nb_gauss_pts);
double diff_shape_fun[12];
phi.resize(nb_gauss_pts, 9 * nb_base_functions,
false);
nb_gauss_pts);
if (cachePhiPtr) {
cachePhiPtr->get<0>() =
order;
cachePhiPtr->get<1>() = nb_gauss_pts;
cachePhiPtr->get<2>().resize(
phi.size1(),
phi.size2(),
false);
noalias(cachePhiPtr->get<2>()) =
phi;
}
}
}
const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
auto bubble_cache =
boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0,
MatrixDouble());
fe->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe->getRuleHook = [](
int,
int,
int) {
return -1; };
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
fe->getOpPtrVector().push_back(
} else {
fe->getOpPtrVector().push_back(
}
if (calc_rates) {
} else {
fe->getOpPtrVector().push_back(
}
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
}
}
var_vec));
var_vec, MBMAXTYPE));
fe->getOpPtrVector().push_back(
} else {
fe->getOpPtrVector().push_back(
}
}
} else {
}
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto rule_contact = [](
int,
int,
int o) {
return -1; };
auto set_rule_contact = [refine](
int order_col, int order_data
) {
auto rule = 2 * order_data;
fe_raw_ptr->
gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
};
auto time_scale = boost::make_shared<DynamicRelaxationTimeScale>();
fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe_rhs);
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
} else {
fe_rhs->getOpPtrVector().push_back(
}
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
auto set_hybridisation = [&](auto &pip) {
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](
int,
int,
int) {
return -1;
};
op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
flux_mat_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
op_loop_skeleton_side->getOpPtrVector().push_back(
hybrid_ptr));
op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact = [&](auto &pip) {
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
};
auto add_ops_contact_rhs = [&](auto &pip) {
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(new OpTreeSearch(
nullptr));
contact_sfd_map_range_ptr));
pip.push_back(
};
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
pip.push_back(op_loop_skeleton_side);
};
CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
CHKERR set_contact(fe_rhs->getOpPtrVector());
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<
GAUSS>;
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
"BODY_FORCE", Sev::inform);
}
fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe_lhs);
if (add_elastic) {
fe_lhs->getOpPtrVector().push_back(
fe_lhs->getOpPtrVector().push_back(
} else {
fe_lhs->getOpPtrVector().push_back(
}
}
auto set_hybridisation = [&](auto &pip) {
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](
int,
int,
int) {
return -1;
};
op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(
boost::make_shared<double>(1.0), true, false));
pip.push_back(op_loop_skeleton_side);
};
auto set_contact = [&](auto &pip) {
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
};
auto add_ops_contact_lhs = [&](auto &pip) {
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(new OpTreeSearch(
nullptr));
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range(
"CONTACT_SDF",
SPACE_DIM - 1));
contact_sfd_map_range_ptr));
pip.push_back(
pip.push_back(
};
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
pip.push_back(op_loop_skeleton_side);
};
CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
CHKERR set_contact(fe_lhs->getOpPtrVector());
}
if (add_material) {
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
fe_rhs->getRuleHook = [](
int,
int,
int) {
return -1; };
fe_lhs->getRuleHook = [](
int,
int,
int) {
return -1; };
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) {
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
boost::shared_ptr<double> piola_scale_ptr(new double);
op_loop_domain_side->getOpPtrVector().push_back(
pip.push_back(op_loop_domain_side);
return std::make_pair(broken_data_ptr, piola_scale_ptr);
};
auto [broken_data_ptr, piola_scale_ptr] =
get_broken_op_side(fe_rhs->getOpPtrVector());
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(
{
boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt")
}));
fe_rhs->getOpPtrVector().push_back(
{
boost::make_shared<DynamicRelaxationTimeScale>(
"rotation_history.txt")
}));
{boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt")}
));
}
}
boost::shared_ptr<ContactTree> &fe_contact_tree
) {
auto get_body_range = [this](auto name, int dim, auto sev) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
MOFEM_LOG(
"EPSYNC", sev) <<
"Meshset: " << m_ptr->getMeshsetId() <<
" "
<< ents.size() << " entities";
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto calcs_side_traction = [&](auto &pip) {
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr(),
boost::make_shared<double>(1.0)));
pip.push_back(op_loop_domain_side);
};
auto add_contact_three = [&]() {
auto tree_moab_ptr = boost::make_shared<moab::Core>();
fe_contact_tree = boost::make_shared<ContactTree>(
get_body_range(
"CONTACT",
SPACE_DIM - 1, Sev::inform));
fe_contact_tree->getOpPtrVector().push_back(
contactDisp, contact_common_data_ptr->contactDispPtr()));
CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
fe_contact_tree->getOpPtrVector().push_back(
fe_contact_tree->getOpPtrVector().push_back(
new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
};
}
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_contact_bc = [&]() {
return op_loop_side;
};
}
boost::shared_ptr<FEMethod> null;
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
null);
null);
null);
} else {
null);
null);
null);
}
}
struct solve_elastic_setup {
bool set_ts_monitor) {
#ifdef PYTHON_SDF
auto setup_sdf = [&]() {
boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
char sdf_file_name[255] = "sdf.py";
sdf_file_name, 255, PETSC_NULL);
if (file_exists(sdf_file_name)) {
MOFEM_LOG(
"EP", Sev::inform) << sdf_file_name <<
" file found";
sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
} else {
MOFEM_LOG(
"EP", Sev::warning) << sdf_file_name <<
" file NOT found";
}
return sdf_python_ptr;
};
auto sdf_python_ptr = setup_sdf();
#endif
auto setup_ts_monitor = [&]() {
boost::shared_ptr<TsCtx>
ts_ctx;
if (set_ts_monitor) {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
}
return std::make_tuple(
ts_ctx);
};
auto setup_snes_monitor = [&]() {
SNES snes;
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
void *))SNESMonitorFields,
};
auto setup_section = [&]() {
PetscSection section_raw;
int num_fields;
CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
}
return section_raw;
};
auto set_vector_on_mesh = [&]() {
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
};
auto setup_schur_block_solver = [&]() {
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
schur_ptr =
}
return schur_ptr;
};
#ifdef PYTHON_SDF
return std::make_tuple(setup_sdf(), setup_ts_monitor(),
setup_snes_monitor(), setup_section(),
set_vector_on_mesh(), setup_schur_block_solver());
#else
return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
setup_section(), set_vector_on_mesh(),
setup_schur_block_solver());
#endif
}
};
PetscBool debug_model = PETSC_FALSE;
PETSC_NULL);
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";
};
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 TSSolve(ts, PETSC_NULL);
TetPolynomialBase::switchCacheBaseOff<HDIV>(
SNES snes;
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
PETSC_NULL);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
if (lin_solver_iterations > expected_lin_solver_iterations)
SETERRQ2(
"Expected number of iterations is different than expected %d > %d",
lin_solver_iterations, expected_lin_solver_iterations);
}
}
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;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Dynamic Relaxation Options",
"none");
CHKERR PetscOptionsScalar(
"-dynamic_final_time",
"dynamic relaxation final time", "", final_time,
&final_time, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-dynamic_delta_time",
"dynamic relaxation final time", "", delta_time,
&delta_time, PETSC_NULL);
CHKERR PetscOptionsInt(
"-dynamic_max_it",
"dynamic relaxation iterations",
"",
max_it, &max_it, PETSC_NULL);
CHKERR PetscOptionsBool(
"-dynamic_h1_update",
"update each ts step",
"",
ts_h1_update, &ts_h1_update, PETSC_NULL);
ierr = PetscOptionsEnd();
auto setup_ts_monitor = [&]() {
auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
return monitor_ptr;
};
auto monitor_ptr = setup_ts_monitor();
TetPolynomialBase::switchCacheBaseOn<HDIV>(
double ts_delta_time;
CHKERR TSGetTimeStep(ts, &ts_delta_time);
if (ts_h1_update) {
}
monitor_ptr->ts_u = PETSC_NULL;
CHKERR TSSetStepNumber(ts, 0);
CHKERR TSSetTimeStep(ts, ts_delta_time);
if (!ts_h1_update) {
}
CHKERR TSSolve(ts, PETSC_NULL);
if (!ts_h1_update) {
}
monitor_ptr->ts_u = PETSC_NULL;
break;
}
TetPolynomialBase::switchCacheBaseOff<HDIV>(
}
auto set_block = [&](auto name, int dim) {
std::map<int, Range> map;
auto set_tag_impl = [&](auto name) {
std::regex((boost::format("%s(.*)") % name).str())
);
for (auto bc : bcs) {
true);
map[bc->getMeshsetId()] =
r;
}
};
return std::make_pair(name, map);
};
auto set_skin = [&](auto &&map) {
for (
auto &
m : map.second) {
}
return map;
};
auto set_tag = [&](auto &&map) {
auto name = map.first;
int def_val[] = {-1};
MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
"create tag");
for (
auto &
m : map.second) {
"clear tag");
}
};
set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
}
Vec f_residual,
Vec var_vector,
std::vector<Tag> tags_to_transfer) {
if (
rval == MB_SUCCESS) {
}
int def_val[] = {0};
"CRACK", 1, MB_TYPE_INTEGER,
th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
int mark[] = {1};
tags_to_transfer.push_back(
th);
}
tags_to_transfer.push_back(
t);
}
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
auto post_proc_ptr =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
auto domain_ops = [&](auto &fe, int sense) {
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
auto piola_scale_ptr = boost::make_shared<double>(1.0);
fe.getOpPtrVector().push_back(
} else {
fe.getOpPtrVector().push_back(
}
if (var_vector) {
boost::make_shared<double>(1), vec));
boost::make_shared<double>(1), vec, MBMAXTYPE));
fe.getOpPtrVector().push_back(
} else {
fe.getOpPtrVector().push_back(
}
}
fe.getOpPtrVector().push_back(
if (auto op =
fe.getOpPtrVector().push_back(op);
}
struct OpSidePPMap :
public OpPPMap {
std::vector<EntityHandle> &map_gauss_pts,
DataMapVec data_map_scalar, DataMapMat data_map_vec,
DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
int sense)
:
OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
data_map_vec, data_map_mat, data_symm_map_mat),
tagSense(sense) {}
if (tagSense != OpPPMap::getSkeletonSense())
}
private:
int tagSense;
};
vec_fields[
"SpatialDisplacementL2"] =
dataAtPts->getSmallWL2AtPts();
vec_fields[
"SpatialDisplacementH1"] =
dataAtPts->getSmallWH1AtPts();
vec_fields[
"Omega"] =
dataAtPts->getRotAxisAtPts();
vec_fields[
"ContactDisplacement"] =
dataAtPts->getContactL2AtPts();
vec_fields[
"AngularMomentum"] =
dataAtPts->getLeviKirchhoffAtPts();
vec_fields[
"X"] =
dataAtPts->getLargeXH1AtPts();
if (var_vector) {
vec_fields[
"VarOmega"] =
dataAtPts->getVarRotAxisPts();
vec_fields["VarSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"VarSpatialDisplacementL2"], vec, MBTET));
}
if (f_residual) {
vec_fields["ResSpatialDisplacementL2"] =
boost::make_shared<MatrixDouble>();
spatialL2Disp, vec_fields[
"ResSpatialDisplacementL2"], vec, MBTET));
vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
rotAxis, vec_fields[
"ResOmega"], vec, MBTET));
}
mat_fields[
"PiolaStress"] =
dataAtPts->getApproxPAtPts();
if (var_vector) {
mat_fields[
"VarPiolaStress"] =
dataAtPts->getVarPiolaPts();
}
if (f_residual) {
mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
boost::make_shared<double>(1), vec));
boost::make_shared<double>(1), vec, MBMAXTYPE));
}
mat_fields_symm["LogSpatialStretch"] =
mat_fields_symm[
"SpatialStretch"] =
dataAtPts->getStretchTensorAtPts();
if (var_vector) {
mat_fields_symm["VarLogSpatialStretch"] =
}
if (f_residual) {
mat_fields_symm["ResLogSpatialStretch"] =
boost::make_shared<MatrixDouble>();
fe.getOpPtrVector().push_back(
MBTET));
}
}
fe.getOpPtrVector().push_back(
new OpSidePPMap(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
vec_fields,
mat_fields,
mat_fields_symm,
sense
)
);
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
};
post_proc_ptr->getOpPtrVector().push_back(
auto X_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
domain_ops(*(op_loop_side->getSideFEPtr()), sense);
post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
return post_proc_ptr;
};
auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
op_loop_domain_side->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr(),
boost::make_shared<double>(1.0)));
pip.push_back(op_loop_domain_side);
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
contactDisp, contact_common_data_ptr->contactDispPtr()));
pip.push_back(new OpTreeSearch(
&post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
};
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
post_proc_ptr->getOpPtrVector());
auto own_tets =
own_faces, moab::Interface::UNION);
auto get_post_negative = [&](auto &&ents) {
auto crack_faces_pos = ents;
auto crack_faces_neg = crack_faces_pos;
auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
for (
auto f : crack_on_proc_skin) {
tet = intersect(tet, own_tets);
int side_number, sense, offset;
offset);
if (sense == 1) {
crack_faces_neg.erase(
f);
} else {
crack_faces_pos.erase(
f);
}
}
return std::make_pair(crack_faces_pos, crack_faces_neg);
};
auto get_crack_faces = [&](auto crack_faces) {
auto get_adj = [&](auto e, auto dim) {
moab::Interface::UNION);
return adj;
};
auto tets = get_adj(crack_faces, 3);
auto faces = subtract(get_adj(tets, 2), crack_faces);
tets = subtract(tets, get_adj(faces, 3));
return subtract(crack_faces, get_adj(tets, 2));
};
auto [crack_faces_pos, crack_faces_neg] =
get_post_negative(intersect(own_faces, get_crack_faces(*
crackFaces)));
auto only_crack_faces = [&](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
return false;
}
return true;
};
auto only_negative_crack_faces = [&](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
return false;
}
return true;
};
auto get_adj_front = [&]() {
moab::Interface::UNION);
adj_front = intersect(adj_front, skeleton_faces);
adj_front = intersect(own_faces, adj_front);
return adj_front;
};
post_proc_ptr->setTagsToTransfer(tags_to_transfer);
post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
auto post_proc_begin =
post_proc_ptr->exeTestHook = only_crack_faces;
post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
post_proc_negative_sense_ptr, 0,
constexpr
bool debug =
false;
auto [adj_front_pos, adj_front_neg] =
auto only_front_faces_pos = [adj_front_pos](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (adj_front_pos.find(ent) == adj_front_pos.end()) {
return false;
}
return true;
};
auto only_front_faces_neg = [adj_front_neg](
FEMethod *fe_method_ptr) {
auto ent = fe_method_ptr->getFEEntityHandle();
if (adj_front_neg.find(ent) == adj_front_neg.end()) {
return false;
}
return true;
};
post_proc_ptr->exeTestHook = only_front_faces_pos;
post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
post_proc_negative_sense_ptr, 0,
}
CHKERR post_proc_end.writeFile(file.c_str());
}
const std::string file,
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, SPACE_DIM>::add(
auto hybrid_disp = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
auto hybrid_res = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"hybrid_disp", hybrid_disp}},
{},
{}
)
);
if (f_residual) {
auto hybrid_res = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
hybridSpatialDisp, hybrid_res,
post_proc_ptr->getOpPtrVector().push_back(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
{},
{{"res_hybrid", hybrid_res}},
{},
{}
)
);
}
auto post_proc_begin =
CHKERR post_proc_end.writeFile(file.c_str());
}
constexpr
bool debug =
true;
auto get_tags_vec = [&]() {
std::vector<Tag> tags(1);
auto create_and_clean = [&]() {
auto rval = moab.tag_get_handle(
"ReleaseEnergy", tags[0]);
if (
rval == MB_SUCCESS) {
moab.tag_delete(tags[0]);
}
double def_val[] = {0};
CHKERR moab.tag_get_handle(
"ReleaseEnergy", 1, MB_TYPE_DOUBLE, tags[0],
MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
};
return tags;
};
auto tags = get_tags_vec();
auto get_adj_front = [&]() {
moab::Interface::UNION);
adj_front = intersect(adj_front, skeleton_faces);
return adj_front;
};
auto get_nb_front_faces_on_proc = [&](auto &&adj_front) {
int send_size = adj_front.size();
MPI_Allgather(&send_size, 1, MPI_INT, recv_data.data(), 1, MPI_INT,
return std::pair(recv_data, adj_front);
};
auto get_snes = [&](TS ts) {
SNES snes;
return snes;
};
auto get_ksp = [&](SNES snes) {
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
return ksp;
};
auto integrate_energy = [&](auto x, auto release_energy_ptr) {
auto set_volume = [&]() {
auto fe_ptr =
boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe_ptr->ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
fe_ptr->getOpPtrVector().push_back(
return fe_ptr;
};
};
auto calc_release_energy = [&](
auto t,
auto &&tuple_of_vectors,
auto adj_front, double &energy, double &area) {
MOFEM_LOG(
"EP", Sev::inform) <<
"Calculate face release energy";
auto copy_is = [&](
auto is,
auto a,
auto b) {
const PetscInt *is_array;
PetscInt is_size;
CHKERR ISGetLocalSize(is, &is_size);
CHKERR ISGetIndices(is, &is_array);
double *pa;
const double *pb;
CHKERR VecGetArrayRead(b, &pb);
for (
int i = 0;
i != is_size; ++
i) {
pa[is_array[
i]] = pb[is_array[
i]];
}
CHKERR VecRestoreArrayRead(b, &pb);
CHKERR ISRestoreIndices(is, &is_array);
};
auto axpy_is = [&](
auto is,
double alpha,
auto a,
auto b) {
const PetscInt *is_array;
PetscInt is_size;
CHKERR ISGetLocalSize(is, &is_size);
CHKERR ISGetIndices(is, &is_array);
double *pa;
const double *pb;
CHKERR VecGetArrayRead(b, &pb);
for (
int i = 0;
i != is_size; ++
i) {
pa[is_array[
i]] += alpha * pb[is_array[
i]];
}
CHKERR VecRestoreArrayRead(b, &pb);
CHKERR ISRestoreIndices(is, &is_array);
};
auto zero_is = [&](
auto is,
auto a) {
const PetscInt *is_array;
PetscInt is_size;
CHKERR ISGetLocalSize(is, &is_size);
CHKERR ISGetIndices(is, &is_array);
double *pa;
for (
int i = 0;
i != is_size; ++
i) {
}
CHKERR ISRestoreIndices(is, &is_array);
};
auto estimate_energy = [&](auto x) {
auto release_energy_ptr = boost::make_shared<double>(0);
"integrate_energy");
double area = 0;
area += t_normal.
l2() / 2;
}
std::array<double, 2> send_data{area, *release_energy_ptr};
std::array<double, 2> recv_data{0, 0};
MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
return std::pair(recv_data[1], recv_data[0]);
};
auto set_tag = [&](auto &&release_energy, auto adj_front) {
MOFEM_LOG(
"EP", Sev::inform) <<
"Energy release " << release_energy;
CHKERR moab.tag_clear_data(tags[0], own, &release_energy);
};
auto [x,
f] = tuple_of_vectors;
Mat mA;
auto ksp = get_ksp(get_snes(ts));
CHKERR KSPGetOperators(ksp, &mA, PETSC_NULL);
std::vector<boost::weak_ptr<NumeredDofEntity>> piola_skelton_dofs_vec;
->getSideDofsOnBrokenSpaceEntities(
->isCreateProblemBrokenFieldAndRankLocal(piola_skelton_dofs_vec,
skeleton_piola_is_local);
->isCreateProblemBrokenFieldAndRank(piola_skelton_dofs_vec,
skeleton_piola_is_global);
skeleton_hybrid_is_local, &faces);
skeleton_hybrid_is_global, &faces);
CHKERR copy_is(skeleton_piola_is_local, x,
t);
CHKERR zero_is(skeleton_piola_is_local,
f);
CHKERR zero_is(skeleton_hybrid_is_local,
f);
std::vector<double> save_mat_storage(*mat_storage);
std::vector<double> save_mat_precon_storage(*mat_precon_storage);
CHKERR MatZeroRowsColumnsIS(mA, skeleton_piola_is_global, 1, PETSC_NULL,
PETSC_NULL);
CHKERR MatZeroRowsColumnsIS(mA, skeleton_hybrid_is_global, 1, PETSC_NULL,
PETSC_NULL);
PetscBool factor_schur = PETSC_FALSE;
&factor_schur, PETSC_NULL);
if (factor_schur == PETSC_TRUE) {
schur_skeleton_hybrid_is, &faces);
CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatZeroRowsColumnsIS(
S, schur_skeleton_hybrid_is, 1, PETSC_NULL,
PETSC_NULL);
}
CHKERR copy_is(skeleton_piola_is_local,
f,
t);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
std::copy(save_mat_storage.begin(), save_mat_storage.end(),
mat_storage->begin());
std::copy(save_mat_precon_storage.begin(), save_mat_precon_storage.end(),
mat_precon_storage->begin());
};
auto [x,
f] = tuple_of_vectors;
std::tie(energy, area) = estimate_energy(x);
CHKERR set_tag(energy, adj_front);
};
auto run_faces = [&](auto &&p) {
auto [recv_data, adj_front] = p;
std::map<int, std::tuple<double, double, Range, SmartPetscObj<Vec>>>
release_by_energy;
auto solve_and_add = [&](auto &&face, auto id) {
double energy, area;
std::make_tuple(face_x, face_f),
face, energy, area);
if (release_by_energy.find(id) == release_by_energy.end()) {
release_by_energy[id] = std::make_tuple(energy, area, face, vec);
} else {
auto [e,
a, faces, vec] = release_by_energy.at(
id);
if (e < energy) {
faces.merge(face);
release_by_energy[id] = std::make_tuple(energy, area, faces, vec);
}
}
};
ParallelComm *pcomm =
for (
auto f = 0;
f < recv_data[p]; ++
f) {
int id;
MOFEM_LOG(
"SYNC", Sev::noisy) <<
"edge id " << id;
solve_and_add(
Range(),
id);
}
}
auto face =
Range(adj_front[
f], adj_front[
f]);
moab::Interface::UNION),
"get adj");
int owner;
CHK_MOAB_THROW(pcomm->get_owner_handle(edges[0], owner, owner_handle),
"get handle");
MOFEM_LOG(
"SYNC", Sev::noisy) <<
"owner " << owner <<
"edge id " << id;
solve_and_add(face, id);
}
for (
auto f = 0;
f < recv_data[p]; ++
f) {
int id;
MOFEM_LOG(
"SYNC", Sev::noisy) <<
"edge id " << id;
solve_and_add(
Range(),
id);
}
}
return release_by_energy;
};
auto release_by_energy = run_faces(
);
std::map<double, std::tuple<double, EntityHandle, Range, SmartPetscObj<Vec>>>
sorted_release_by_energy;
for (
auto &
m : release_by_energy) {
auto [energy, area, faces, vec] =
m.second;
sorted_release_by_energy[energy] =
std::make_tuple(area,
m.first, faces, vec);
}
double total_area = 0;
double total_energy_sum = 0;
for (
auto m = sorted_release_by_energy.rbegin();
m != sorted_release_by_energy.rend(); ++
m) {
auto [area, id, faces, vec] =
m->second;
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
auto release_energy_ptr = boost::make_shared<double>(0);
CHKERR integrate_energy(x, release_energy_ptr);
std::array<double, 2> send_data{area, *release_energy_ptr};
std::array<double, 2> recv_data{0, 0};
MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
if (recv_data[1] < total_energy_sum) {
double release_energy = 0;
} else {
total_energy_sum = recv_data[1];
total_area += recv_data[0];
}
<< "Aggregated energy release " << total_area << " " << recv_data[1];
}
if (tet.size() != 1) {
MOFEM_LOG(
"EP", Sev::error) <<
"tet.size()!=1";
continue;
}
double release_energy;
release_energy /= total_energy_sum;
}
}
}
bool set_orientation) {
constexpr
bool debug =
true;
constexpr auto sev = Sev::inform;
Range geometry_edges_verts;
geometry_edges_verts, true);
true);
*
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
auto get_tags_vec = [&](auto tag_name, int dim) {
std::vector<Tag> tags(1);
if (dim > 3)
auto create_and_clean = [&]() {
auto rval = moab.tag_get_handle(tag_name, tags[0]);
if (
rval == MB_SUCCESS) {
moab.tag_delete(tags[0]);
}
double def_val[] = {0., 0., 0.};
CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
};
return tags;
};
auto get_adj_front = [&](bool subtract_crack) {
adj_front, moab::Interface::UNION);
if (subtract_crack)
return adj_front;
};
auto th_front_position = get_tags_vec("FrontPosition", 3);
auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
auto get_crack_adj_tets = [&](
auto r) {
Range crack_faces_conn_tets;
true, crack_faces_conn_tets,
moab::Interface::UNION);
return crack_faces_conn_tets;
};
auto get_layers_for_sides = [&](auto &side) {
std::vector<Range> layers;
auto get = [&]() {
auto get_adj = [&](
auto &
r,
int dim) {
moab::Interface::UNION);
return adj;
};
auto get_tets = [&](
auto r) {
return get_adj(
r,
SPACE_DIM); };
true);
Range front_faces = get_adj(front_nodes, 2);
auto front_tets = get_tets(front_nodes);
auto front_side = intersect(side, front_tets);
layers.push_back(front_side);
for (;;) {
adj_faces = intersect(adj_faces, front_faces);
auto adj_faces_tets = get_tets(adj_faces);
adj_faces_tets = intersect(adj_faces_tets, front_tets);
layers.push_back(unite(layers.back(), adj_faces_tets));
if (layers.back().size() == layers[layers.size() - 2].size()) {
break;
}
}
};
return layers;
};
auto layers_top = get_layers_for_sides(sides_pair.first);
auto layers_bottom = get_layers_for_sides(sides_pair.second);
#ifndef NDEBUG
"crack_tets_" +
sides_pair.second);
MOFEM_LOG(
"SELF", sev) <<
"Nb. layers " << layers_top.size();
for (
auto &
r : layers_top) {
MOFEM_LOG(
"SELF", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
"layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk",
r);
}
for (
auto &
r : layers_bottom) {
MOFEM_LOG(
"SELF", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
"layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk",
r);
}
}
#endif
auto get_cross = [&](
auto t_dir,
auto f) {
return t_cross;
};
auto get_sense = [&](
auto f,
auto e) {
int side, sense, offset;
"get sense");
return std::make_tuple(side, sense, offset);
};
auto calculate_edge_direction = [&](auto e) {
int num_nodes;
std::array<double, 6> coords;
&coords[0], &coords[1], &coords[2]};
&coords[3], &coords[4], &coords[5]};
t_dir(
i) = t_p1(
i) - t_p0(
i);
return t_dir;
};
auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
auto front_faces,
auto &sides_pair,
auto th_position) {
Tag th_face_energy;
auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
auto &edge_face_max_energy_map) {
for (auto e : front_edges) {
faces = intersect(faces, front_faces);
std::vector<double> face_energy(faces.size());
face_energy.data());
auto max_energy_it =
std::max_element(face_energy.begin(), face_energy.end());
double max_energy =
max_energy_it != face_energy.end() ? *max_energy_it : 0;
&max_energy);
edge_face_max_energy_map[e] =
std::make_tuple(faces[max_energy_it - face_energy.begin()],
max_energy, static_cast<double>(0));
<< "Edge " << e << " max energy " << max_energy;
};
};
auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
auto up_down_face = [&](
auto &face_angle_map_up,
auto &face_angle_map_down
) {
for (
auto &
m : edge_face_max_energy_map) {
auto [max_face, energy, opt_angle] =
m.second;
faces = intersect(faces, front_faces);
false, adj_tets,
moab::Interface::UNION);
if (adj_tets.size()) {
false, adj_tets,
moab::Interface::UNION);
if (adj_tets.size()) {
adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
moab::Interface::UNION);
adj_tets_faces = intersect(adj_tets_faces, faces);
if (adj_tets_faces.size() != 3 && adj_tets_faces.size() != 2) {
<< "Wrong number of crack faces " << adj_tets_faces.size()
<< " " << adj_tets.size();
"Wrong number of crack faces");
}
auto t_cross_max =
get_cross(calculate_edge_direction(e), max_face);
auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
t_cross_max(
i) *= sense_max;
for (
auto t : adj_tets) {
adj_tets_faces = intersect(adj_tets_faces, faces);
adj_tets_faces =
subtract(adj_tets_faces,
Range(max_face, max_face));
if (adj_tets_faces.size() == 1) {
auto t_cross = get_cross(calculate_edge_direction(e),
adj_tets_faces[0]);
auto [side, sense, offset] =
get_sense(adj_tets_faces[0], e);
double dot = t_cross(
i) * t_cross_max(
i);
auto angle = std::acos(dot);
double energy;
th_face_energy, adj_tets_faces, &energy);
auto [side_face, sense_face, offset_face] =
if (sense_face > 0) {
face_angle_map_up[e] =
std::make_tuple(energy, angle, adj_tets_faces[0]);
} else {
face_angle_map_down[e] =
std::make_tuple(energy, -angle, adj_tets_faces[0]);
}
}
}
}
}
}
};
auto calc_optimal_angle_impl = [&](
double a0,
double e0,
double a1,
double e1,
double a2,
double e2) {
b[0] = e1;
b[1] = e2;
b[2] = e0;
if (b[0] > -std::numeric_limits<float>::epsilon()) {
return std::make_pair(e0, 0.);
} else {
auto optimal_angle = -b[1] / (2 * b[0]);
auto optimal_energy = b[0] * optimal_angle * optimal_angle +
b[1] * optimal_angle + b[2];
#ifndef NDEBUG
<< "Optimal_energy " << optimal_energy << " ( " << e0 << " "
<< (optimal_energy - e0) / e0 << "% ) " << optimal_angle
<<
" e1 : a1 " << e1 <<
" : " <<
a1 <<
" e2 : a2 " << e2
double angle = -M_PI;
for (
double a = -M_PI;
a < M_PI;
a += M_PI / 20.) {
double energy = b[0] *
a *
a + b[1] *
a + b[2];
MOFEM_LOG(
"SELF", sev) <<
"PlotAngles " <<
a <<
" " << energy;
}
MOFEM_LOG(
"SELF", sev) <<
"PlotAngles " << endl;
MOFEM_LOG(
"SELF", sev) <<
"AnglePts " <<
a1 <<
" " << e1;
MOFEM_LOG(
"SELF", sev) <<
"AnglePts " <<
a0 <<
" " << e0;
MOFEM_LOG(
"SELF", sev) <<
"AnglePts " <<
a2 <<
" " << e2;
MOFEM_LOG(
"SELF", sev) <<
"AnglePts " << endl;
}
#endif // NDEBUG
return std::make_pair(optimal_energy, optimal_angle);
}
};
#ifndef NDEBUG
auto [opt_e0, opt_a0] = calc_optimal_angle_impl(1, 1, 0, 0, 3, -3);
if (std::abs(opt_e0 - 1) > std::numeric_limits<double>::epsilon()) {
}
if (std::abs(opt_a0 - 1) > std::numeric_limits<double>::epsilon()) {
}
#endif // NDEBUG
}
auto calc_optimal_angle = [&](
auto &face_angle_map_up,
auto &face_angle_map_down
) {
for (
auto &
m : edge_face_max_energy_map) {
auto &[max_face, e0,
a0] =
m.second;
if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
face_angle_map_down.find(e) == face_angle_map_down.end()) {
} else {
auto [e1,
a1, face_up] = face_angle_map_up.at(e);
auto [e2,
a2, face_down] = face_angle_map_down.at(e);
auto [optimal_energy, optimal_angle] =
calc_optimal_angle_impl(
a0, e0,
a1, e1,
a2, e2);
e0 = optimal_energy;
}
}
}
};
std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
face_angle_map_up;
std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
face_angle_map_down;
CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
auto th_angle = get_tags_vec("Angle", 1);
for (
auto &
m : face_angle_map_up) {
auto [e,
a, face] =
m.second;
up.insert(face);
}
for (
auto &
m : face_angle_map_down) {
auto [e,
a, face] =
m.second;
down.insert(face);
}
for (
auto &
m : edge_face_max_energy_map) {
auto [face, e, angle] =
m.second;
max_energy_faces.insert(face);
&angle);
}
max_energy_faces);
}
}
};
auto sort_by_energy = [&](auto &edge_face_max_energy_map) {
std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
sort_by_energy;
for (
auto m : edge_face_max_energy_map) {
auto [max_face, energy, opt_angle] =
m.second;
sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
}
return sort_by_energy;
};
auto set_face_orientation = [&](auto &sort_by_energy, auto &layers,
auto &set_vertices, double &all_quality) {
body_skin, 1, false, body_skin_edges, moab::Interface::UNION);
Range boundary_skin_verts;
boundary_skin_verts, true);
auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
auto t_edge_dir = calculate_edge_direction(e);
auto [side, sense, offset] = get_sense(
f, e);
t_edge_dir.normalize();
t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
return std::make_tuple(t_normal, t_rotated_normal);
};
for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
++it) {
auto energy = it->first;
auto [e, max_face, opt_angle] = it->second;
auto [t_normal, t_rotated_normal] =
get_rotated_normal(e, max_face, opt_angle);
rotated_normals.insert(max_face);
adj_tets);
#ifndef NDEBUG
if (adj_tets.size() != 2) {
<< "Wrong number of tets adjacent to max face "
<< adj_tets.size();
"Wrong number of crack faces");
}
#endif // NDEBUG
adj_front_faces);
adj_front_faces = subtract(adj_front_faces, *
crackFaces);
adj_tets, 2, false, adj_tet_faces, moab::Interface::UNION);
adj_tet_faces =
intersect(adj_tet_faces, adj_front_faces);
#ifndef NDEBUG
if (adj_tet_faces.size() != 3) {
MOFEM_LOG(
"SELF", sev) <<
"Wrong number of faces adj. to test "
<< adj_tet_faces.size();
"Wrong number of crack faces");
}
#endif // NDEBUG
Range adj_front_faces_edges;
adj_front_faces_edges,
moab::Interface::UNION);
std::array<EntityHandle, 3> processed_faces{0, max_face, 0};
int set_index = 0;
auto adj_tets_layer = intersect(adj_tets,
l);
if (adj_tets_layer.empty()) {
continue;
}
adj_tets_layer, 2, false, faces_layer, moab::Interface::UNION);
faces_layer = intersect(faces_layer, adj_front_faces);
faces_layer = subtract(faces_layer,
Range(max_face, max_face));
if (set_index > 0) {
faces_layer = subtract(
faces_layer,
Range(processed_faces[0], processed_faces[0]));
}
#ifndef NDEBUG
if (faces_layer.size() != 1) {
<< "Wrong number of faces to move " << faces_layer.size();
"Wrong number of crack faces");
}
#endif // NDEBUG
processed_faces[set_index] = faces_layer[0];
set_index = 2;
}
for (
auto f : {0, 1, 2}) {
if (processed_faces[
f] == 0)
break;
1, vertex, true);
if (
rval != MB_SUCCESS) {
<<
"get_connectivity failed " <<
f <<
" " "can noy get connectivity");
}
vertex = subtract(vertex, front_vertex);
if (set_vertices.find(vertex[0]) != set_vertices.end()) {
processed_faces[1] = 0;
processed_faces[2] = 0;
processed_faces[0] = 0;
processed_faces[2] = 0;
} else {
processed_faces[0] = 0;
processed_faces[1] = 0;
}
}
}
sort_by_quality;
int face_idx = 0;
for (
auto f : processed_faces) {
continue;
if (
rval != MB_SUCCESS) {
<<
"get_connectivity failed " <<
f <<
" " "can noy get connectivity");
}
vertex = subtract(vertex, front_vertex);
#ifndef NDEBUG
if (vertex.size() != 1) {
<< "Wrong number of vertex to move " << vertex.size();
}
#endif // NDEBUG
t_vertex_coords;
vertex_edges);
vertex_edges = subtract(vertex_edges, adj_front_faces_edges);
if (boundary_skin_verts.size() &&
boundary_skin_verts.find(vertex[0]) !=
boundary_skin_verts.end()) {
vertex_edges = intersect(vertex_edges, body_skin_edges);
}
if (geometry_edges_verts.size() &&
geometry_edges_verts.find(vertex[0]) !=
geometry_edges_verts.end()) {
MOFEM_LOG(
"SELF", sev) <<
"Geometry edge vertex";
vertex_edges = intersect(vertex_edges, geometry_edges);
}
if (crack_faces_verts.size() &&
crack_faces_verts.find(vertex_edges[0]) !=
crack_faces_verts.end()) {
MOFEM_LOG(
"SELF", sev) <<
"Crack face vertex";
vertex_edges = intersect(vertex_edges, crack_faces_edges);
}
std::map<EntityHandle, FTensor::Tensor1<double, SPACE_DIM>>
node_positions_map;
for (auto e : vertex_edges) {
edge_conn = subtract(edge_conn, vertex);
#ifndef NDEBUG
if (edge_conn.size() != 1) {
<< "Wrong number of edge conn " << edge_conn.size();
}
#endif // NDEBUG
auto it = set_vertices.find(vertex[0]);
if (it != set_vertices.end()) {
node_positions_map[e](
i) = std::get<2>(it->second)(
i);
} else {
t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
auto a = (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
auto b = (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
<< face_idx << " edge: " << e << " gamma: " << gamma;
if (std::isnormal(gamma) &&
gamma > -std::numeric_limits<double>::epsilon() &&
gamma < 1 - std::numeric_limits<double>::epsilon()) {
for (
auto i : {0, 1, 2}) {
node_positions_map[e](
i) =
gamma * (t_v3(
i) - t_vertex_coords(
i));
}
}
}
}
if (vertex_edges.size() == 0) {
node_positions_map[0](
i) = 0;
}
adj_vertex_tets);
for (auto &p : node_positions_map) {
double min_quality = 1;
for (
auto t : adj_vertex_tets) {
int num_nodes;
true);
std::array<double, 12> coords;
coords.data());
int idx = 0;
for (; idx < num_nodes; ++idx) {
if (conn[idx] == vertex[0]) {
break;
}
}
if (set_vertices.find(vertex[0]) == set_vertices.end()) {
for (auto ii : {0, 1, 2}) {
coords[3 * idx + ii] += p.second(ii);
}
}
for (auto vv : {0, 1, 2, 3}) {
auto it = set_vertices.find(conn[vv]);
if (it != set_vertices.end()) {
auto &[
f, energy, t_position] = it->second;
for (auto ii : {0, 1, 2})
coords[3 * vv + ii] += t_position(ii);
}
}
double q = Tools::volumeLengthQuality(coords.data());
if (!std::isnormal(q))
q = -2;
min_quality = std::min(min_quality, q);
}
sort_by_quality[min_quality] =
std::make_tuple(vertex[0],
f,
p.second(0), p.second(1), p.second(2)});
}
++face_idx;
}
for (auto s : sort_by_quality) {
auto &[vertex, face, t_position] = s.second;
<< "Quality: " << s.first << " vertex " << vertex << " face "
<< face << " point: " << t_position;
}
}
for (auto it = sort_by_quality.rbegin(); it != sort_by_quality.rend();
++it) {
auto &[vertex, face, t_position] = it->second;
all_quality = std::min(all_quality, it->first);
<< "Set quality: " << it->first << " vertex " << vertex
<< " face " << face << " point: " << t_position;
set_vertices[vertex] = std::make_tuple(face, energy, t_position);
break;
}
}
};
std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
edge_face_max_energy_map;
CHKERR find_maximal_face_energy(front_edges, front_faces,
edge_face_max_energy_map);
CHKERR calculate_face_orientation(edge_face_max_energy_map);
auto edge_map_sort_by_energy = sort_by_energy(edge_face_max_energy_map);
double top_quality = 1;
set_vertices_top;
CHKERR set_face_orientation(edge_map_sort_by_energy, layers_top,
set_vertices_top, top_quality);
double bottom_quality = 1;
set_vertices_bottom;
CHKERR set_face_orientation(edge_map_sort_by_energy, layers_bottom,
set_vertices_bottom, bottom_quality);
MOFEM_LOG(
"SELF", sev) <<
"Top quality " << top_quality;
MOFEM_LOG(
"SELF", sev) <<
"Bottom quality " << bottom_quality;
auto set_tag = [&](auto &set_vertices, auto th_position) {
for (
auto m : set_vertices) {
auto &[
f, energy, t_p] =
m.second;
&energy);
}
};
if (top_quality > bottom_quality) {
CHKERR set_tag(set_vertices_top, th_position);
} else {
CHKERR set_tag(set_vertices_bottom, th_position);
}
auto th_front_top = get_tags_vec("FrontPositionTop", 3);
auto th_front_bottom = get_tags_vec("FrontPositionBottom", 3);
for (
auto m : set_vertices_top) {
auto &[
f, energy, t_p] =
m.second;
top_moved_faces.insert(
f);
&t_p(0));
}
Range bottom_moved_faces;
for (
auto m : set_vertices_bottom) {
auto &[
f, energy, t_p] =
m.second;
bottom_moved_faces.insert(
f);
&t_p(0));
}
for (
auto m : set_vertices_bottom) {
auto &[
f, energy, t_p] =
m.second;
bottom_moved_faces.insert(
f);
}
top_moved_faces);
bottom_moved_faces);
auto front_faces = get_adj_front(false);
}
};
CHKERR evaluate_face_energy_and_set_orientation(
*
frontEdges, get_adj_front(
false), sides_pair, th_front_position);
}
SCATTER_FORWARD);
SCATTER_FORWARD);
SCATTER_FORWARD);
auto get_max_moved_faces = [&]() {
auto adj_front = get_adj_front(false);
std::vector<double> face_energy(adj_front.size());
face_energy.data());
for (
int i = 0;
i != adj_front.size(); ++
i) {
if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
max_moved_faces.insert(adj_front[
i]);
}
}
return boost::make_shared<Range>(max_moved_faces);
};
"projected_tets_" +
tets);
"max_moved_faces_" +
}
}
Tag th_front_position;
std::vector<double> coords(3 * verts.size());
std::vector<double> pos(3 * verts.size());
for (
int i = 0;
i != 3 * verts.size(); ++
i) {
}
double zero[] = {0., 0., 0.};
}
constexpr
bool debug =
true;
"set_coords_faces_" +
}
}
constexpr bool potential_crack_debug = false;
if constexpr (potential_crack_debug) {
true);
crack_front_verts);
true, crack_front_faces,
moab::Interface::UNION);
crack_front_faces = intersect(crack_front_faces, add_ents);
crack_front_faces);
}
auto get_crack_faces = [&]() {
} else {
}
};
auto get_extended_crack_faces = [&]() {
auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
ParallelComm *pcomm =
if (!pcomm->rank()) {
auto get_nodes = [&](auto &&e) {
"get connectivity");
return nodes;
};
auto get_adj = [&](auto &&e, auto dim,
auto t = moab::Interface::UNION) {
"get adj");
return adj;
};
body_ents);
auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
size_t s;
do {
s = crack_faces.size();
auto crack_face_nodes = get_nodes(crack_faces_org);
auto crack_faces_edges =
get_adj(crack_faces_org, 1, moab::Interface::UNION);
auto crack_skin =
auto crack_skin_nodes = get_nodes(crack_skin);
auto crack_skin_faces =
get_adj(crack_skin, 2, moab::Interface::UNION);
crack_skin_faces = subtract(crack_skin_faces, crack_faces_org);
crack_faces = crack_faces_org;
for (
auto f : crack_skin_faces) {
auto edges = intersect(
get_adj(
Range(
f,
f), 1, moab::Interface::UNION), crack_skin);
auto edge_conn = get_nodes(
Range(edges));
if (edges.size() == 2) {
auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
crack_faces_org);
if (faces.size() == 2) {
auto edge0_conn = get_nodes(
Range(edges[0], edges[0]));
auto edge1_conn = get_nodes(
Range(edges[1], edges[1]));
auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
crack_skin_nodes);
auto node_edges = subtract(intersect(
get_adj(
edges_conn, 1, moab::Interface::INTERSECT),
crack_faces_edges), crack_skin);
if (node_edges.size()) {
auto get_t_dir = [&](auto e_conn) {
&t_dir(0));
return t_dir;
};
get_t_dir(edge0_conn)(
i) + get_t_dir(edge1_conn)(
i);
t_crack_surface_ave_dir(
i) = 0;
for (auto e : node_edges) {
auto e_conn = get_nodes(
Range(e, e));
auto t_dir = get_t_dir(e_conn);
t_crack_surface_ave_dir(
i) += t_dir(
i);
}
auto dot = t_ave_dir(
i) * t_crack_surface_ave_dir(
i);
if (dot < -std::numeric_limits<double>::epsilon()) {
}
}
}
}
} else if (edges.size() == 3) {
}
}
crack_faces_org = crack_faces;
} while (s != crack_faces.size());
};
};
return get_faces_of_crack_front_verts(get_crack_faces());
};
#ifndef NDEBUG
constexpr
bool debug =
false;
"new_crack_surface_" +
get_extended_crack_faces());
}
#endif // NDEBUG
auto new_crack_faces = subtract(get_extended_crack_faces(), *
crackFaces);
};
auto crack_faces =
verts = subtract(verts, conn);
std::vector<double> coords(3 * verts.size());
double def_coords[] = {0., 0., 0.};
Tag th_org_coords;
"ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
}
};
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
auto post_proc_norm_rule_hook = [](
int,
int,
int p) ->
int {
return 2 * (p);
};
post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
post_proc_norm_fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
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(
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) {
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);
}
}
}
false, false);
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto force_bc = bc.second->forceBcPtr) {
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);
}
}
}
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 << " ) ";
};
}
}
