Eshelbian plasticity implementation.
#include <boost/math/constants/constants.hpp>
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#pragma message "With PYTHON_SDF"
#else
#pragma message "Without PYTHON_SDF"
#endif
struct VolUserDataOperatorStabAssembly
: public VolumeElementForcesAndSourcesCore::UserDataOperator {
using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
: public FaceElementForcesAndSourcesCore::UserDataOperator {
using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
}
const std::string block_name) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
bc->getMeshsetIdEntitiesByDimension(m_field.
get_moab(), 2, faces,
true),
"get meshset ents");
r.merge(faces);
}
};
template <>
matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPB(), row_data, col_data,
m, ADD_VALUES);
};
return 0;
}
}
: mField(m_field), piolaStress("P"), eshelbyStress("S"),
spatialL2Disp("wL2"), spatialH1Disp("wH1"), contactDisp("contactDisp"),
materialL2Disp("W"), stretchTensor("u"), rotAxis("omega"),
materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
bubbleField("BUBBLE"), elementVolumeName("EP"),
naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC"),
skinElement("SKIN_ELEMENT"), contactElement("CONTACT_ELEMENT") {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
EshelbianCore::~EshelbianCore() = default;
const char *list_rots[] = {"small", "moderate", "large"};
PetscInt choice_rot = EshelbianCore::rotSelector;
PetscInt choice_grad = EshelbianCore::gradApperoximator;
const char *list_stretches[] = {"linear", "log"};
PetscInt choice_stretch = StretchSelector::LOG;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
"none");
spaceOrder = 2;
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
spaceOrder, &spaceOrder, PETSC_NULL);
spaceH1Order = -1;
CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
spaceH1Order, &spaceH1Order, PETSC_NULL);
materialOrder = 1;
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
"", materialOrder, &materialOrder, PETSC_NULL);
alphaU = 0;
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"", alphaU,
&alphaU, PETSC_NULL);
alphaW = 0;
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"", alphaW,
&alphaW, PETSC_NULL);
alphaRho = 0;
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"", alphaRho,
&alphaRho, PETSC_NULL);
precEps = 0;
CHKERR PetscOptionsScalar(
"-preconditioner_eps",
"preconditioner_eps",
"",
precEps, &precEps, PETSC_NULL);
precEpsOmega = 0;
CHKERR PetscOptionsScalar(
"-preconditioner_eps_omega",
"preconditioner_eps",
"", precEpsOmega, &precEpsOmega, PETSC_NULL);
precEpsW = 0;
CHKERR PetscOptionsScalar(
"-preconditioner_eps_w",
"preconditioner_eps",
"",
precEpsW, &precEpsW, PETSC_NULL);
precEpsContactDisp = 0;
CHKERR PetscOptionsScalar(
"-preconditioner_eps_contact_disp",
"preconditioner_eps", "", precEpsContactDisp,
&precEpsContactDisp, 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, LARGE_ROT + 1, list_rots[choice_grad],
&choice_grad, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-exponent_base",
"exponent_base",
"", exponentBase,
&EshelbianCore::exponentBase, PETSC_NULL);
"-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
ierr = PetscOptionsEnd();
EshelbianCore::rotSelector =
static_cast<RotSelector>(choice_rot);
EshelbianCore::gradApperoximator =
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:
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;
break;
default:
break;
};
precEpsOmega += precEps;
precEpsW += precEps;
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) <<
"precEps " << precEps;
MOFEM_LOG(
"EP", Sev::inform) <<
"precEpsOmega " << precEpsOmega;
MOFEM_LOG(
"EP", Sev::inform) <<
"precEpsW " << precEpsW;
<< "Rotations " << list_rots[EshelbianCore::rotSelector];
MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
<< list_rots[EshelbianCore::gradApperoximator];
if (exponentBase != exp(1))
<< "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;
}
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
if (bcSpatialDispVecPtr)
for (
auto &
v : *bcSpatialDispVecPtr) {
tets_skin = subtract(tets_skin,
v.faces);
}
if (bcSpatialRotationVecPtr)
for (
auto &
v : *bcSpatialRotationVecPtr) {
tets_skin = subtract(tets_skin,
v.faces);
}
if (bcSpatialTraction)
for (
auto &
v : *bcSpatialTraction) {
tets_skin = subtract(tets_skin,
v.faces);
}
auto subtract_faces_where_displacements_are_applied =
[&](const std::string block_name) {
tets_skin = subtract(tets_skin, contact_range);
};
CHKERR subtract_faces_where_displacements_are_applied(
"CONTACT");
CHKERR mField.get_moab().get_adjacencies(tets, 2,
true, faces,
moab::Interface::UNION);
Range faces_not_on_the_skin = subtract(faces, tets_skin);
};
auto add_l2_field = [
this, meshset](
const std::string
field_name,
};
auto add_h1_field = [
this, meshset](
const std::string
field_name,
};
auto add_l2_field_by_range = [
this](
const std::string
field_name,
const int field_dim,
Range &&
r) {
};
auto add_bubble_field = [
this, meshset](
const std::string
field_name,
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) {
};
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 add_hdiv_field(piolaStress, spaceOrder, 3);
CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
CHKERR add_l2_field(stretchTensor, spaceOrder, 6);
CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
}
EshelbianCore::addVolumeFiniteElement(
const EntityHandle meshset) {
auto add_field_to_fe = [this](const std::string fe,
};
if (!mField.check_finite_element(elementVolumeName)) {
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);
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, stretchTensor);
CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
CHKERR add_field_to_fe(elementVolumeName, contactDisp);
}
CHKERR mField.build_finite_elements(elementVolumeName);
}
EshelbianCore::addBoundaryFiniteElement(
const EntityHandle meshset) {
auto add_field_to_fe = [this](const std::string fe,
};
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);
}
}
Range essentail_bc_elements;
if (bcSpatialTraction) {
for (
auto &
v : *bcSpatialTraction) {
essentail_bc_elements.merge(
v.faces);
}
}
CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
naturalBcElement);
CHKERR add_field_to_fe(naturalBcElement, piolaStress);
CHKERR mField.build_finite_elements(naturalBcElement);
CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
essentialBcElement);
CHKERR add_field_to_fe(essentialBcElement, piolaStress);
CHKERR mField.build_finite_elements(essentialBcElement);
auto get_skin = [&]() {
CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto &&skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto skin = filter_true_skin(get_skin());
CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
CHKERR add_field_to_fe(skinElement, piolaStress);
CHKERR add_field_to_fe(skinElement, spatialH1Disp);
CHKERR add_field_to_fe(skinElement, contactDisp);
CHKERR mField.build_finite_elements(skinElement);
MOFEM_LOG(
"EP", Sev::inform) <<
"Contact elements " << contact_range.size();
CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
contactElement);
CHKERR add_field_to_fe(contactElement, piolaStress);
CHKERR add_field_to_fe(contactElement, spatialH1Disp);
CHKERR add_field_to_fe(contactElement, contactDisp);
CHKERR mField.build_finite_elements(contactElement);
}
dM =
createDM(mField.get_comm(),
"DMMOFEM");
BitRefLevel().set());
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
auto remove_dofs_on_essential_spatial_stress_boundary =
[&](const std::string prb_name) {
for (int d : {0, 1, 2})
prb_name, piolaStress, (*bcSpatialFreeTraction)[
d],
d,
d, 0,
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary(
"ESHELBY_PLASTICITY");
dmElastic =
createDM(mField.get_comm(),
"DMMOFEM");
"ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
"ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
"ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
"ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
"ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
"ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
"ELASTIC_PROBLEM", bubbleField, bubbleField);
"ELASTIC_PROBLEM", bubbleField, piolaStress);
"ELASTIC_PROBLEM", piolaStress, bubbleField);
{
PetscSection section;
§ion);
CHKERR DMSetSection(dmElastic, section);
CHKERR DMSetGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(§ion);
}
dmPrjSpatial =
createDM(mField.get_comm(),
"DMMOFEM");
->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
"PROJECT_SPATIAL", spatialH1Disp, true, false);
}
BcDisp::BcDisp(std::string name, std::vector<double> &attr,
Range &faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
<< "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " << faces.size();
}
BcRot::BcRot(std::string name, std::vector<double> &attr,
Range &faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
}
theta = attr[3];
}
TractionBc::TractionBc(std::string name, std::vector<double> &attr,
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
<< "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " << faces.size();
}
EshelbianCore::getTractionFreeBc(
const EntityHandle meshset,
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
bc_ptr->resize(3);
for (
int dd = 0;
dd != 3; ++
dd)
(*bc_ptr)[
dd] = tets_skin;
if (bcSpatialDispVecPtr)
for (
auto &
v : *bcSpatialDispVecPtr) {
(*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 (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);
}
std::regex((boost::format(
"%s(.*)") % contact_set_name).str()))) {
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);
}
}
int operator()(
int p_row,
int p_col,
int p_data)
const {
if (EshelbianCore::gradApperoximator > MODERATE_ROT)
return p_data + p_data + (p_data - 1);
else
return 2 * p_data;
}
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data; }
};
CGGUserPolynomialBase() {}
~CGGUserPolynomialBase() {}
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
CHKERR getValueHdivForCGGBubble(pts);
break;
default:
}
}
private:
MatrixDouble shapeFun;
boost::tuple<int, int, MatrixDouble> cachePhi = {0, 0, MatrixDouble()};
"Wrong base, should be USER_BASE");
}
const int nb_gauss_pts = pts.size2();
if(cachePhi.get<0>() ==
order && cachePhi.get<1>() == nb_gauss_pts) {
auto &mat = cachePhi.get<2>();
phi.resize(mat.size1(), mat.size2(),
false);
} else {
shapeFun.resize(nb_gauss_pts, 4, false);
&pts(2, 0), nb_gauss_pts);
double diff_shape_fun[12];
phi.resize(nb_gauss_pts, 9 * nb_base_functions,
false);
nb_gauss_pts);
cachePhi.get<0>() =
order;
cachePhi.get<1>() = nb_gauss_pts;
cachePhi.get<2>().resize(
phi.size1(),
phi.size2(),
false);
noalias(cachePhi.get<2>()) =
phi;
}
}
};
const int tag, const bool do_rhs, const bool do_lhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
if (!dataAtPts) {
dataAtPts =
dataAtPts->physicsPtr = physicalEquations;
}
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
piolaStress, dataAtPts->getDivPAtPts()));
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
fe->getOpPtrVector().push_back(
stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
}
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto time_scale = boost::make_shared<TimeScale>();
CHKERR setBaseVolumeElementOps(tag,
true,
false, fe_rhs);
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialRotation(rotAxis, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialPhysical(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));
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<
GAUSS>;
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
fe_rhs->getOpPtrVector(), mField, "w", {time_scale}, "BODY_FORCE",
Sev::inform);
}
CHKERR setBaseVolumeElementOps(tag,
true,
true, fe_lhs);
if (add_elastic) {
const bool symmetric_system =
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_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));
}
if constexpr (A == AssemblyType::SCHUR) {
using OpMassStab =
if (precEpsOmega > std::numeric_limits<double>::epsilon()) {
fe_lhs->getOpPtrVector().push_back(
new OpMassStab(rotAxis, rotAxis, [this](double, double, double) {
return precEpsOmega;
}));
}
if (std::abs(alphaRho + alphaW) <
std::numeric_limits<double>::epsilon()) {
fe_lhs->getOpPtrVector().push_back(new OpMassStab(
spatialL2Disp, spatialL2Disp,
[this](double, double, double) { return precEpsW; }));
}
}
if (add_material) {
}
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
{HDIV});
{HDIV});
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr,
{
boost::make_shared<TimeScale>("disp_history.txt")
}));
fe_rhs->getOpPtrVector().push_back(
new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
}
}
boost::shared_ptr<ContactTree> &fe_contact_tree,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_lhs
) {
auto get_body_range = [
this](
auto name,
int dim) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
fe_face_rhs = boost::make_shared<BoundaryEle>(mField);
fe_face_lhs = boost::make_shared<BoundaryEle>(mField);
auto rule = [](int, int,
int o) {
return -1; };
int levels = 0;
&levels, PETSC_NULL);
auto refine = Tools::refineTriangle(levels);
auto set_rule = [levels, refine](
int order_col, int order_data
) {
auto rule = 2 * order_data;
fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
};
fe_face_rhs->getRuleHook = rule;
fe_face_lhs->getRuleHook = rule;
fe_face_rhs->setRuleHook = set_rule;
fe_face_lhs->setRuleHook = set_rule;
fe_face_rhs->getOpPtrVector(), {HDIV});
fe_face_lhs->getOpPtrVector(), {HDIV});
auto add_contact_three = [&]() {
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)));
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
fe_contact_tree->getOpPtrVector(), {HDIV});
fe_contact_tree->getOpPtrVector().push_back(
contactDisp, contact_common_data_ptr->contactDispPtr()));
fe_contact_tree->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
fe_contact_tree->getOpPtrVector().push_back(
fe_contact_tree->getOpPtrVector().push_back(
new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
};
auto add_ops_face_lhs = [&](auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
contactDisp, contact_common_data_ptr->contactDispPtr()));
piolaStress, contact_common_data_ptr->contactTractionPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", 2));
contactDisp, contactDisp, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
contactDisp, piolaStress, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
piolaStress, contactDisp, contact_common_data_ptr, fe_contact_tree));
if constexpr (A == AssemblyType::SCHUR) {
using OpMassStab =
if (precEpsContactDisp > std::numeric_limits<double>::epsilon()) {
pip.push_back(new OpMassStab(
contactDisp, contactDisp,
[this](double, double, double) { return precEpsContactDisp; }));
}
}
};
auto add_ops_face_rhs = [&](auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
piolaStress, contact_common_data_ptr->contactTractionPtr()));
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", 2));
contactDisp, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
piolaStress, contact_common_data_ptr, fe_contact_tree));
};
CHKERR add_ops_face_lhs(fe_face_lhs->getOpPtrVector());
CHKERR add_ops_face_rhs(fe_face_rhs->getOpPtrVector());
}
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 = [&]() {
mField, contactElement,
SPACE_DIM - 1, Sev::noisy, adj_cache);
return op_loop_side;
};
auto op_contact_bc = get_op_contact_bc();
elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
contactTreeRhs, contactRhs, op_contact_bc->getSideFEPtr()
);
}
boost::shared_ptr<FEMethod> null;
boost::shared_ptr<FeTractionBc> spatial_traction_bc(
new FeTractionBc(mField, piolaStress, bcSpatialTraction));
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
null);
null);
null);
spatial_traction_bc);
null);
null);
} else {
null);
null);
null);
spatial_traction_bc);
null);
null);
}
}
#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;
};
if (file_exists("sdf.py")) {
MOFEM_LOG(
"EP", Sev::inform) <<
"sdf.py file found";
sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdf_python_ptr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
} else {
MOFEM_LOG(
"EP", Sev::warning) <<
"sdf.py file NOT found";
}
#else
#endif
boost::shared_ptr<TsCtx>
ts_ctx;
CHKERR TSMonitorSet(ts, TsMonitorSet,
ts_ctx.get(), PETSC_NULL);
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
CHKERR TSSetDM(ts, dmElastic);
SNES snes;
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
PetscSection section;
CHKERR DMGetSection(dmElastic, §ion);
int num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
}
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
boost::shared_ptr<SetUpSchur> schur_ptr;
if constexpr (A == AssemblyType::SCHUR) {
ts_ctx_ptr->zeroMatrix = false;
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetI2Jacobian(ts,
m,
p, PETSC_NULL, PETSC_NULL);
} else {
CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetIJacobian(ts,
m,
p, PETSC_NULL, PETSC_NULL);
}
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
elasticFeLhs->preProcessHook = [&]() {
*(elasticFeLhs->matAssembleSwitch) = false;
};
elasticFeLhs->postProcessHook = [&]() {
};
elasticBcLhs->preProcessHook = [&]() {
};
elasticBcLhs->preProcessHook = [&]() {
*(elasticBcLhs->matAssembleSwitch) = false;
};
}
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
Vec xx;
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
CHKERR TSSolve(ts, PETSC_NULL);
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
PETSC_NULL);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
}
}
const std::string file) {
if (!dataAtPts) {
dataAtPts =
}
{HDIV});
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto domain_ops = [&](auto &fe) {
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
fe.getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, true, false, dataAtPts, physicalEquations));
spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
};
CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
post_proc.getOpPtrVector().push_back(op_loop_side);
contactDisp, dataAtPts->getContactL2AtPts()));
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), dataAtPts));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc.getOpPtrVector().push_back(
post_proc.getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
post_proc.getOpPtrVector().push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
&post_proc.getMapGaussPts()));
CHKERR post_proc.writeFile(file.c_str());
}
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
auto post_proc_norm_rule_hook = [](int, int,
int p) ->
int {
};
post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
post_proc_norm_fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
auto norms_vec =
CHKERR VecZeroEntries(norms_vec);
auto u_l2_ptr = boost::make_shared<MatrixDouble>();
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
u_h1_ptr));
auto piola_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
CHKERR mField.loop_finite_elements(
"ELASTIC_PROBLEM", elementVolumeName,
*post_proc_norm_fe);
CHKERR VecAssemblyBegin(norms_vec);
CHKERR VecAssemblyEnd(norms_vec);
const double *norms;
CHKERR VecGetArrayRead(norms_vec, &norms);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
<< "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
<< "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
}
auto bc_mng = mField.getInterface<
BcManager>();
"", piolaStress, false, false);
bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto disp_bc = bc.second->dispBcPtr) {
std::vector<double> block_attributes(6, 0.);
if (disp_bc->data.flag1 == 1) {
block_attributes[0] = disp_bc->data.value1;
block_attributes[3] = 1;
}
if (disp_bc->data.flag2 == 1) {
block_attributes[1] = disp_bc->data.value2;
block_attributes[4] = 1;
}
if (disp_bc->data.flag3 == 1) {
block_attributes[2] = disp_bc->data.value3;
block_attributes[5] = 1;
}
auto faces = bc.second->bcEnts.subset_by_dimension(2);
bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
}
}
CHKERR getBc(bcSpatialDispVecPtr,
"SPATIAL_DISP_BC", 6);
}
auto bc_mng = mField.getInterface<
BcManager>();
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);
}
}
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Contains definition of EshelbianMonitor class.
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name)
Eshelbian plasticity interface.
#define MOFEM_LOG_C(channel, severity, format,...)
DomainNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
static PetscErrorCode ierr
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
#define MOFEM_LOG(channel, severity)
Log.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
const double v
phase velocity of light in medium (cm/ns)
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.
enum RotSelector EshelbianCore
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
constexpr auto field_name
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
static boost::function< double(const double)> dd_f
static enum RotSelector gradApperoximator
static boost::function< double(const double)> inv_f
static boost::function< double(const double)> inv_dd_f
static boost::function< double(const double)> inv_d_f
static double inv_f_log(const double v)
static double inv_d_f_log(const double v)
static double dd_f_log(const double v)
static enum StretchSelector stretchSelector
static double inv_dd_f_log(const double v)
static boost::function< double(const double)> d_f
static enum RotSelector rotSelector
EshelbianCore(MoFEM::Interface &m_field)
static double d_f_log(const double v)
static boost::function< double(const double)> f
static double exponentBase
static double f_log(const double v)
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Set integration rule to boundary elements.
int operator()(int, int, int) const
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark block DOFs.
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Definition of the displacement bc data structure.
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
default operator for TRI element
Provide data structure for (tensor) field approximation.
structure to get information form mofem into EntitiesFieldData
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Problem manager is used to build and partition problems.
intrusive_ptr for managing petsc objects
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)
int operator()(int, int, int) const