v0.9.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>
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);
alpha_u = 0;
CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alpha_u,
&alpha_u, PETSC_NULL);
alpha_w = 0;
CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alpha_w,
&alpha_w, PETSC_NULL);
preconditioner_eps = 1e-3;
CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
preconditioner_eps, &preconditioner_eps,
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 add_field_to_fe(elementVolumeName, materialDisp);
// CHKERR add_field_to_fe(elementVolumeName, tauField);
// CHKERR add_field_to_fe(elementVolumeName, lambdaField);
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());
// set integration rule
fe->getRuleHook = VolRule();
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->divPAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->wDotAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->logStreachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->diffStreachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->logStreachDotTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
dataAtPts->physicsPtr = physicalEquations;
}
// calculate fields values
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
piolaStress, dataAtPts->approxPAtPts));
fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
piolaStress, dataAtPts->divPAtPts));
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
eshelbyStress, dataAtPts->approxSigmaAtPts));
fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
eshelbyStress, dataAtPts->divSigmaAtPts));
fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
streachTensor, dataAtPts->logStreachTensorAtPts, MBTET));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
rotAxis, dataAtPts->rotAxisAtPts, MBTET));
fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
materialGradient, dataAtPts->GAtPts, MBTET));
fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
materialGradient + "0", dataAtPts->G0AtPts, MBTET));
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
spatialDisp, dataAtPts->wAtPts, MBTET));
// velocities
fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
spatialDisp, dataAtPts->wDotAtPts, MBTET));
fe->getOpPtrVector().push_back(
streachTensor, dataAtPts->logStreachDotTensorAtPts, 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, alpha_w));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialRotation(rotAxis, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialPhysical(streachTensor, dataAtPts, alpha_u));
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, alpha_u));
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, alpha_w));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
bubbleField, rotAxis, dataAtPts, false));
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));
// fe_lhs->getOpPtrVector().push_back(
// new OpSpatialRotation_domega_du(rotAxis, streachTensor, dataAtPts));
// Schur
dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
fe_lhs->getOpPtrVector().push_back(
new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
if (alpha_w < 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, preconditioner_eps));
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));
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);
schurAssembly = boost::make_shared<EpFEMethod>();
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);
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO aoS) {
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, Mat S, AO ao) {
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) {
Mat Suu;
CHKERR DMCreateMatrix(dmElasticSchurStreach, &Suu);
AO aoSuu;
CHKERR sub_data->getRowMap(&aoSuu);
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);
CHKERR MatDestroy(&Suu);
CHKERR AODestroy(&aoSuu);
const MoFEM::Problem *schur_bubble_prb_ptr;
CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
Mat SBubble;
CHKERR DMCreateMatrix(dmElasticSchurBubble, &SBubble);
AO aoSBubble;
CHKERR bubble_data->getRowMap(&aoSBubble);
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);
CHKERR MatDestroy(&SBubble);
CHKERR AODestroy(&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) {
Mat SOmega;
CHKERR DMCreateMatrix(dmElasticSchurOmega, &SOmega);
AO aoSOmega;
CHKERR tet_stress_data->getRowMap(&aoSOmega);
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) {
Mat Sw;
CHKERR DMCreateMatrix(dmElasticSchurSpatialDisp, &Sw);
AO aoSw;
CHKERR spatial_disp_data->getRowMap(&aoSw);
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);
CHKERR MatDestroy(&Sw);
CHKERR AODestroy(&aoSw);
} else
SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
"Problem does not have sub-problem data");
CHKERR MatDestroy(&SOmega);
CHKERR AODestroy(&aoSOmega);
} 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> volPostProcEle;
boost::shared_ptr<double> gEnergy;
Monitor(EshelbianCore &ep)
: eP(ep),
dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
->getData<VolEle>()),
volPostProcEle(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(
eP.dataAtPts->approxPAtPts));
vol_ele->getOpPtrVector().push_back(
eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
vol_ele->getOpPtrVector().push_back(
eP.streachTensor, eP.dataAtPts->logStreachTensorAtPts, MBTET));
vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
vol_ele->getOpPtrVector().push_back(
eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
eP.spatialDisp, eP.dataAtPts->wAtPts, MBTET));
vol_ele->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
eP.dataAtPts));
};
auto set_element_for_post_process = [&]() {
volPostProcEle->getRuleHook = VolRule();
volPostProcEle->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
volPostProcEle->getOpPtrVector().push_back(
eP.dataAtPts->approxPAtPts));
volPostProcEle->getOpPtrVector().push_back(
eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
volPostProcEle->getOpPtrVector().push_back(
eP.streachTensor, eP.dataAtPts->logStreachTensorAtPts, MBTET));
volPostProcEle->getOpPtrVector().push_back(
eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
volPostProcEle->getOpPtrVector().push_back(
eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
volPostProcEle->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(eP.spatialDisp,
eP.dataAtPts->wAtPts, MBTET));
volPostProcEle->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
eP.dataAtPts));
volPostProcEle->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;
};
PetscViewer viewer;
CHKERR PetscViewerBinaryOpen(
PETSC_COMM_WORLD, ("restart_" + get_str_time(ts_t) + ".dat").c_str(),
FILE_MODE_WRITE, &viewer);
CHKERR VecView(ts_u, viewer);
CHKERR PetscViewerDestroy(&viewer);
CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_str_time(ts_t) +
".h5m");
// Loop boundary elements with traction boundary conditions
*gEnergy = 0;
CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
*volPostProcEle,nullptr);
double body_energy;
MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
eP.mField.get_comm());
CHKERR PetscPrintf(eP.mField.get_comm(),
"Step %d time %3.4g strain energy %3.6e\n", 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();
};
std::ostringstream print;
print << add() << "comm rank " << eP.mField.get_comm_rank()
<< std::endl;
print << add() << "point " << getVectorAdaptor(point.data(), 3)
<< std::endl;
print << add() << "coords at gauss pts " << getCoordsAtGaussPts()
<< std::endl;
print << add() << "w " << (*eP.dataAtPts->wAtPts) << std::endl;
print << add() << "Piola " << (*eP.dataAtPts->approxPAtPts)
<< std::endl;
print << add() << "Cauchy " << t_cauchy << std::endl;
print << std::endl;
CHKERR PetscSynchronizedPrintf(eP.mField.get_comm(), "%s",
print.str().c_str());
}
}
}
};
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");
PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
CHKERR post_proc_at_points(pointB, "Point B");
PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
CHKERR post_proc_at_points(pointC, "Point C");
PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
}
};
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);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Field %d name %s\n", ff, field_name);
}
CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
// PetscRandom rctx;
// PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
// VecSetRandom(x, rctx);
// VecScale(x, -1e-1);
// PetscRandomDestroy(&rctx);
// CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
// CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
// 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);
}
}
CHKERR TSSolve(ts, x);
}
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
const std::string file) {
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->logStreachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->diffStreachTensorAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
}
PostProcVolumeOnRefinedMesh post_proc(mField);
post_proc.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
piolaStress, dataAtPts->approxPAtPts));
bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
post_proc.getOpPtrVector().push_back(
streachTensor, dataAtPts->logStreachTensorAtPts, MBTET));
rotAxis, dataAtPts->rotAxisAtPts, MBTET));
materialGradient, dataAtPts->GAtPts, MBTET));
spatialDisp, dataAtPts->wAtPts, 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