v0.14.0
Loading...
Searching...
No Matches
EshelbianPlasticity.cpp

Eshelbian plasticity implementation.

Eshelbian plasticity implementation

/**
* \file EshelbianPlasticity.cpp
* \example EshelbianPlasticity.cpp
*
* \brief Eshelbian plasticity implementation
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <boost/math/constants/constants.hpp>
#include <cholesky.hpp>
#ifdef PYTHON_SFD
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
#pragma message "With PYTHON_SFD"
#else
#pragma message "Without PYTHON_SFD"
#endif
struct VolUserDataOperatorStabAssembly
using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
} // namespace EshelbianPlasticity
template <>
SCHUR,
matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
const EntitiesFieldData::EntData &row_data,
const EntitiesFieldData::EntData &col_data,
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
};
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log;
boost::function<double(const double)> EshelbianCore::d_f =
boost::function<double(const double)> EshelbianCore::dd_f =
boost::function<double(const double)> EshelbianCore::inv_f =
boost::function<double(const double)> EshelbianCore::inv_d_f =
boost::function<double(const double)> EshelbianCore::inv_dd_f =
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
UnknownInterface **iface) const {
*iface = const_cast<EshelbianCore *>(this);
return 0;
}
if (evalRhs)
if (evalLhs)
}
: mField(m_field), piolaStress("P"), eshelbyStress("S"),
spatialL2Disp("wL2"), spatialH1Disp("wH1"), 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") {
ierr = getOptions();
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
EshelbianCore::~EshelbianCore() = default;
MoFEMErrorCode EshelbianCore::getOptions() {
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);
CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
PETSC_NULL);
CHKERR PetscOptionsEList("-grad", "gradient of defomation approximator", "",
list_rots, LARGE_ROT + 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);
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:
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
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;
MOFEM_LOG("EP", Sev::inform)
<< "Rotations " << list_rots[EshelbianCore::rotSelector];
MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
<< list_rots[EshelbianCore::gradApperoximator];
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;
}
MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
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);
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) {
auto msh_mng = mField.getInterface<MeshsetsManager>();
auto reg_exp = std::regex((boost::format("%s(.*)") % block_name).str());
for (auto m : msh_mng->getCubitMeshsetPtr(reg_exp)) {
Range faces;
CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(),
SPACE_DIM - 1, faces, true);
tets_skin = subtract(tets_skin, faces);
MOFEM_LOG("EP", Sev::inform)
<< "Subtracting " << m->getName() << " " << faces.size();
}
};
CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
Range faces;
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_hdiv_field = [&](const std::string field_name, const int order,
const int dim) {
MB_TAG_SPARSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
CHKERR mField.set_field_order(tets_skin, field_name, 0);
};
auto add_l2_field = [this, meshset](const std::string field_name,
const int order, const int 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);
};
auto add_h1_field = [this, meshset](const std::string field_name,
const int order, const int 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);
};
auto add_bubble_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
// 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) {
};
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);
};
// spatial fields
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);
// material fields
// CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
// CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
// CHKERR add_l2_field(materialL2Disp, materialOrder - 1, 3);
// CHKERR add_l2_field(tauField, materialOrder - 1, 1);
// CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
// Add history filedes
// CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
// spatial displacement
CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
CHKERR mField.build_fields();
}
EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
// set finite element fields
auto add_field_to_fe = [this](const std::string fe,
const std::string field_name) {
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);
};
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);
// CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
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, materialGradient);
// CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
// materialGradient +
// "0");
}
// build finite elements data structures
CHKERR mField.build_finite_elements(elementVolumeName);
}
EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
auto bc_elements_add_to_range = [&](const std::string block_name, Range &r) {
auto mesh_mng = mField.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
Range faces;
CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
true);
r.merge(faces);
}
};
// set finite element fields
auto add_field_to_fe = [this](const std::string fe,
const std::string field_name) {
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);
};
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_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, piolaStress);
// CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
CHKERR mField.build_finite_elements(naturalBcElement);
CHKERR mField.add_finite_element(essentialBcElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
essentialBcElement);
CHKERR add_field_to_fe(essentialBcElement, piolaStress);
// CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
CHKERR mField.build_finite_elements(essentialBcElement);
auto get_skin = [&]() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, 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 skin = filter_true_skin(get_skin());
CHKERR mField.add_finite_element(skinElement, MF_ZERO);
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, eshelbyStress);
CHKERR mField.build_finite_elements(skinElement);
Range contact_range;
CHKERR bc_elements_add_to_range("CONTACT", contact_range);
MOFEM_LOG("EP", Sev::inform) << "Contact elements " << contact_range.size();
CHKERR mField.add_finite_element(contactElement, MF_ZERO);
CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
contactElement);
CHKERR add_field_to_fe(contactElement, piolaStress);
CHKERR mField.build_finite_elements(contactElement);
}
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit) {
// 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 DMMoFEMAddElement(dM, elementVolumeName);
CHKERR DMMoFEMAddElement(dM, naturalBcElement);
CHKERR DMMoFEMAddElement(dM, essentialBcElement);
CHKERR DMMoFEMAddElement(dM, contactElement);
CHKERR DMMoFEMAddElement(dM, skinElement);
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
CHKERR DMSetUp(dM);
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})
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0,
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary("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.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp.c_str());
CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement);
CHKERR DMMoFEMAddElement(dmElastic, contactElement);
CHKERR DMMoFEMAddElement(dmElastic, skinElement);
CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
CHKERR DMSetUp(dmElastic);
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);
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);
{
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
&section);
CHKERR DMSetSection(dmElastic, section);
CHKERR DMSetGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(&section);
}
dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp.c_str());
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);
}
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();
}
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
// 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;
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);
}
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);
}
// 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);
}
// for (int dd = 0; dd != 3; ++dd) {
// EntityHandle meshset;
// CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
// CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
// std::string file_name = disp_block_set_name +
// "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
// CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
// &meshset, 1);
// CHKERR mField.get_moab().delete_entities(&meshset, 1);
// }
}
/**
* @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::gradApperoximator > MODERATE_ROT)
return p_data + p_data + (p_data - 1);
else
return 2 * p_data;
}
};
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const { return 2 * p_data; }
};
struct CGGUserPolynomialBase : public BaseFunction {
CGGUserPolynomialBase() {}
~CGGUserPolynomialBase() {}
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
case HDIV:
CHKERR getValueHdivForCGGBubble(pts);
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
}
}
private:
MatrixDouble shapeFun;
boost::tuple<int, int, MatrixDouble> cachePhi = {0, 0, MatrixDouble()};
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts) {
// This should be used only in case USER_BASE is selected
if (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 = 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();
if(cachePhi.get<0>() == order && cachePhi.get<1>() == nb_gauss_pts) {
auto &mat = cachePhi.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);
// direvatives 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
&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);
cachePhi.get<0>() = order;
cachePhi.get<1>() = nb_gauss_pts;
cachePhi.get<2>().resize(phi.size1(), phi.size2(), false);
noalias(cachePhi.get<2>()) = phi;
}
}
};
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
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});
// 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()));
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));
fe->getOpPtrVector().push_back(
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
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
}
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) {
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));
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));
// Body forces
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<GAUSS>;
using OpBodyForce =
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
fe_rhs->getOpPtrVector(), mField, "w", {time_scale}, "BODY_FORCE",
Sev::inform);
}
// Left hand side
CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
// elastic
if (add_elastic) {
const bool symmetric_system =
(gradApperoximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
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));
}
// Stabilisation
if constexpr (A == AssemblyType::SCHUR) {
// Note that we assemble to AMat, however Assembly<SCHUR> assemble by
// default to PMat
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) {
}
}
}
MoFEMErrorCode EshelbianCore::setFaceElementOps(
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);
// set integration rule
fe_rhs->getRuleHook = FaceRule();
fe_lhs->getRuleHook = FaceRule();
{HDIV});
{HDIV});
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr));
fe_rhs->getOpPtrVector().push_back(
new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
}
}
MoFEMErrorCode EshelbianCore::setContactElementOps(
boost::shared_ptr<ContactTree> &fe_contact_tree,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
/** Contact requires that body is marked */
auto get_body_range = [this](auto name) {
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(),
3, 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. */
fe_rhs = boost::make_shared<BoundaryEle>(mField);
fe_lhs = boost::make_shared<BoundaryEle>(mField);
fe_rhs->getRuleHook = FaceRule();
fe_lhs->getRuleHook = FaceRule();
{HDIV});
{HDIV});
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_side = [&](auto disp_ptr, auto col_ent_data_ptr) {
auto op_loop_side = new OpLoopSide<SideEle>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
op_loop_side->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, disp_ptr));
if (col_ent_data_ptr != nullptr) {
op_loop_side->getOpPtrVector().push_back(
spatialL2Disp, col_ent_data_ptr));
}
return op_loop_side;
};
auto add_ops_lhs = [&](auto &pip) {
auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(get_op_side(contact_common_data_ptr->contactDispPtr(),
col_ent_data_ptr));
piolaStress, contact_common_data_ptr->contactTractionPtr()));
piolaStress, contact_common_data_ptr, col_ent_data_ptr));
piolaStress, contact_common_data_ptr));
};
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")));
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto op_loop_side = new OpLoopSide<SideEle>(
mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
op_loop_side->getOpPtrVector().push_back(
spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
op_loop_side->getOpPtrVector().push_back(
col_ent_data_ptr));
fe_contact_tree->getOpPtrVector().push_back(op_loop_side);
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, col_ent_data_ptr));
};
auto add_ops_rhs = [&](auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
piolaStress, contact_common_data_ptr->contactTractionPtr()));
pip.push_back(
get_op_side(contact_common_data_ptr->contactDispPtr(), nullptr));
piolaStress, contact_common_data_ptr));
// pip.push_back(new OpTreeSearch(fe_contact_tree, contact_common_data_ptr,
// nullptr, nullptr));
};
CHKERR add_contact_three();
CHKERR add_ops_lhs(fe_lhs->getOpPtrVector());
CHKERR add_ops_rhs(fe_rhs->getOpPtrVector());
}
MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
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 = [&]() {
auto op_loop_side = new OpLoopSide<SideEle>(
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);
CHKERR setContactElementOps(contactTreeRhs, contactRhs,
op_contact_bc->getSideFEPtr());
}
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
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()) {
CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
null);
CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
null);
CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
null);
CHKERR DMMoFEMTSSetI2Function(dm, skinElement, contactTreeRhs, null, null);
CHKERR DMMoFEMTSSetI2Function(dm, contactElement, contactRhs, null, null);
spatial_traction_bc);
CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
null);
} else {
CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
null);
CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
null);
CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
null);
CHKERR DMMoFEMTSSetIFunction(dm, skinElement, contactTreeRhs, null, null);
CHKERR DMMoFEMTSSetIFunction(dm, contactElement, contactRhs, null, null);
spatial_traction_bc);
CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
null);
}
}
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Mat m, Vec f, Vec x) {
#ifdef PYTHON_SFD
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);
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::SCHUR) {
auto p = matDuplicate(m, MAT_DO_NOT_COPY_VALUES);
auto ts_ctx_ptr = getDMTsCtx(dmElastic);
ts_ctx_ptr->zeroMatrix = false;
// If density is larger than zero, use dynamic time solver
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;
CHKERR KSPGetPC(ksp, &pc);
mField, SmartPetscObj<Mat>(m, true), p, &*this);
CHKERR schur_ptr->setUp(ksp);
elasticFeLhs->preProcessHook = [&]() {
*(elasticFeLhs->matAssembleSwitch) = false;
CHKERR schur_ptr->preProc();
};
elasticFeLhs->postProcessHook = [&]() {
};
elasticBcLhs->preProcessHook = [&]() {
};
elasticBcLhs->preProcessHook = [&]() {
*(elasticBcLhs->matAssembleSwitch) = false;
CHKERR schur_ptr->postProc();
};
}
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);
}
CHKERR TSSetUp(ts);
CHKERR TSSolve(ts, PETSC_NULL);
// 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();
}
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
const std::string file) {
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
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});
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 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));
// contact
fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
};
mField, elementVolumeName, SPACE_DIM);
CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
post_proc.getOpPtrVector().push_back(op_loop_side);
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), spatialL2Disp,
dataAtPts));
// contact
post_proc.getOpPtrVector().push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, &post_proc.getPostProcMesh(),
&post_proc.getMapGaussPts()));
CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(),
contactTreeRhs.get());
CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(), &post_proc);
CHKERR post_proc.writeFile(file.c_str());
}
//! [Getting norms]
MoFEMErrorCode EshelbianCore::gettingNorms() {
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());
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 =
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));
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]);
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);
}
//! [Getting norms]
MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
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);
}
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
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);
}
} // namespace EshelbianPlasticity
static Index< 'p', 3 > p
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
Contains definition of EshelbianMonitor class.
Eshelbian plasticity interface.
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
DomainNaturalBC::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
static PetscErrorCode ierr
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
cholesky decomposition
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
@ MF_ZERO
Definition: definitions.h:98
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
const int dim
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
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
FTensor::Index< 'm', SPACE_DIM > m
static double phi
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
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:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
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:786
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:118
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:509
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
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:572
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:839
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:1003
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:961
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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.
auto bit
set bit
const double v
phase velocity of light in medium (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
const FTensor::Tensor2< T, Dim, Dim > Vec
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.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:424
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.
float d
Definition: sdf_hertz.py:5
int r
Definition: sdf.py:8
constexpr AssemblyType A
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 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
EshelbianCore(MoFEM::Interface &m_field)
static double d_f_log(const double v)
static boost::function< double(const double)> f
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.
Definition: BcManager.hpp:25
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
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
Provide data structure for (tensor) field approximation.
Definition of the force bc data structure.
Definition: BCData.hpp:139
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Definition: Natural.hpp:57
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.
Definition: TsCtx.hpp:98
base class for all interface classes
Base volume element used to integrate on skeleton.
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.
Definition: plastic.cpp:1482
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)
Set integration rule.
int operator()(int, int, int) const