v0.14.0
EshelbianPlasticity.cpp

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_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
};
struct FaceUserDataOperatorStabAssembly
};
} // namespace EshelbianPlasticity
const std::string block_name) {
auto mesh_mng = m_field.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
Range faces;
bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 2, faces, true),
"get meshset ents");
r.merge(faces);
}
return r;
};
template <>
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);
};
enum RotSelector EshelbianCore::rotSelector = LARGE_ROT;
enum RotSelector EshelbianCore::gradApperoximator = LARGE_ROT;
enum StretchSelector EshelbianCore::stretchSelector = LOG;
double EshelbianCore::exponentBase = exp(1);
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log;
boost::function<double(const double)> EshelbianCore::d_f =
EshelbianCore::d_f_log;
boost::function<double(const double)> EshelbianCore::dd_f =
EshelbianCore::dd_f_log;
boost::function<double(const double)> EshelbianCore::inv_f =
EshelbianCore::inv_f_log;
boost::function<double(const double)> EshelbianCore::inv_d_f =
EshelbianCore::inv_d_f_log;
boost::function<double(const double)> EshelbianCore::inv_dd_f =
EshelbianCore::inv_dd_f_log;
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
UnknownInterface **iface) const {
*iface = const_cast<EshelbianCore *>(this);
return 0;
}
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
if (evalRhs)
CHKERR evaluateRhs(data);
if (evalLhs)
CHKERR evaluateLhs(data);
}
: 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") {
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);
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);
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) {
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;
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 contact_range = get_range_from_block(mField, block_name);
tets_skin = subtract(tets_skin, contact_range);
};
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) {
CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, 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) {
CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
};
auto add_h1_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
};
auto add_l2_field_by_range = [this](const std::string field_name,
const int order, const int dim,
const int field_dim, Range &&r) {
CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
MB_TAG_SPARSE, MF_ZERO);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
CHKERR mField.set_field_order(r, field_name, order);
};
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);
CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
get_range_from_block(mField, "CONTACT"));
// 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, contactDisp);
// 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) {
// 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, contactDisp);
// CHKERR add_field_to_fe(skinElement, eshelbyStress);
CHKERR mField.build_finite_elements(skinElement);
auto contact_range = get_range_from_block(mField, "CONTACT");
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 add_field_to_fe(contactElement, spatialH1Disp);
CHKERR add_field_to_fe(contactElement, contactDisp);
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");
auto contact_range = get_range_from_block(mField, "CONTACT");
// Create elastic sub-problem
dmElastic = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
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);
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,
{
boost::make_shared<TimeScale>("disp_history.txt")
}));
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_face_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_lhs
) {
/** Contact requires that body is marked */
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % name).str()
))
) {
Range ents;
CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Skinner skin(&mField.get_moab());
for(auto &m : map) {
Range skin_faces;
CHKERR skin.find_skin(0, m.second, false, skin_faces);
CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
/* The above code is written in C++ and it appears to be defining and using
various operations on boundary elements and side elements. */
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;
CHKERR PetscOptionsGetInt(PETSC_NULL, "-contact_max_post_proc_ref_level",
&levels, PETSC_NULL);
auto refine = Tools::refineTriangle(levels);
auto set_rule = [levels, refine](
ForcesAndSourcesCore *fe_raw_ptr, int order_row,
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(
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));
};
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 OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
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));
// Stabilisation
if constexpr (A == AssemblyType::SCHUR) {
// Note that we assemble to AMat, however Assembly<SCHUR> assemble by
// default to PMat
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 OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
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_contact_three();
CHKERR add_ops_face_lhs(fe_face_lhs->getOpPtrVector());
CHKERR add_ops_face_rhs(fe_face_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, 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, 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_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);
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());
}
{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});
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 OpCalculateVectorFieldValues<3>(
contactDisp, dataAtPts->getContactL2AtPts()));
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), dataAtPts));
// contact
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc.getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
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,
get_range_from_block(mField, "CONTACT"), &post_proc.getPostProcMesh(),
&post_proc.getMapGaussPts()));
CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
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
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:21
H1
@ H1
continuous field
Definition: definitions.h:85
TSElasticPostStep.cpp
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NOISY
@ NOISY
Definition: definitions.h:211
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::BaseFunction::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, MoFEM::UnknownInterface **iface) const
Definition: BaseFunction.cpp:17
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1072
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::TsCtx::getLoopsMonitor
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:98
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
TSElasticPostStep::postStepFun
static MoFEMErrorCode postStepFun(TS ts)
Definition: TSElasticPostStep.cpp:135
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2635
MoFEM::NaturalBC
Natural boundary conditions.
Definition: Natural.hpp:57
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
FaceRule::operator()
int operator()(int, int, int) const
Definition: simple_interface.cpp:92
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
EshelbianMonitor.cpp
Contains definition of EshelbianMonitor class.
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name)
Definition: EshelbianPlasticity.cpp:46
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
VolRule
Set integration rule.
Definition: simple_interface.cpp:88
MoFEM::ForceCubitBcData
Definition of the force bc data structure.
Definition: BCData.hpp:139
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::BaseFunction
Base class if inherited used to calculate base functions.
Definition: BaseFunction.hpp:40
SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: adolc_plasticity.cpp:98
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
SetUpSchur::createSetUpSchur
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.
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
MoFEM::OpCalcNormL2Tensor2
Get norm of input MatrixDouble for Tensor2.
Definition: NormsOperators.hpp:72
MoFEM::BaseFunction::getValue
virtual MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: BaseFunction.cpp:24
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:853
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
TSElasticPostStep::preStepFun
static MoFEMErrorCode preStepFun(TS ts)
Definition: TSElasticPostStep.cpp:87
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
MoFEM::DMMoFEMTSSetI2Jacobian
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:1017
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::DMMoFEMTSSetI2Function
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:975
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::ProblemsManager::removeDofsOnEntities
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.
Definition: ProblemsManager.cpp:2605
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU
Definition: EshelbianContact.hpp:177
EshelbianPlasticity::OpConstrainBoundaryHDivRhs
Definition: EshelbianContact.hpp:129
EshelbianPlasticity::VolUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:35
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2576
EshelbianPlasticity::StretchSelector
StretchSelector
Definition: EshelbianPlasticity.hpp:21
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU
Definition: EshelbianContact.hpp:143
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EshelbianPlasticity::CGG_BubbleBase_MBTET
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
Definition: CGGTonsorialBubbleBase.cpp:20
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
VolRule::operator()
int operator()(int, int, int) const
Definition: simple_interface.cpp:89
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:234
EshelbianPlasticity::EshelbianCore
enum RotSelector EshelbianCore
Definition: EshelbianPlasticity.cpp:84
MoFEM::BaseFunctionUnknownInterface
Definition: BaseFunction.hpp:13
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
BiLinearForm
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:2764
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
TSElasticPostStep::postStepDestroy
static MoFEMErrorCode postStepDestroy()
Definition: TSElasticPostStep.cpp:77
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
cholesky.hpp
cholesky decomposition
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2696
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
EshelbianPlasticity::OpConstrainBoundaryL2Rhs
Definition: EshelbianContact.hpp:114
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:21
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:20
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EshelbianContact.cpp
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:800
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:978
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EshelbianContact.hpp
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
EshelbianMonitor
Definition: EshelbianMonitor.cpp:6
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
TSElasticPostStep::postStepInitialise
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
Definition: TSElasticPostStep.cpp:11
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
SetUpSchurImpl.cpp
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:575
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::SmartPetscObj< Mat >
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP
Definition: EshelbianContact.hpp:160
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1103
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1109
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1942
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182