v0.14.0
EshelbianPlasticity.cpp

Eshelbian plasticity implementation

/**
* \file EshelbianPlasticity.cpp
* \example EshelbianPlasticity.cpp
*
* \brief Eshelbian plasticity implementation
*/
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
#include <CGGTonsorialBubbleBase.hpp>
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
#include <cholesky.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
#include <EshelbianContact.hpp>
extern "C" {
#include <phg-quadrule/quad.h>
}
namespace EshelbianPlasticity {
struct VolUserDataOperatorStabAssembly
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
: public FaceElementForcesAndSourcesCore::UserDataOperator {
using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
} // namespace EshelbianPlasticity
static auto get_range_from_block(MoFEM::Interface &m_field,
const std::string block_name, int dim) {
Range r;
auto mesh_mng = m_field.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
Range faces;
CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
faces, true),
"get meshset ents");
r.merge(faces);
}
return r;
};
static auto save_range(moab::Interface &moab, const std::string name,
const Range r) {
MoFEMFunctionBegin;
auto out_meshset = get_temp_meshset_ptr(moab);
CHKERR moab.add_entities(*out_meshset, r);
if (r.size()) {
CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
} else {
MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
}
MoFEMFunctionReturn(0);
};
template <>
typename MoFEM::OpBaseImpl<
SCHUR,
EshelbianPlasticity::VolUserDataOperatorStabAssembly>::MatSetValuesHook
MoFEM::OpBaseImpl<SCHUR,
EshelbianPlasticity::VolUserDataOperatorStabAssembly>::
matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
const EntitiesFieldData::EntData &row_data,
const EntitiesFieldData::EntData &col_data,
MatrixDouble &m) {
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
};
namespace EshelbianPlasticity {
struct SetIntegrationAtFrontVolume {
SetIntegrationAtFrontVolume() = default;
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
int order_col, int order_data) {
MoFEMFunctionBegin;
auto &m_field = fe_raw_ptr->mField;
auto vol_fe_ptr = dynamic_cast<VolFe *>(fe_raw_ptr);
auto fe_handle = vol_fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
MoFEMFunctionBegin;
int rule = 2 * order_data + 1;
if (rule < QUAD_3D_TABLE_SIZE) {
if (QUAD_3D_TABLE[rule]->dim != 3) {
SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
"wrong dimension");
}
if (QUAD_3D_TABLE[rule]->order < rule) {
SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
"wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
}
const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
auto &gauss_pts = vol_fe_ptr->gaussPts;
gauss_pts.resize(4, nb_gauss_pts, false);
cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
&gauss_pts(0, 0), 1);
cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
&gauss_pts(1, 0), 1);
cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
&gauss_pts(2, 0), 1);
cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
&gauss_pts(3, 0), 1);
auto &data = vol_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 {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"rule > quadrature order %d < %d", rule,
QUAD_3D_TABLE_SIZE);
}
MoFEMFunctionReturn(0);
};
CHKERR set_base_quadrature();
if (frontNodes->size() && EshelbianCore::setSingularity) {
auto get_singular_nodes = [&]() {
int num_nodes;
const EntityHandle *conn;
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++) {
EntityHandle edge;
CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
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) {
MoFEMFunctionBegin;
vol_fe_ptr->gaussPts.swap(ref_gauss_pts);
const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
// auto &dataH1 = vol_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();
// CHKERR ShapeMBTET(shape_ptr, &vol_fe_ptr->gaussPts(0, 0),
// &vol_fe_ptr->gaussPts(1, 0), &gaussPts(2, 0),
// nb_gauss_pts);
MoFEMFunctionReturn(0);
};
auto singular_nodes = get_singular_nodes();
if (singular_nodes.count()) {
auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
if (it_map_ref_coords != mapRefCoords.end()) {
CHKERR set_gauss_pts(it_map_ref_coords->second);
MoFEMFunctionReturnHot(0);
} else {
auto refine_quadrature = [&]() {
MoFEMFunctionBegin;
const int max_level = refinementLevels;
EntityHandle tet;
moab::Core moab_ref;
double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
EntityHandle nodes[4];
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);
MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
MoFEM::Interface &m_field_ref = mofem_ref_core;
CHKERR m_field_ref.getInterface<BitRefManager>()
->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
Range nodes_at_front;
for (int nn = 0; nn != numNodes; nn++) {
if (singular_nodes[nn]) {
EntityHandle ent;
CHKERR moab_ref.side_element(tet, 0, nn, ent);
nodes_at_front.insert(ent);
}
}
auto singular_edges = get_singular_edges();
EntityHandle meshset;
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
for (int ee = 0; ee != 6; ee++) {
if (singular_edges[ee]) {
EntityHandle ent;
CHKERR moab_ref.side_element(tet, 1, ee, ent);
CHKERR moab_ref.add_entities(meshset, &ent, 1);
}
}
// refine mesh
auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
for (int ll = 0; ll != max_level; ll++) {
Range edges;
CHKERR m_field_ref.getInterface<BitRefManager>()
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
BitRefLevel().set(), MBEDGE,
edges);
Range ref_edges;
CHKERR moab_ref.get_adjacencies(
nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
ref_edges = intersect(ref_edges, edges);
Range ents;
CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
ref_edges = intersect(ref_edges, ents);
Range tets;
CHKERR m_field_ref.getInterface<BitRefManager>()
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
ref_edges, BitRefLevel().set(ll + 1));
CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
CHKERR m_field_ref.getInterface<BitRefManager>()
->updateMeshsetByEntitiesChildren(meshset,
BitRefLevel().set(ll + 1),
meshset, MBEDGE, true);
}
// get ref coords
Range tets;
CHKERR m_field_ref.getInterface<BitRefManager>()
->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
BitRefLevel().set(), MBTET,
tets);
MatrixDouble ref_coords(tets.size(), 12, false);
int tt = 0;
for (Range::iterator tit = tets.begin(); tit != tets.end();
tit++, tt++) {
int num_nodes;
const EntityHandle *conn;
CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
}
auto &data = vol_fe_ptr->dataOnElement[H1];
const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
MatrixDouble &shape_n =
data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
int gg = 0;
for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
double *tet_coords = &ref_coords(tt, 0);
double det = Tools::tetVolume(tet_coords);
det *= 6;
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (int dd = 0; dd != 3; dd++) {
ref_gauss_pts(dd, gg) =
shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
}
ref_gauss_pts(3, gg) = vol_fe_ptr->gaussPts(3, ggg) * det;
}
}
mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
MoFEMFunctionReturn(0);
};
}
}
}
MoFEMFunctionReturn(0);
}
struct VolFe : public VolumeElementForcesAndSourcesCore {
using VolumeElementForcesAndSourcesCore::dataOnElement;
using VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore;
};
static inline constexpr int numNodes = 4;
static inline constexpr int numEdges = 6;
static inline constexpr int refinementLevels = 3;
static inline boost::shared_ptr<Range> frontNodes;
static inline boost::shared_ptr<Range> frontEdges;
static inline std::map<long int, MatrixDouble> mapRefCoords;
};
double EshelbianCore::exponentBase = exp(1);
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log_e;
boost::function<double(const double)> EshelbianCore::d_f =
EshelbianCore::d_f_log_e;
boost::function<double(const double)> EshelbianCore::dd_f =
EshelbianCore::dd_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_f =
EshelbianCore::inv_f_log;
boost::function<double(const double)> EshelbianCore::inv_d_f =
EshelbianCore::inv_d_f_log;
boost::function<double(const double)> EshelbianCore::inv_dd_f =
EshelbianCore::inv_dd_f_log;
MoFEMErrorCode
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
UnknownInterface **iface) const {
*iface = const_cast<EshelbianCore *>(this);
return 0;
}
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
if (evalRhs)
CHKERR evaluateRhs(data);
if (evalLhs)
CHKERR evaluateLhs(data);
MoFEMFunctionReturn(0);
}
EshelbianCore::EshelbianCore(MoFEM::Interface &m_field) : mField(m_field) {
CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
}
EshelbianCore::~EshelbianCore() = default;
MoFEMErrorCode EshelbianCore::getOptions() {
MoFEMFunctionBegin;
const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
PetscInt choice_rot = EshelbianCore::rotSelector;
PetscInt choice_grad = EshelbianCore::gradApproximator;
const char *list_stretches[] = {"linear", "log"};
PetscInt choice_stretch = StretchSelector::LOG;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
spaceOrder, &spaceOrder, PETSC_NULL);
CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
spaceH1Order, &spaceH1Order, PETSC_NULL);
CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
"", materialOrder, &materialOrder, PETSC_NULL);
CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
&alphaU, PETSC_NULL);
CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
&alphaW, PETSC_NULL);
CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
&alphaRho, PETSC_NULL);
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, NO_H1_CONFIGURATION + 1,
list_rots[choice_grad], &choice_grad, PETSC_NULL);
CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
&EshelbianCore::exponentBase, PETSC_NULL);
CHKERR PetscOptionsEList(
"-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
noStretch, &noStretch, PETSC_NULL);
CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
setSingularity, &setSingularity, PETSC_NULL);
// contact parameters
CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
"", contactRefinementLevels, &contactRefinementLevels,
PETSC_NULL);
ierr = PetscOptionsEnd();
CHKERRG(ierr);
EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
switch (EshelbianCore::stretchSelector) {
case StretchSelector::LINEAR:
EshelbianCore::f = f_linear;
EshelbianCore::d_f = d_f_linear;
EshelbianCore::dd_f = dd_f_linear;
EshelbianCore::inv_f = inv_f_linear;
EshelbianCore::inv_d_f = inv_d_f_linear;
EshelbianCore::inv_dd_f = inv_dd_f_linear;
break;
case StretchSelector::LOG:
if (EshelbianCore::exponentBase != exp(1)) {
EshelbianCore::f = f_log;
EshelbianCore::d_f = d_f_log;
EshelbianCore::dd_f = dd_f_log;
EshelbianCore::inv_f = inv_f_log;
EshelbianCore::inv_d_f = inv_d_f_log;
EshelbianCore::inv_dd_f = inv_dd_f_log;
} else {
EshelbianCore::f = f_log_e;
EshelbianCore::d_f = d_f_log_e;
EshelbianCore::dd_f = dd_f_log_e;
EshelbianCore::inv_f = inv_f_log;
EshelbianCore::inv_d_f = inv_d_f_log;
EshelbianCore::inv_dd_f = inv_dd_f_log;
}
break;
default:
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
break;
};
MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
MOFEM_LOG("EP", Sev::inform)
<< "Rotations " << list_rots[EshelbianCore::rotSelector];
MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
<< list_rots[EshelbianCore::gradApproximator];
if (exponentBase != exp(1))
MOFEM_LOG("EP", Sev::inform)
<< "Base exponent " << EshelbianCore::exponentBase;
else
MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
if (spaceH1Order == -1)
spaceH1Order = spaceOrder;
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
MoFEMFunctionBegin;
auto get_tets = [&]() {
Range tets;
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
return tets;
};
auto get_tets_skin = [&]() {
Range tets_skin_part;
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Range tets_skin;
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) {
// That mean, that hybrid field on all faces on which traction is applied,
// on other faces, or enforcing displacements as
// natural boundary condition.
if (bcSpatialTraction)
for (auto &v : *bcSpatialTraction) {
tets_skin = subtract(tets_skin, v.faces);
}
return tets_skin;
};
auto add_blockset = [&](auto block_name, auto &&tets_skin) {
auto crack_faces =
get_range_from_block(mField, "block_name", SPACE_DIM - 1);
tets_skin.merge(crack_faces);
return tets_skin;
};
auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
auto contact_range =
get_range_from_block(mField, block_name, SPACE_DIM - 1);
tets_skin = subtract(tets_skin, contact_range);
return tets_skin;
};
auto get_stress_trace_faces = [&](auto &&tets_skin) {
Range faces;
CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
faces, moab::Interface::UNION);
Range trace_faces = subtract(faces, tets_skin);
return trace_faces;
};
auto tets = get_tets();
// remove also contact faces, i.e. that is also kind of hybrid field but
// named but used to enforce contact conditions
auto trace_faces = get_stress_trace_faces(
subtract_blockset("CONTACT",
subtract_boundary_conditions(get_tets_skin()))
);
contactFaces = boost::make_shared<Range>(intersect(
trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
skeletonFaces =
boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
materialSkeletonFaces =
boost::make_shared<Range>(get_stress_trace_faces(Range()));
#ifndef NDEBUG
if (contactFaces->size())
CHKERR save_range(mField.get_moab(), "contact_faces.vtk", *contactFaces);
if (skeletonFaces->size())
CHKERR save_range(mField.get_moab(), "skeleton_faces.vtk", *skeletonFaces);
if (skeletonFaces->size())
CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
*materialSkeletonFaces);
#endif
auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
const int order) {
MoFEMFunctionBegin;
constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
auto get_side_map_hdiv = [&]() {
return std::vector<
std::pair<EntityType,
std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
>>{
{MBTET,
[&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
dofs_side_map);
}}
};
};
CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
MoFEMFunctionReturn(0);
};
auto add_l2_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MoFEMFunctionBegin;
CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
MoFEMFunctionReturn(0);
};
auto add_h1_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MoFEMFunctionBegin;
CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
MoFEMFunctionReturn(0);
};
auto add_l2_field_by_range = [this](const std::string field_name,
const int order, const int dim,
const int field_dim, Range &&r) {
MoFEMFunctionBegin;
CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
CHKERR mField.set_field_order(r, field_name, order);
MoFEMFunctionReturn(0);
};
auto add_bubble_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MoFEMFunctionBegin;
CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
MF_ZERO);
// Modify field
auto field_ptr = mField.get_field_structure(field_name);
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) {
return NBVOLUMETET_CCG_BUBBLE(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;
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
MoFEMFunctionReturn(0);
};
auto add_user_l2_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MoFEMFunctionBegin;
CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
MF_ZERO);
// Modify field
auto field_ptr = mField.get_field_structure(field_name);
auto field_order_table =
const_cast<Field *>(field_ptr)->getFieldOrderTable();
auto zero_dofs = [](int p) { return 0; };
auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
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;
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
MoFEMFunctionReturn(0);
};
// spatial fields
CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
CHKERR add_user_l2_field(rotAxis, spaceOrder - 1, 3);
CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
if (!skeletonFaces)
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
if (!contactFaces)
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
Range(*skeletonFaces));
CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
Range(*contactFaces));
// spatial displacement
CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
// material positions
CHKERR add_h1_field(materialH1Positions, 2, 3);
// Eshelby stress
CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
Range(*materialSkeletonFaces));
CHKERR mField.build_fields();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::projectGeometry(const EntityHandle meshset) {
MoFEMFunctionBegin;
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
return mField.loop_dofs(materialH1Positions, ent_method);
};
CHKERR project_ho_geometry();
auto get_skin = [&](auto &body_ents) {
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto &&skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto get_crack_front_edges = [&]() {
auto &moab = mField.get_moab();
auto crack_skin = get_skin(*crackFaces);
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
auto body_skin = filter_true_skin(get_skin(body_ents));
Range body_skin_edges;
CHKERR moab.get_adjacencies(body_skin, SPACE_DIM - 2, true, body_skin_edges,
moab::Interface::UNION);
crack_skin = subtract(crack_skin, body_skin_edges);
std::map<int, Range> received_ents;
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
crack_skin, &received_ents);
Range to_remove;
for (auto &m : received_ents) {
if (m.first != mField.get_comm_rank()) {
// if the same edges appear on two or more processors, they are not at
// the front
to_remove.merge(intersect(crack_skin, m.second));
}
}
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(to_remove);
crack_skin = subtract(crack_skin, to_remove);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
crack_skin);
return boost::make_shared<Range>(crack_skin);
};
auto get_adj_front_edges = [&](auto &front_edges) {
auto &moab = mField.get_moab();
Range front_crack_nodes;
CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
front_crack_nodes);
Range crack_front_edges;
CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
crack_front_edges, moab::Interface::UNION);
crack_front_edges = subtract(crack_front_edges, front_edges);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
crack_front_edges);
Range crack_front_edges_nodes;
CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
true);
crack_front_edges_nodes =
subtract(crack_front_edges_nodes, front_crack_nodes);
Range crack_front_edges_with_both_nodes_not_at_front;
CHKERR moab.get_adjacencies(crack_front_edges_nodes, SPACE_DIM - 2, false,
crack_front_edges_with_both_nodes_not_at_front,
moab::Interface::UNION);
crack_front_edges_with_both_nodes_not_at_front = intersect(
crack_front_edges_with_both_nodes_not_at_front, 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));
};
crackFaces = boost::make_shared<Range>(
get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
frontEdges = get_crack_front_edges();
auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
frontVertices = front_vertices;
frontAdjEdges = front_adj_edges;
auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
MoFEMFunctionBegin;
auto &moab = mField.get_moab();
double eps = 1;
double beta = 0;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "-singularity_eps", &beta,
PETSC_NULL);
MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
eps -= beta;
for (auto edge : front_adj_edges) {
int num_nodes;
const EntityHandle *conn;
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++) {
dof[dd] = -dir[dd] * eps;
}
} else if (front_vertices.find(conn[1]) != front_vertices.end()) {
for (int dd = 0; dd != 3; dd++) {
dof[dd] = +dir[dd] * eps;
}
}
for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(
mField, materialH1Positions, edge, dit)) {
const int idx = dit->get()->getEntDofIdx();
if (idx > 2) {
dit->get()->getFieldData() = 0;
} else {
dit->get()->getFieldData() = dof[idx];
}
}
}
MoFEMFunctionReturn(0);
};
if (setSingularity)
CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode
EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
MoFEMFunctionBegin;
// set finite element fields
auto add_field_to_fe = [this](const std::string fe,
const std::string field_name) {
MoFEMFunctionBegin;
CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
MoFEMFunctionReturn(0);
};
if (!mField.check_finite_element(elementVolumeName)) {
CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
elementVolumeName);
CHKERR add_field_to_fe(elementVolumeName, piolaStress);
CHKERR add_field_to_fe(elementVolumeName, bubbleField);
if (!noStretch)
CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
CHKERR add_field_to_fe(elementVolumeName, rotAxis);
CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
CHKERR add_field_to_fe(elementVolumeName, contactDisp);
CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
materialH1Positions);
// build finite elements data structures
CHKERR mField.build_finite_elements(elementVolumeName);
}
if (!mField.check_finite_element(materialVolumeElement)) {
Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
Range front_elements;
for (auto l = 0; l != frontLayers; ++l) {
Range front_elements_layer;
CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
front_elements_layer,
moab::Interface::UNION);
front_elements.merge(front_elements_layer);
front_edges.clear();
CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
SPACE_DIM - 2, true, front_edges,
moab::Interface::UNION);
}
CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
materialVolumeElement);
CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
CHKERR mField.build_finite_elements(materialVolumeElement);
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode
EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
MoFEMFunctionBegin;
auto set_fe_adjacency = [&](auto fe_name) {
MoFEMFunctionBegin;
parentAdjSkeletonFunctionDim2 =
boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
CHKERR mField.modify_finite_element_adjacency_table(
fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
MoFEMFunctionReturn(0);
};
// set finite element fields
auto add_field_to_fe = [this](const std::string fe,
const std::string field_name) {
MoFEMFunctionBegin;
CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
MoFEMFunctionReturn(0);
};
if (!mField.check_finite_element(naturalBcElement)) {
Range natural_bc_elements;
if (bcSpatialDispVecPtr) {
for (auto &v : *bcSpatialDispVecPtr) {
natural_bc_elements.merge(v.faces);
}
}
if (bcSpatialRotationVecPtr) {
for (auto &v : *bcSpatialRotationVecPtr) {
natural_bc_elements.merge(v.faces);
}
}
if (bcSpatialTraction) {
for (auto &v : *bcSpatialTraction) {
natural_bc_elements.merge(v.faces);
}
}
CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
naturalBcElement);
CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
CHKERR mField.build_finite_elements(naturalBcElement);
}
auto get_skin = [&](auto &body_ents) {
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto &&skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
if (!mField.check_finite_element(skinElement)) {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
auto skin = filter_true_skin(get_skin(body_ents));
CHKERR mField.add_finite_element(skinElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
CHKERR mField.modify_finite_element_add_field_data(skinElement,
spatialH1Disp);
CHKERR mField.modify_finite_element_add_field_data(skinElement,
materialH1Positions);
CHKERR mField.modify_finite_element_add_field_data(skinElement,
hybridSpatialDisp);
CHKERR mField.modify_finite_element_add_field_data(skinElement,
hybridMaterialDisp);
CHKERR mField.modify_finite_element_add_field_data(skinElement,
contactDisp);
CHKERR mField.build_finite_elements(skinElement);
}
if (!mField.check_finite_element(contactElement)) {
if (contactFaces) {
MOFEM_LOG("EP", Sev::inform)
<< "Contact elements " << contactFaces->size();
CHKERR mField.add_finite_element(contactElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(*contactFaces, MBTRI,
contactElement);
CHKERR add_field_to_fe(contactElement, piolaStress);
CHKERR add_field_to_fe(contactElement, spatialH1Disp);
CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
CHKERR add_field_to_fe(contactElement, contactDisp);
CHKERR set_fe_adjacency(contactElement);
CHKERR mField.build_finite_elements(contactElement);
}
}
if (!mField.check_finite_element(skeletonElement)) {
if (!skeletonFaces)
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
MOFEM_LOG("EP", Sev::inform)
<< "Skeleton elements " << skeletonFaces->size();
CHKERR mField.add_finite_element(skeletonElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(*skeletonFaces, MBTRI,
skeletonElement);
CHKERR add_field_to_fe(skeletonElement, piolaStress);
CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
CHKERR add_field_to_fe(contactElement, contactDisp);
CHKERR set_fe_adjacency(skeletonElement);
CHKERR mField.build_finite_elements(skeletonElement);
}
if (!mField.check_finite_element(materialSkeletonElement)) {
Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
Range front_elements;
for (auto l = 0; l != frontLayers; ++l) {
Range front_elements_layer;
CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
front_elements_layer,
moab::Interface::UNION);
front_elements.merge(front_elements_layer);
front_edges.clear();
CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
SPACE_DIM - 2, true, front_edges,
moab::Interface::UNION);
}
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
Range front_elements_faces;
CHKERR mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
true, front_elements_faces,
moab::Interface::UNION);
auto body_skin = filter_true_skin(get_skin(body_ents));
auto front_elements_skin = filter_true_skin(get_skin(front_elements));
Range material_skeleton_faces =
subtract(front_elements_faces, front_elements_skin);
material_skeleton_faces.merge(intersect(front_elements_skin, body_skin));
#ifndef NDEBUG
CHKERR save_range(mField.get_moab(), "front_elements.vtk", front_elements);
CHKERR save_range(mField.get_moab(), "front_elements_skin.vtk",
front_elements_skin);
CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
material_skeleton_faces);
#endif
CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(
material_skeleton_faces, MBTRI, materialSkeletonElement);
CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
CHKERR set_fe_adjacency(materialSkeletonElement);
CHKERR mField.build_finite_elements(materialSkeletonElement);
}
if (!mField.check_finite_element(crackElement)) {
auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
CHKERR mField.add_finite_element(crackElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(crack_faces, MBTRI,
crackElement);
CHKERR add_field_to_fe(crackElement, eshelbyStress);
CHKERR add_field_to_fe(crackElement, hybridMaterialDisp);
CHKERR mField.modify_finite_element_add_field_data(crackElement,
spatialH1Disp);
CHKERR mField.modify_finite_element_add_field_data(crackElement,
materialH1Positions);
CHKERR mField.modify_finite_element_add_field_data(crackElement,
hybridSpatialDisp);
CHKERR mField.modify_finite_element_add_field_data(crackElement,
hybridMaterialDisp);
CHKERR mField.modify_finite_element_add_field_data(crackElement,
contactDisp);
CHKERR set_fe_adjacency(crackElement);
CHKERR mField.build_finite_elements(crackElement);
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit,
const EntityHandle meshset) {
MoFEMFunctionBegin;
// find adjacencies between finite elements and dofs
CHKERR mField.build_adjacencies(bit, QUIET);
// Create coupled problem
dM = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
BitRefLevel().set());
CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dM, elementVolumeName);
CHKERR DMMoFEMAddElement(dM, naturalBcElement);
CHKERR DMMoFEMAddElement(dM, skeletonElement);
CHKERR DMMoFEMAddElement(dM, contactElement);
CHKERR DMMoFEMAddElement(dM, crackElement);
CHKERR DMMoFEMAddElement(dM, skinElement);
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
CHKERR DMSetUp(dM);
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
MoFEMFunctionBegin;
for (int d : {0, 1, 2}) {
std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
CHKERR mField.getInterface<ProblemsManager>()
->getSideDofsOnBrokenSpaceEntities(
dofs_to_remove, prb_name, ROW, piolaStress,
(*bcSpatialFreeTraction)[d], SPACE_DIM, d, d);
// remove piola dofs, i.e. traction free boundary
CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
dofs_to_remove);
CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
dofs_to_remove);
}
MoFEMFunctionReturn(0);
};
CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
// Create elastic sub-problem
dmElastic = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
if (!noStretch) {
CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
}
CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
CHKERR DMMoFEMAddElement(dmElastic, contactElement);
CHKERR DMMoFEMAddElement(dmElastic, skinElement);
CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
CHKERR DMSetUp(dmElastic);
// dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
// CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
// CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
// CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
// CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
// CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
// CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
// CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
// CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
// CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
// CHKERR DMSetUp(dmMaterial);
auto set_zero_block = [&]() {
MoFEMFunctionBegin;
if (!noStretch) {
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
}
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
if (!noStretch) {
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", bubbleField, bubbleField);
CHKERR
mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", piolaStress, piolaStress);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", bubbleField, piolaStress);
CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
"ELASTIC_PROBLEM", piolaStress, bubbleField);
}
MoFEMFunctionReturn(0);
};
auto set_section = [&]() {
MoFEMFunctionBegin;
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
&section);
CHKERR DMSetSection(dmElastic, section);
CHKERR DMSetGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(&section);
MoFEMFunctionReturn(0);
};
CHKERR set_zero_block();
CHKERR set_section();
dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
CHKERR DMSetUp(dmPrjSpatial);
CHKERR mField.getInterface<BcManager>()
->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
"PROJECT_SPATIAL", spatialH1Disp, true, false);
MoFEMFunctionReturn(0);
}
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;
MOFEM_LOG("EP", Sev::inform)
<< "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
MOFEM_LOG("EP", Sev::inform)
<< "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,
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;
MOFEM_LOG("EP", Sev::inform)
<< "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
MOFEM_LOG("EP", Sev::inform)
<< "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
}
MoFEMErrorCode
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
MoFEMFunctionBegin;
// get skin from all tets
Range tets;
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
Range tets_skin_part;
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, tets, false, tets_skin_part);
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Range tets_skin;
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;
// Do not remove dofs on which traction is applied
if (bcSpatialDispVecPtr)
for (auto &v : *bcSpatialDispVecPtr) {
if (v.flags[0])
(*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
if (v.flags[1])
(*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
if (v.flags[2])
(*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
}
// Do not remove dofs on which rotation is applied
if (bcSpatialRotationVecPtr)
for (auto &v : *bcSpatialRotationVecPtr) {
(*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);
}
if (bcSpatialTraction)
for (auto &v : *bcSpatialTraction) {
(*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);
}
// remove contact
for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
Range faces;
CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
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);
}
MoFEMFunctionReturn(0);
}
/**
* @brief Set integration rule on element
* \param order on row
* \param order on column
* \param order on data
*
* Use maximal oder on data in order to determine integration rule
*
*/
struct VolRule {
int operator()(int p_row, int p_col, int p_data) const {
// if (EshelbianCore::gradApproximator > MODERATE_ROT)
// return p_data + p_data + (p_data - 1);
// else
return 2 * p_data + 1;
}
};
struct FaceRule {
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)
: TetPolynomialBase(), cachePhiPtr(cache_phi_otr) {}
MoFEMErrorCode CGGUserPolynomialBase::query_interface(
boost::typeindex::type_index type_index,
BaseFunctionUnknownInterface **iface) const {
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
MoFEMErrorCode
CGGUserPolynomialBase::getValue(MatrixDouble &pts,
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
MoFEMFunctionBeginHot;
this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
MoFEMFunctionReturnHot(0);
}
if (pts.size1() < 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
const auto base = this->cTx->bAse;
EntitiesFieldData &data = this->cTx->dAta;
switch (this->cTx->sPace) {
case HDIV:
CHKERR getValueHdivForCGGBubble(pts);
break;
case L2:
data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
CHKERR Tools::shapeFunMBTET(
&*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
&pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
this->cTx->basePolynomialsType0 = Legendre_polynomials;
CHKERR getValueL2AinsworthBase(pts);
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
}
MoFEMFunctionReturnHot(0);
}
MoFEMErrorCode
CGGUserPolynomialBase::getValueHdivForCGGBubble(MatrixDouble &pts) {
MoFEMFunctionBegin;
// This should be used only in case USER_BASE is selected
if (this->cTx->bAse != USER_BASE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong base, should be USER_BASE");
}
// get access to data structures on element
EntitiesFieldData &data = this->cTx->dAta;
// Get approximation order on element. Note that bubble functions are only
// on tetrahedron.
const int order = data.dataOnEntities[MBTET][0].getOrder();
/// number of integration points
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>();
auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
phi.resize(mat.size1(), mat.size2(), false);
noalias(phi) = mat;
} else {
// calculate shape functions, i.e. barycentric coordinates
shapeFun.resize(nb_gauss_pts, 4, false);
CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
&pts(2, 0), nb_gauss_pts);
// derivatives of shape functions
double diff_shape_fun[12];
CHKERR ShapeDiffMBTET(diff_shape_fun);
const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
// get base functions and set size
MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
// finally calculate base functions
FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
&phi(0, 0), &phi(0, 1), &phi(0, 2),
&phi(0, 3), &phi(0, 4), &phi(0, 5),
&phi(0, 6), &phi(0, 7), &phi(0, 8));
CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
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;
}
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
const int tag, const bool do_rhs, const bool do_lhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
MoFEMFunctionBegin;
fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
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->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
// set integration rule
fe->getRuleHook = VolRule();
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
dataAtPts->physicsPtr = physicalEquations;
}
// calculate fields values
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
piolaStress, dataAtPts->getApproxPAtPts()));
fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
piolaStress, dataAtPts->getDivPAtPts()));
if (noStretch) {
fe->getOpPtrVector().push_back(
physicalEquations->returnOpCalculateStretchFromStress(
dataAtPts, physicalEquations));
} else {
fe->getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValues<3>(
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
}
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
// velocities
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
if (noStretch) {
} else {
fe->getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValuesDot<3>(
stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
}
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
// acceleration
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
}
// H1 displacements
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
// calculate other derived quantities
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
// evaluate integration points
if (noStretch) {
} else {
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setVolumeElementOps(
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
MoFEMFunctionBegin;
/** Contact requires that body is marked */
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % name).str()
))
) {
Range ents;
CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto rule_contact = [](int, int, int o) { return -1; };
auto refine = Tools::refineTriangle(contactRefinementLevels);
auto set_rule_contact = [refine](
ForcesAndSourcesCore *fe_raw_ptr, int order_row,
int order_col, int order_data
) {
MoFEMFunctionBegin;
auto rule = 2 * order_data;
fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
MoFEMFunctionReturn(0);
};
auto time_scale = boost::make_shared<TimeScale>();
// Right hand side
CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
// elastic
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialRotation(rotAxis, dataAtPts));
if (noStretch) {
// do nothing - no stretch approximation
} else {
fe_rhs->getOpPtrVector().push_back(
physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
alphaU));
}
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyP(piolaStress, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
auto set_hybridisation = [&](auto &pip) {
MoFEMFunctionBegin;
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
// First: Iterate over skeleton FEs adjacent to Domain FEs
// Note: BoundaryEle, i.e. uses skeleton interation rule
auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
CHKERR EshelbianPlasticity::
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
materialH1Positions, frontAdjEdges);
// Second: Iterate over domain FEs adjacent to skelton, particularly one
// domain element.
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
// Note: EleOnSide, i.e. uses on domain projected skeleton rule
auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(piolaStress,
flux_mat_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
// Assemble on skeleton
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dHybrid(hybridSpatialDisp, broken_data_ptr, 1.));
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
hybrid_ptr));
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
// Add skeleton to domain pipeline
pip.push_back(op_loop_skeleton_side);
MoFEMFunctionReturn(0);
};
auto set_contact = [&](auto &pip) {
MoFEMFunctionBegin;
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
// First: Iterate over skeleton FEs adjacent to Domain FEs
// Note: BoundaryEle, i.e. uses skeleton interation rule
auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
mField, contactElement, SPACE_DIM - 1, Sev::noisy);
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
CHKERR EshelbianPlasticity::
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
materialH1Positions, frontAdjEdges);
// Second: Iterate over domain FEs adjacent to skelton, particularly
// one domain element.
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
// Data storing contact fields
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
MoFEMFunctionBegin;
// Note: EleOnSide, i.e. uses on domain projected skeleton rule
auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
MoFEMFunctionReturn(0);
};
auto add_ops_contact_rhs = [&](auto &pip) {
MoFEMFunctionBegin;
// get body id and SDF range
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", SPACE_DIM - 1));
pip.push_back(new OpCalculateVectorFieldValues<3>(
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
pip.push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
nullptr));
pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
contactDisp, contact_common_data_ptr, contactTreeRhs,
contact_sfd_map_range_ptr));
pip.push_back(
new EshelbianPlasticity::OpConstrainBoundaryHDivRhs<A, GAUSS>(
broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
MoFEMFunctionReturn(0);
};
// push ops to face/side pipeline
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
// Add skeleton to domain pipeline
pip.push_back(op_loop_skeleton_side);
MoFEMFunctionReturn(0);
};
CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
CHKERR set_contact(fe_rhs->getOpPtrVector());
// Body forces
using BodyNaturalBC =
NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
Assembly<PETSC>::LinearForm<GAUSS>;
using OpBodyForce =
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {time_scale},
"BODY_FORCE", Sev::inform);
}
// Left hand side
CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
// elastic
if (add_elastic) {
const bool symmetric_system =
(gradApproximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
if (noStretch) {
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
bubbleField, piolaStress, dataAtPts));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
dataAtPts));
} else {
fe_lhs->getOpPtrVector().push_back(
physicalEquations->returnOpSpatialPhysical_du_du(
stretchTensor, stretchTensor, dataAtPts, alphaU));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
stretchTensor, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
stretchTensor, bubbleField, dataAtPts, true));
if (!symmetric_system) {
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
stretchTensor, rotAxis, dataAtPts, false));
}
}
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
spatialL2Disp, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
piolaStress, rotAxis, dataAtPts, symmetric_system));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
bubbleField, rotAxis, dataAtPts, symmetric_system));
if (!symmetric_system) {
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
rotAxis, piolaStress, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
rotAxis, bubbleField, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
}
auto set_hybridisation = [&](auto &pip) {
MoFEMFunctionBegin;
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
// First: Iterate over skeleton FEs adjacent to Domain FEs
// Note: BoundaryEle, i.e. uses skeleton interation rule
auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
CHKERR EshelbianPlasticity::
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
materialH1Positions, frontAdjEdges);
// Second: Iterate over domain FEs adjacent to skelton, particularly one
// domain element.
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
// Note: EleOnSide, i.e. uses on domain projected skeleton rule
auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC(hybridSpatialDisp, broken_data_ptr, 1., true, false));
pip.push_back(op_loop_skeleton_side);
MoFEMFunctionReturn(0);
};
auto set_contact = [&](auto &pip) {
MoFEMFunctionBegin;
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
using BdyEleOp = BoundaryEle::UserDataOperator;
// First: Iterate over skeleton FEs adjacent to Domain FEs
// Note: BoundaryEle, i.e. uses skeleton interation rule
auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
mField, contactElement, SPACE_DIM - 1, Sev::noisy);
op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
CHKERR EshelbianPlasticity::
AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
op_loop_skeleton_side->getOpPtrVector(), {L2},
materialH1Positions, frontAdjEdges);
// Second: Iterate over domain FEs adjacent to skelton, particularly
// one domain element.
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
// Data storing contact fields
auto contact_common_data_ptr =
boost::make_shared<ContactOps::CommonData>();
auto add_ops_domain_side = [&](auto &pip) {
MoFEMFunctionBegin;
// Note: EleOnSide, i.e. uses on domain projected skeleton rule
auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
broken_data_ptr));
op_loop_domain_side->getOpPtrVector().push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
MoFEMFunctionReturn(0);
};
auto add_ops_contact_lhs = [&](auto &pip) {
MoFEMFunctionBegin;
pip.push_back(new OpCalculateVectorFieldValues<3>(
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
pip.push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
nullptr));
// get body id and SDF range
auto contact_sfd_map_range_ptr =
boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", SPACE_DIM - 1));
pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
contact_sfd_map_range_ptr));
pip.push_back(
new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP<A, GAUSS>(
contactDisp, broken_data_ptr, contact_common_data_ptr,
contactTreeRhs, contact_sfd_map_range_ptr));
pip.push_back(
new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU<A, GAUSS>(
broken_data_ptr, contactDisp, contact_common_data_ptr,
contactTreeRhs));
MoFEMFunctionReturn(0);
};
// push ops to face/side pipeline
CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
// Add skeleton to domain pipeline
pip.push_back(op_loop_skeleton_side);
MoFEMFunctionReturn(0);
};
CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
CHKERR set_contact(fe_lhs->getOpPtrVector());
}
if (add_material) {
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setFaceElementOps(
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
MoFEMFunctionBegin;
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
// set integration rule
fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
if (add_elastic) {
auto get_broken_op_side = [this](auto &pip) {
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
// Iterate over domain FEs adjacent to boundary.
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
// Note: EleOnSide, i.e. uses on domain projected skeleton rule
auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
boost::make_shared<CGGUserPolynomialBase>();
CHKERR
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
pip.push_back(op_loop_domain_side);
return broken_data_ptr;
};
auto broken_data_ptr = get_broken_op_side(fe_rhs->getOpPtrVector());
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr,
{
boost::make_shared<TimeScale>("disp_history.txt")
}));
fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
broken_data_ptr, bcSpatialRotationVecPtr,
{
boost::make_shared<TimeScale>("rotation_history.txt")
}));
fe_rhs->getOpPtrVector().push_back(
new OpBrokenTractionBc(hybridSpatialDisp, bcSpatialTraction));
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setContactElementRhsOps(
boost::shared_ptr<ContactTree> &fe_contact_tree
) {
MoFEMFunctionBegin;
/** Contact requires that body is marked */
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % name).str()
))
) {
Range ents;
CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Skinner skin(&mField.get_moab());
for (auto &m : map) {
Range skin_faces;
CHKERR skin.find_skin(0, m.second, false, skin_faces);
CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
/* The above code is written in C++ and it appears to be defining and using
various operations on boundary elements and side elements. */
using BoundaryEle = FaceElementForcesAndSourcesCore;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto calcs_side_traction = [&](auto &pip) {
MoFEMFunctionBegin;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
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},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
MoFEMFunctionReturn(0);
};
auto add_contact_three = [&]() {
MoFEMFunctionBegin;
auto tree_moab_ptr = boost::make_shared<moab::Core>();
fe_contact_tree = boost::make_shared<ContactTree>(
mField, tree_moab_ptr, spaceOrder,
get_map_skin(get_body_range("BODY", 3)));
fe_contact_tree->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
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(
new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
fe_contact_tree->getOpPtrVector().push_back(
new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
MoFEMFunctionReturn(0);
};
CHKERR add_contact_three();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
MoFEMFunctionBegin;
// Add contact operators. Note that only for rhs. THe lhs is assembled with
// volume element, to enable schur complement evaluation.
CHKERR setContactElementRhsOps(contactTreeRhs);
CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_contact_bc = [&]() {
using SideEle = FaceElementForcesAndSourcesCore;
auto op_loop_side = new OpLoopSide<SideEle>(
mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
return op_loop_side;
};
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
MoFEMFunctionBegin;
boost::shared_ptr<FEMethod> null;
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
null);
CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
null);
CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
// CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
// null);
} else {
CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
null);
CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
null);
CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
// CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
// null);
}
MoFEMFunctionReturn(0);
}
#include "impl/EshelbianMonitor.cpp"
#include "impl/SetUpSchurImpl.cpp"
#include "impl/TSElasticPostStep.cpp"
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
MoFEMFunctionBegin;
#ifdef PYTHON_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";
CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
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";
}
#else
#endif
boost::shared_ptr<TsCtx> ts_ctx;
CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
ts_ctx->getLoopsMonitor().push_back(
TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
CHKERR TSAppendOptionsPrefix(ts, "elastic_");
CHKERR TSSetFromOptions(ts);
CHKERR TSSetDM(ts, dmElastic);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
CHKERR SNESMonitorSet(
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
PetscSection section;
CHKERR DMGetSection(dmElastic, &section);
int num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
const char *field_name;
CHKERR PetscSectionGetFieldName(section, ff, &field_name);
MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
}
CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
// Adding field split solver
boost::shared_ptr<SetUpSchur> schur_ptr;
if constexpr (A == AssemblyType::BLOCK_SCHUR) {
schur_ptr = SetUpSchur::createSetUpSchur(mField, &*this);
CHKERR schur_ptr->setUp(ts);
}
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
Vec xx;
CHKERR VecDuplicate(x, &xx);
CHKERR VecZeroEntries(xx);
CHKERR TS2SetSolution(ts, x, xx);
CHKERR VecDestroy(&xx);
} else {
CHKERR TSSetSolution(ts, x);
}
TetPolynomialBase::switchCacheBaseOn<HDIV>(
{elasticFeLhs.get(), elasticFeRhs.get()});
CHKERR TSSetUp(ts);
CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
CHKERR TSElasticPostStep::postStepInitialise(this);
CHKERR TSSolve(ts, PETSC_NULL);
CHKERR TSElasticPostStep::postStepDestroy();
TetPolynomialBase::switchCacheBaseOff<HDIV>(
{elasticFeLhs.get(), elasticFeRhs.get()});
// CHKERR TSGetSNES(ts, &snes);
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
MOFEM_LOG("EP", Sev::inform)
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
PETSC_NULL);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
if (lin_solver_iterations > expected_lin_solver_iterations)
SETERRQ2(
PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Expected number of iterations is different than expected %d > %d",
lin_solver_iterations, expected_lin_solver_iterations);
}
CHKERR gettingNorms();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
const std::string file) {
MoFEMFunctionBegin;
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_begin =
PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto get_post_proc = [&](auto sense) {
using FaceEle = FaceElementForcesAndSourcesCore;
auto post_proc_ptr =
boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
mField, post_proc_mesh);
EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
frontAdjEdges);
auto domain_ops = [&](auto &fe, int sense) {
MoFEMFunctionBegin;
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
frontAdjEdges);
fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
piolaStress, dataAtPts->getApproxPAtPts()));
fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
if (noStretch) {
fe.getOpPtrVector().push_back(
physicalEquations->returnOpCalculateStretchFromStress(
dataAtPts, physicalEquations));
} else {
fe.getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValues<3>(
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
}
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
// evaluate derived quantities
fe.getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
// evaluate integration points
fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, true, false, dataAtPts, physicalEquations));
if (auto op =
physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
fe.getOpPtrVector().push_back(op);
fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
}
// contact
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
// post-proc
fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
dataAtPts, sense));
MoFEMFunctionReturn(0);
};
post_proc_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(contactDisp,
dataAtPts->getContactL2AtPts()));
auto X_h1_ptr = boost::make_shared<MatrixDouble>();
// H1 material positions
post_proc_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(materialH1Positions,
dataAtPts->getLargeXH1AtPts()));
// domain
auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
mField, elementVolumeName, SPACE_DIM);
domain_ops(*(op_loop_side->getSideFEPtr()), sense);
post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
return post_proc_ptr;
};
auto post_proc_ptr = get_post_proc(1);
auto post_proc_negative_sense_ptr = get_post_proc(-1);
// contact
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
using BoundaryEle = FaceElementForcesAndSourcesCore;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
auto calcs_side_traction = [&](auto &pip) {
MoFEMFunctionBegin;
using EleOnSide =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy);
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},
materialH1Positions, frontAdjEdges);
op_loop_domain_side->getOpPtrVector().push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_domain_side);
MoFEMFunctionReturn(0);
};
CHKERR calcs_side_traction(post_proc_ptr->getOpPtrVector());
post_proc_ptr->getOpPtrVector().push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
&post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
CHKERR DMoFEMLoopFiniteElements(dM, skinElement, post_proc_ptr);
CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, crackElement, post_proc_ptr,
0, mField.get_comm_size());
CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, crackElement,
post_proc_negative_sense_ptr, 0,
mField.get_comm_size());
CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
CHKERR post_proc_end.writeFile(file.c_str());
MoFEMFunctionReturn(0);
}
//! [Getting norms]
MoFEMErrorCode EshelbianCore::gettingNorms() {
MoFEMFunctionBegin;
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(
post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
frontAdjEdges);
enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
auto norms_vec =
createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
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(
new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
u_h1_ptr));
auto piola_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
post_proc_norm_fe->getOpPtrVector().push_back(
new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
*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]);
MOFEM_LOG("EP", Sev::inform)
<< "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
MOFEM_LOG("EP", Sev::inform)
<< "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
MoFEMFunctionReturn(0);
}
//! [Getting norms]
MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
MoFEMFunctionBegin;
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
"", piolaStress, false, false);
bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto disp_bc = bc.second->dispBcPtr) {
MOFEM_LOG("EP", Sev::noisy) << *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);
bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
}
}
// old way of naming blocksets for displacement BCs
CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
MoFEMFunctionBegin;
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
false, false);
bcSpatialTraction = boost::make_shared<TractionBcVec>();
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);
bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
}
}
CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
MoFEMFunctionReturn(0);
}
} // namespace EshelbianPlasticity
#include <impl/EshelbianContact.cpp>
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:44
H1
@ H1
continuous field
Definition: definitions.h:85
TSElasticPostStep.cpp
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
EleOnSide
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
Definition: scalar_check_approximation.cpp:27
MoFEM::BaseFunction::DofsSideMap
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.
Definition: BaseFunction.hpp:73
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1071
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::TsCtx::getLoopsMonitor
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:98
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
TSElasticPostStep::postStepFun
static MoFEMErrorCode postStepFun(TS ts)
Definition: TSElasticPostStep.cpp:145
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2631
MoFEM::NaturalBC
Natural boundary conditions.
Definition: Natural.hpp:57
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
FaceRule::operator()
int operator()(int, int, int) const
Definition: simple_interface.cpp:92
_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(MFIELD, NAME, ENT, IT)
loop over all dofs from a moFEM field and particular field
Definition: Interface.hpp:1917
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianPlasticity::AddHOOps
Definition: EshelbianPlasticity.hpp:1150
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
EshelbianMonitor.cpp
Contains definition of EshelbianMonitor class.
QUAD_::order
int order
Definition: quad.h:28
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
MoFEM::ForceCubitBcData
Definition of the force bc data structure.
Definition: BCData.hpp:139
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::OpSetFlux
Definition: FormsBrokenSpaceConstraintImpl.hpp:139
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
EshelbianPlasticity::NO_H1_CONFIGURATION
@ NO_H1_CONFIGURATION
Definition: EshelbianPlasticity.hpp:43
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
SideEleOp
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
MoFEM::OpCalcNormL2Tensor2
Get norm of input MatrixDouble for Tensor2.
Definition: NormsOperators.hpp:72
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
MoFEM::OpGetBrokenBaseSideData
Definition: FormsBrokenSpaceConstraintImpl.hpp:68
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::OpBrokenLoopSide
Definition: FormsBrokenSpaceConstraintImpl.hpp:15
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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
Definition: DMMoFEM.cpp:853
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
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.
Definition: DMMoFEM.cpp:567
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
TSElasticPostStep::preStepFun
static MoFEMErrorCode preStepFun(TS ts)
Definition: TSElasticPostStep.cpp:95
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
save_range
static auto save_range(moab::Interface &moab, const std::string name, const Range r)
Definition: EshelbianPlasticity.cpp:74
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
MoFEM::DMMoFEMTSSetI2Jacobian
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
Definition: DMMoFEM.cpp:1017
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::DMMoFEMTSSetI2Function
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
Definition: DMMoFEM.cpp:975
MoFEM::CoreTmp
Definition: Core.hpp:36
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:461
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU
Definition: EshelbianContact.hpp:198
OpCalculateEshelbyStress
Definition: HookeElement.hpp:177
COL
@ COL
Definition: definitions.h:136
EshelbianPlasticity::OpConstrainBoundaryHDivRhs
Definition: EshelbianContact.hpp:134
EshelbianPlasticity::VolUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:41
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2572
EshelbianPlasticity::StretchSelector
StretchSelector
Definition: EshelbianPlasticity.hpp:44
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU
Definition: EshelbianContact.hpp:155
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MoFEM::PostProcBrokenMeshInMoabBaseEnd
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:952
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EshelbianPlasticity::CGG_BubbleBase_MBTET
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
Definition: CGGTonsorialBubbleBase.cpp:20
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
VolRule::operator()
int operator()(int, int, int) const
Definition: simple_interface.cpp:89
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MoFEM::BaseFunctionUnknownInterface
Definition: BaseFunction.hpp:13
MoFEM::DMMoFEMCreateMoFEM
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.
Definition: DMMoFEM.cpp:114
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
BiLinearForm
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FractureMechanics::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: CrackFrontElement.cpp:53
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:2760
TSElasticPostStep::postStepDestroy
static MoFEMErrorCode postStepDestroy()
Definition: TSElasticPostStep.cpp:85
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
cholesky.hpp
cholesky decomposition
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2692
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:30
Range
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:43
FTensor::dd
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)
Definition: ddTensor0.hpp:33
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
EshelbianPlasticity::OpConstrainBoundaryL2Rhs
Definition: EshelbianContact.hpp:118
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:44
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:43
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition: definitions.h:101
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EshelbianContact.cpp
MoFEM::DMMoFEMTSSetIFunction
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
Definition: DMMoFEM.cpp:800
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:977
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EshelbianContact.hpp
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
EshelbianMonitor
Definition: EshelbianMonitor.cpp:6
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
TSElasticPostStep::postStepInitialise
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
Definition: TSElasticPostStep.cpp:11
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition: test_broken_space.cpp:37
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
SetUpSchurImpl.cpp
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:768
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP
Definition: EshelbianContact.hpp:172
get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition: EshelbianPlasticity.cpp:52
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:43
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1104
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::PostProcBrokenMeshInMoabBaseBegin
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:934
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:17
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1977
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
quad.h
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow