v0.10.0
EshelbianPlasticity.cpp

Eshelbian plasticity implementation

/**
* \file EshelbianPlasticity.cpp
* \example EshelbianPlasticity.cpp
*
* \brief Eshelbian plasticity implementation
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <BasicFiniteElements.hpp>
#include <boost/math/constants/constants.hpp>
MoFEMErrorCode EshelbianCore::query_interface(const MOFEMuuid &uuid,
UnknownInterface **iface) const {
*iface = NULL;
*iface = const_cast<EshelbianCore *>(this);
}
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
}
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
if (type != MBVERTEX)
if (evalRhs)
CHKERR evaluateRhs(data);
if (evalLhs)
CHKERR evaluateLhs(data);
}
EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
: mField(m_field), piolaStress("P"), eshelbyStress("S"), spatialDisp("w"),
materialDisp("W"), streachTensor("u"), rotAxis("omega"),
materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
bubbleField("BUBBLE"), elementVolumeName("EP"),
naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC") {
ierr = getOptions();
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
EshelbianCore::~EshelbianCore() {}
MoFEMErrorCode EshelbianCore::getOptions() {
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
"none");
spaceOrder = 1;
CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
spaceOrder, &spaceOrder, 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 = 1e-3;
CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
precEps, &precEps, PETSC_NULL);
ierr = PetscOptionsEnd();
}
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);
auto subtract_faces_where_displacements_are_applied =
[&](const std::string disp_block_set_name) {
if (it->getName().compare(0, disp_block_set_name.length(),
disp_block_set_name) == 0) {
Range faces;
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
faces, true);
tets_skin = subtract(tets_skin, faces);
}
}
};
CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
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_hdiv_rt_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_DENSE, MF_ZERO);
CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
CHKERR mField.set_field_order(tets_skin, field_name, order);
};
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_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(spatialDisp, spaceOrder - 1, 3);
CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
CHKERR add_l2_field(streachTensor, spaceOrder, 6);
// material fields
CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
// CHKERR add_l2_field(materialDisp, 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);
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, streachTensor);
CHKERR add_field_to_fe(elementVolumeName, rotAxis);
CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
CHKERR add_field_to_fe(elementVolumeName, streachTensor);
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 disp_block_set_name,
Range &r) {
if (it->getName().compare(0, disp_block_set_name.length(),
disp_block_set_name) == 0) {
Range faces;
CHKERR it->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;
CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
Range essentail_bc_elements;
CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
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);
}
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit) {
// find adjacencies between finite elements and dofs
CHKERR mField.build_adjacencies(bit, QUIET);
// Create coupled problem
dM = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
BitRefLevel().set());
CHKERR DMMoFEMAddElement(dM, elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dM, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dM, essentialBcElement.c_str());
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, NOISY,
true);
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
// Create elastic sub-problem
dmElastic = createSmartDM(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, streachTensor.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialDisp.c_str());
CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement.c_str());
CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
CHKERR DMSetUp(dmElastic);
// Create elastic streach-problem
dmElasticSchurStreach = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElasticSchurStreach, dmElastic,
"ELASTIC_PROBLEM_STREACH_SCHUR");
CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurStreach, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, piolaStress.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, bubbleField.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, rotAxis.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, spatialDisp.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurStreach, elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurStreach, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurStreach, essentialBcElement.c_str());
CHKERR DMMoFEMSetSquareProblem(dmElasticSchurStreach, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurStreach, PETSC_TRUE);
CHKERR DMSetUp(dmElasticSchurStreach);
// Create elastic bubble-problem
dmElasticSchurBubble = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElasticSchurBubble, dmElasticSchurStreach,
"ELASTIC_PROBLEM_BUBBLE_SCHUR");
CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurBubble, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, piolaStress.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, rotAxis.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, spatialDisp.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurBubble, elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurBubble, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurBubble, essentialBcElement.c_str());
CHKERR DMMoFEMSetSquareProblem(dmElasticSchurBubble, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurBubble, PETSC_TRUE);
CHKERR DMSetUp(dmElasticSchurBubble);
// Create elastic omega-problem
dmElasticSchurOmega = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElasticSchurOmega, dmElasticSchurBubble,
"ELASTIC_PROBLEM_OMEGA_SCHUR");
CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurOmega, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, piolaStress.c_str());
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, spatialDisp.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurOmega, elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurOmega, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurOmega, essentialBcElement.c_str());
CHKERR DMMoFEMSetSquareProblem(dmElasticSchurOmega, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurOmega, PETSC_TRUE);
CHKERR DMSetUp(dmElasticSchurOmega);
// Create elastic tet_stress-problem
dmElasticSchurSpatialDisp = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dmElasticSchurSpatialDisp, dmElasticSchurOmega,
"ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurSpatialDisp, piolaStress.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
elementVolumeName.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
essentialBcElement.c_str());
CHKERR DMMoFEMSetSquareProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurSpatialDisp, PETSC_TRUE);
CHKERR DMSetUp(dmElasticSchurSpatialDisp);
{
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
&section);
CHKERR DMSetDefaultSection(dmElastic, section);
CHKERR DMSetDefaultGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(&section);
}
}
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]);
}
}
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]);
}
}
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string disp_block_set_name,
const std::string rot_block_set_name) {
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 (it->getName().compare(0, disp_block_set_name.length(),
disp_block_set_name) == 0) {
std::vector<double> block_attributes;
CHKERR it->getAttributes(block_attributes);
if (block_attributes.size() != 6) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"In block %s six attributes are required for given BC "
"blockset (3 values + "
"3 flags) != %d",
it->getName().c_str(), block_attributes.size());
}
Range faces;
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
true);
if (block_attributes[3] != 0)
(*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
if (block_attributes[4] != 0)
(*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
if (block_attributes[5] != 0)
(*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
}
if (it->getName().compare(0, rot_block_set_name.length(),
rot_block_set_name) == 0) {
Range faces;
CHKERR it->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 {
return 2 * (p_data + 1);
}
};
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(const MOFEMuuid &uuid,
*iface = NULL;
if (uuid == IDD_TET_BASE_FUNCTION) {
*iface = const_cast<CGGUserPolynomialBase *>(this);
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
}
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
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;
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
// Get approximation order on element. Note that bubble functions are only
// on tetrahedron.
const int order = data.dataOnEntities[MBTET][0].getDataOrder();
/// number of integration points
const int nb_gauss_pts = pts.size2();
// number of base functions
// 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);
}
};
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
const int tag, const bool do_rhs, const bool do_lhs,
boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
fe->getOpPtrVector().push_back(new OpL2Transform());
// 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 OpCalculateHVecTensorField<3, 3>(
eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
eshelbyStress, dataAtPts->getDivSigmaAtPts()));
fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
materialGradient, dataAtPts->getBigGAtPts(), MBTET));
fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
// velocities
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
fe->getOpPtrVector().push_back(
streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), 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>(
spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
}
// calculate other derived quantities
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
// evaluate integration points
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
}
MoFEMErrorCode EshelbianCore::setGenericVolumeElementOps(
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
// Right hand side
CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
// elastic
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpSpatialEquilibrium(spatialDisp, dataAtPts, alphaW, alphaRho));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialRotation(rotAxis, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialPhysical(streachTensor, 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));
}
// Left hand side
CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
// elastic
if (add_elastic) {
// Schur
fe_lhs->getOpPtrVector().push_back(
new OpSpatialSchurBegin(spatialDisp, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
streachTensor, streachTensor, dataAtPts, alphaU));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
streachTensor, bubbleField, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
spatialDisp, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
spatialDisp, spatialDisp, dataAtPts, alphaW, alphaRho));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
bubbleField, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
streachTensor, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
streachTensor, rotAxis, dataAtPts, false));
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));
// Schur
dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
fe_lhs->getOpPtrVector().push_back(
new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
if (alphaW < std::numeric_limits<double>::epsilon() &&
alphaRho < std::numeric_limits<double>::epsilon()) {
dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
fe_lhs->getOpPtrVector().push_back(
new OpSpatialPreconditionMass(spatialDisp, dataAtPts->wwMatPtr));
} else {
dataAtPts->wwMatPtr.reset();
}
fe_lhs->getOpPtrVector().push_back(
new OpSpatialSchurEnd(spatialDisp, dataAtPts, precEps));
if (add_material) {
}
}
}
MoFEMErrorCode EshelbianCore::setGenericFaceElementOps(
const bool add_elastic, const bool add_material,
boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
fe_rhs =
boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
fe_lhs =
boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
// set integration rule
fe_rhs->getRuleHook = FaceRule();
fe_lhs->getRuleHook = FaceRule();
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::setElasticElementOps(const int tag) {
CHKERR setGenericVolumeElementOps(tag, true, false, elasticFeRhs,
elasticFeLhs);
CHKERR setGenericFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
}
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
boost::shared_ptr<FEMethod> null;
boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
new EpElement<FeTractionBc>(mField, piolaStress, bcSpatialTraction));
schurAssembly = boost::make_shared<EpFEMethod>();
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);
spatial_traction_bc);
CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
null);
CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
} 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);
spatial_traction_bc);
CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
null);
CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
null);
CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
}
}
MoFEMErrorCode EshelbianCore::setUpTSElastic(TS ts, Mat m, Vec f) {
boost::shared_ptr<TsCtx> ts_ctx;
CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetI2Jacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
} else {
CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
}
CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
auto add_schur_streach_op = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addStreachSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_streach_pre = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addStreachSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_bubble_op = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addBubbleSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_bubble_pre = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addBubbleSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_spatial_disp_op = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_spatial_disp_pre = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_omega_op = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addOmegaSchurMatrix(S, aoS);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
auto add_schur_omega_pre = [this](auto &list, SmartPetscObj<Mat> S,
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addOmegaSchurMatrix(S, ao);
else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
};
const MoFEM::Problem *schur_streach_prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
CHKERR DMCreateMatrix_MoFEM(dmElasticSchurStreach, Suu);
SmartPetscObj<AO> aoSuu = sub_data->getSmartRowMap();
CHKERR add_schur_streach_op(ts_ctx->loops_to_do_IJacobian, Suu, aoSuu);
CHKERR add_schur_streach_pre(ts_ctx->preProcess_IJacobian, Suu, aoSuu);
CHKERR add_schur_streach_pre(ts_ctx->postProcess_IJacobian, Suu, aoSuu);
const MoFEM::Problem *schur_bubble_prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
CHKERR DMCreateMatrix_MoFEM(dmElasticSchurBubble, SBubble);
SmartPetscObj<AO> aoSBubble = bubble_data->getSmartRowMap();
CHKERR add_schur_bubble_op(ts_ctx->loops_to_do_IJacobian, SBubble,
aoSBubble);
CHKERR add_schur_bubble_pre(ts_ctx->preProcess_IJacobian, SBubble,
aoSBubble);
CHKERR add_schur_bubble_pre(ts_ctx->postProcess_IJacobian, SBubble,
aoSBubble);
const MoFEM::Problem *schur_omega_prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
CHKERR DMCreateMatrix_MoFEM(dmElasticSchurOmega, SOmega);
SmartPetscObj<AO> aoSOmega = tet_stress_data->getSmartRowMap();
CHKERR add_schur_omega_op(ts_ctx->loops_to_do_IJacobian, SOmega,
aoSOmega);
CHKERR add_schur_omega_pre(ts_ctx->preProcess_IJacobian, SOmega,
aoSOmega);
CHKERR add_schur_omega_pre(ts_ctx->postProcess_IJacobian, SOmega,
aoSOmega);
const MoFEM::Problem *schur_spatial_disp_prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
&schur_spatial_disp_prb_ptr);
if (auto spatial_disp_data =
schur_spatial_disp_prb_ptr->subProblemData) {
CHKERR DMCreateMatrix_MoFEM(dmElasticSchurSpatialDisp, Sw);
SmartPetscObj<AO> aoSw = spatial_disp_data->getSmartRowMap();
CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
aoSw);
} else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
"Problem does not have sub-problem data");
} else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
"Problem does not have sub-problem data");
} else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
"Problem does not have sub-problem data");
} else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
"Problem does not have sub-problem data");
struct Monitor : public FEMethod {
EshelbianCore &eP;
boost::shared_ptr<SetPtsData> dataFieldEval;
boost::shared_ptr<VolEle> volPostProcEnergy;
boost::shared_ptr<double> gEnergy;
Monitor(EshelbianCore &ep)
: eP(ep),
dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
->getData<VolEle>()),
volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
dataFieldEval, "EP");
CHKERRABORT(PETSC_COMM_WORLD, ierr);
auto no_rule = [](int, int, int) { return -1; };
auto set_element_for_field_eval = [&]() {
boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
vol_ele->getRuleHook = no_rule;
vol_ele->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
vol_ele->getOpPtrVector().push_back(new OpL2Transform());
vol_ele->getOpPtrVector().push_back(
eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
vol_ele->getOpPtrVector().push_back(
eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
vol_ele->getOpPtrVector().push_back(
eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
MBTET));
vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
vol_ele->getOpPtrVector().push_back(
eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
vol_ele->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
eP.dataAtPts));
};
auto set_element_for_post_process = [&]() {
volPostProcEnergy->getRuleHook = VolRule();
volPostProcEnergy->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
volPostProcEnergy->getOpPtrVector().push_back(
eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
volPostProcEnergy->getOpPtrVector().push_back(
eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
volPostProcEnergy->getOpPtrVector().push_back(
eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
MBTET));
volPostProcEnergy->getOpPtrVector().push_back(
eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
volPostProcEnergy->getOpPtrVector().push_back(
eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
volPostProcEnergy->getOpPtrVector().push_back(
eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
volPostProcEnergy->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
eP.dataAtPts));
volPostProcEnergy->getOpPtrVector().push_back(
new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
};
set_element_for_field_eval();
set_element_for_post_process();
}
MoFEMErrorCode preProcess() { return 0; }
MoFEMErrorCode operator()() { return 0; }
MoFEMErrorCode postProcess() {
// auto get_str_time = [](auto ts_t) {
// std::ostringstream ss;
// ss << boost::str(boost::format("%d") %
// static_cast<int>(std::ceil(ts_t * 1e6)));
// std::string s = ss.str();
// return s;
// };
auto get_step = [](auto ts_step) {
std::ostringstream ss;
ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
std::string s = ss.str();
return s;
};
PetscViewer viewer;
CHKERR PetscViewerBinaryOpen(
PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
FILE_MODE_WRITE, &viewer);
CHKERR VecView(ts_u, viewer);
CHKERR PetscViewerDestroy(&viewer);
CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
".h5m");
// Loop boundary elements with traction boundary conditions
*gEnergy = 0;
CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
*volPostProcEnergy);
double body_energy;
MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
eP.mField.get_comm());
MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
ts_step, ts_t, body_energy);
auto post_proc_at_points = [&](std::array<double, 3> point,
std::string str) {
dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
struct OpPrint : public VolOp {
EshelbianCore &eP;
std::array<double, 3> point;
std::string str;
OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
std::string &str)
: VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
str(str) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBTET) {
if (getGaussPts().size2()) {
auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
auto t_approx_P =
getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
const double jac = dEterminant(t_h);
t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
auto add = [&]() {
std::ostringstream s;
s << str << " elem " << getFEEntityHandle() << " ";
return s.str();
};
auto print_tensor = [](auto &t) {
std::ostringstream s;
s << t;
return s.str();
};
std::ostringstream print;
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "comm rank " << eP.mField.get_comm_rank();
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "point " << getVectorAdaptor(point.data(), 3);
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "coords at gauss pts " << getCoordsAtGaussPts();
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "w " << eP.dataAtPts->wAtPts;
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "Piola " << eP.dataAtPts->approxPAtPts;
MOFEM_LOG("EPSYNC", Sev::inform)
<< add() << "Cauchy " << print_tensor(t_cauchy);
}
}
}
};
if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
"EP", dataFieldEval,
eP.mField.get_comm_rank(),
eP.mField.get_comm_rank(), MF_EXIST, QUIET);
fe_ptr->getOpPtrVector().pop_back();
}
};
// Points for Cook beam
std::array<double, 3> pointA = {48., 60., 5.};
CHKERR post_proc_at_points(pointA, "Point A");
MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
CHKERR post_proc_at_points(pointB, "Point B");
MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
CHKERR post_proc_at_points(pointC, "Point C");
MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
}
};
boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
ts_ctx->get_loops_to_do_Monitor().push_back(
TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
CHKERR TSAppendOptionsPrefix(ts, "elastic_");
CHKERR TSSetFromOptions(ts);
}
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
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 DMGetDefaultSection(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
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_uu_field_split;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
if (is_uu_field_split) {
const MoFEM::Problem *schur_uu_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_uu_ptr);
if (auto uu_data = schur_uu_ptr->subProblemData) {
const MoFEM::Problem *prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElastic, &prb_ptr);
map<std::string, IS> is_map;
for (int ff = 0; ff != num_fields; ff++) {
const char *field_name;
CHKERR PetscSectionGetFieldName(section, ff, &field_name);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
&is_map[field_name]);
}
// CHKERR mField.getInterface<ISManager>()
// ->isCreateProblemFieldAndEntityType(
// prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
// MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
// CHKERR mField.getInterface<ISManager>()
// ->isCreateProblemFieldAndEntityType(
// prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
// MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
schurAssembly->Suu);
CHKERR PCSetUp(pc);
PetscInt n;
KSP *uu_ksp;
CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
PC bubble_pc;
CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
PetscBool is_bubble_field_split;
PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
&is_bubble_field_split);
if (is_bubble_field_split) {
const MoFEM::Problem *schur_bubble_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
if (auto bubble_data = schur_bubble_ptr->subProblemData) {
CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
AO uu_ao;
CHKERR uu_data->getRowMap(&uu_ao);
CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
CHKERR PCFieldSplitSetSchurPre(
bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
CHKERR PCSetUp(bubble_pc);
PetscInt bubble_n;
KSP *bubble_ksp;
CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
PC omega_pc;
CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
PetscBool is_omega_field_split;
PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
&is_omega_field_split);
if (is_omega_field_split) {
const MoFEM::Problem *schur_omega_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
if (auto omega_data = schur_omega_ptr->subProblemData) {
AO bubble_ao;
CHKERR bubble_data->getRowMap(&bubble_ao);
CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
CHKERR PCFieldSplitSetSchurPre(omega_pc,
PC_FIELDSPLIT_SCHUR_PRE_USER,
schurAssembly->SOmega);
CHKERR PCSetUp(omega_pc);
PetscInt omega_n;
KSP *omega_ksp;
CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
PC w_pc;
CHKERR KSPGetPC(omega_ksp[1], &w_pc);
PetscBool is_w_field_split;
PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
&is_w_field_split);
if (is_w_field_split) {
const MoFEM::Problem *schur_w_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
&schur_w_ptr);
if (auto w_data = schur_w_ptr->subProblemData) {
AO omega_ao;
CHKERR omega_data->getRowMap(&omega_ao);
CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
CHKERR PCFieldSplitSetSchurPre(
w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
CHKERR AODestroy(&omega_ao);
}
}
CHKERR PetscFree(omega_ksp);
}
}
CHKERR PetscFree(bubble_ksp);
CHKERR AODestroy(&uu_ao);
}
}
CHKERR PetscFree(uu_ksp);
for (auto &m : is_map)
CHKERR ISDestroy(&m.second);
}
}
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 TSSolve(ts, PETSC_NULL);
}
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
const std::string file) {
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
PostProcVolumeOnRefinedMesh post_proc(mField);
post_proc.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
post_proc.getOpPtrVector().push_back(new OpL2Transform());
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
post_proc.getOpPtrVector().push_back(
streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
materialGradient, dataAtPts->getBigGAtPts(), MBTET));
spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
// evaluate derived quantities
post_proc.getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
// evaluate integration points
post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
spatialDisp, tag, true, false, dataAtPts, physicalEquations));
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
CHKERR post_proc.writeFile(file.c_str());
}
} // namespace EshelbianPlasticity
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field valuse for given petsc vector.
Definition: UserDataOperators.hpp:358
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:34
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEM::IDD_TET_BASE_FUNCTION
static const MOFEMuuid IDD_TET_BASE_FUNCTION
Definition: EntPolynomialBaseCtx.hpp:28
FTensor::Tensor2
Definition: Tensor2_value.hpp:17
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
EntityHandle
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
H1
@ H1
continuous field
Definition: definitions.h:177
MoFEM::VolumeElementForcesAndSourcesCoreSwitch::UserDataOperator
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:344
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into DataForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:34
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
NOISY
@ NOISY
Definition: definitions.h:280
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:130
n
static Index< 'n', 3 > n
Definition: BasicFeTools.hpp:78
MoFEM::BaseFunction
Base class if inherited used to calculate base functions.
Definition: BaseFunction.hpp:53
MoFEM::EntPolynomialBaseCtx::dAta
DataForcesAndSourcesCore & dAta
Definition: EntPolynomialBaseCtx.hpp:57
FTensor::d
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMMoFEM.cpp:378
L2
@ L2
field with C-1 continuity
Definition: definitions.h:180
MoFEM::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
Definition: VolumeElementForcesAndSourcesCore.hpp:354
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
EshelbianPlasticity::IDD_MOFEMEshelbianCrackInterface
static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
Definition: EshelbianPlasticity.hpp:17
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:607
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
MoFEM::VolumeElementForcesAndSourcesCoreSwitch
Volume finite element with switches.
Definition: VolumeElementForcesAndSourcesCore.hpp:339
VolEle
VolumeElementForcesAndSourcesCore VolEle
Definition: test_cache_on_entities.cpp:34
MoFEM::Problem::getName
std::string getName() const
Definition: ProblemsMultiIndices.hpp:527
MoFEM::BaseFunctionUnknownInterface
Definition: BaseFunction.hpp:28
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:1605
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: DMMMoFEM.cpp:445
dim
const int dim
Definition: elec_phys_2D.cpp:12
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:70
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:160
ROW
@ ROW
Definition: definitions.h:192
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:315
MoFEM::ForcesAndSourcesCore::getUserPolynomialBase
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:439
m
static Index< 'm', 3 > m
Definition: BasicFeTools.hpp:77
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1053
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEM::BaseFunction::query_interface
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::BaseFunctionUnknownInterface **iface) const
Definition: BaseFunction.cpp:35
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:198
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
MoFEM::FieldEvaluatorInterface::getData
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
Definition: FieldEvaluator.hpp:118
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: DMMMoFEM.cpp:772
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:776
HenkyOps::f
auto f
Definition: HenkyOps.hpp:19
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:239
Monitor
[Push operators to pipeline]
Definition: dynamic_elastic.cpp:295
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: DMMMoFEM.cpp:866
MoFEM::createSmartDM
auto createSmartDM
Creates smart DM object.
Definition: PetscSmartObj.hpp:129
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: DMMMoFEM.cpp:844
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: DMMMoFEM.cpp:177
convert.type
type
Definition: convert.py:66
MoFEM.hpp
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:360
FTensor::Index< 'i', 3 >
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:714
StdRDOperators::order
const int order
Definition: rd_stdOperators.hpp:28
MoFEM::FieldEvaluatorInterface::SetPtsData
Definition: FieldEvaluator.hpp:43
MoFEM::MOFEMuuid
MoFEM interface unique ID.
Definition: UnknownInterface.hpp:26
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:1664
MoFEM::SmartPetscObj< Mat >
VolRule
Set integration rule.
Definition: simple_interface.cpp:103
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:34
FaceElementForcesAndSourcesCore
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:106
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:642
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: DMMMoFEM.cpp:105
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:36
VolOp
VolEle::UserDataOperator VolOp
Definition: test_cache_on_entities.cpp:35
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:131
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:425
MoFEM::DataForcesAndSourcesCore
data structure for finite element entity
Definition: DataStructures.hpp:52
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:305
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:281
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:64
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:349
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
Definition: ForcesAndSourcesCore.hpp:99
MF_ZERO
@ MF_ZERO
Definition: definitions.h:189
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:91
MoFEM::DataForcesAndSourcesCore::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: DataStructures.hpp:1040
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:80
BLOCKSET
@ BLOCKSET
Definition: definitions.h:217
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:1724
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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: DMMMoFEM.cpp:719
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
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
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1001
point
bg::model::point< double, 4, bg::cs::cartesian > point
Definition: r_tree_testing.hpp:13
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:209
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:342
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:24
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:44
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:60
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
QUIET
@ QUIET
Definition: definitions.h:277
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:79
MoFEM::VolumeElementForcesAndSourcesCoreBase::UserDataOperator
default operator for TET element
Definition: VolumeElementForcesAndSourcesCore.hpp:45
MoFEM::dEterminant
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:875
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MF_EXIST
@ MF_EXIST
Definition: definitions.h:189
convert.int
int
Definition: convert.py:66
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:179
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
MoFEM::Problem::subProblemData
boost::shared_ptr< SubProblemData > subProblemData
Definition: ProblemsMultiIndices.hpp:228
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
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, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2856
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:43
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:84
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
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
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:318
DataAtIntegrationPts
Definition: HookeElement.hpp:92
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:59