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

Eshelbian plasticity implementation.

Eshelbian plasticity implementation

/**
* \file EshelbianPlasticity.cpp
* \example EshelbianPlasticity.cpp
*
* \brief Eshelbian plasticity implementation
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <boost/math/constants/constants.hpp>
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
UnknownInterface **iface) const {
*iface = const_cast<EshelbianCore *>(this);
return 0;
}
if (type != MBVERTEX)
if (evalRhs)
if (evalLhs)
}
: 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);
}
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
spaceOrder, &spaceOrder, PETSC_NULL);
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();
}
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) {
MB_TAG_SPARSE, MF_ZERO);
CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
};
auto add_hdiv_rt_field = [&](const std::string field_name, const int order,
const int dim) {
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
};
auto add_l2_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MB_TAG_DENSE, MF_ZERO);
};
auto add_h1_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
MB_TAG_DENSE, MF_ZERO);
CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
};
auto add_bubble_field = [this, meshset](const std::string field_name,
const int order, const int dim) {
// Modify field
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;
};
// 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);
}
// set finite element fields
auto add_field_to_fe = [this](const std::string fe,
const std::string field_name) {
};
CHKERR add_field_to_fe(elementVolumeName, rotAxis);
}
// build finite elements data structures
}
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) {
};
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_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
CHKERR add_field_to_fe(naturalBcElement, piolaStress);
CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
}
// find adjacencies between finite elements and dofs
// Create coupled problem
dM = createSmartDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
BitRefLevel().set());
CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
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})
prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0, 100,
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);
// Create elastic bubble-problem
dmElasticSchurBubble = createSmartDM(mField.get_comm(), "DMMOFEM");
"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);
// Create elastic omega-problem
dmElasticSchurOmega = createSmartDM(mField.get_comm(), "DMMOFEM");
"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);
// Create elastic tet_stress-problem
dmElasticSchurSpatialDisp = createSmartDM(mField.get_comm(), "DMMOFEM");
"ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurSpatialDisp, piolaStress.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp, naturalBcElement.c_str());
CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
CHKERR DMMoFEMSetSquareProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurSpatialDisp, PETSC_TRUE);
{
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]);
}
}
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 {
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
case HDIV:
break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
}
}
private:
// 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
DataForcesAndSourcesCore &data = cTx->dAta;
// Get approximation order on element. Note that bubble functions are only
// on tetrahedron.
const int order = data.dataOnEntities[MBTET][0].getOrder();
/// number of integration points
const int nb_gauss_pts = pts.size2();
// 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);
}
};
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) {
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
// 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));
}
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(
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
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) {
}
}
}
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();
fe_rhs->getOpPtrVector().push_back(
fe_lhs->getOpPtrVector().push_back(
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
}
}
}
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);
CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, 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);
CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, 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);
}
}
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->loopsIJacobian, Suu, aoSuu);
CHKERR add_schur_streach_pre(ts_ctx->preProcessIJacobian, Suu, aoSuu);
CHKERR add_schur_streach_pre(ts_ctx->postProcessIJacobian, 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->loopsIJacobian, SBubble,
aoSBubble);
CHKERR add_schur_bubble_pre(ts_ctx->preProcessIJacobian, SBubble,
aoSBubble);
CHKERR add_schur_bubble_pre(ts_ctx->postProcessIJacobian, 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->loopsIJacobian, SOmega,
aoSOmega);
CHKERR add_schur_omega_pre(ts_ctx->preProcessIJacobian, SOmega,
aoSOmega);
CHKERR add_schur_omega_pre(ts_ctx->postProcessIJacobian, 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->loopsIJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcessIJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcessIJacobian, Sw,
aoSw);
} else
"Problem does not have sub-problem data");
} else
"Problem does not have sub-problem data");
} else
"Problem does not have sub-problem data");
} else
"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,
DataForcesAndSourcesCore::EntData &data) {
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 = determinantTensor3by3(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(), nullptr, 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->getLoopsMonitor().push_back(
CHKERR TSAppendOptionsPrefix(ts, "elastic_");
CHKERR TSSetFromOptions(ts);
}
CHKERR TSSetDM(ts, dmElastic);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
CHKERR SNESMonitorSet(
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
PetscSection section;
CHKERR 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(
&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,
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);
// CHKERR TSGetSNES(ts, &snes);
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
MOFEM_LOG("EP", Sev::inform)
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
PETSC_NULL);
if (test_cook_flg)
if (lin_solver_iterations != 2)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Expected number of iterations is different than expected");
}
const std::string file) {
if (!dataAtPts) {
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
post_proc.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
post_proc.getOpPtrVector().push_back(new OpL2Transform());
post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
piolaStress, dataAtPts->getApproxPAtPts()));
post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
post_proc.getOpPtrVector().push_back(
streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
materialGradient, dataAtPts->getBigGAtPts(), MBTET));
post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
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
static Index< 'p', 3 > p
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
Eshelbian plasticity interface.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:338
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
static PetscErrorCode ierr
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ MF_EXIST
Definition: definitions.h:100
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
static double phi
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual Field * get_field_structure(const std::string &name)=0
get field structure
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#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.
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
double f(const double v)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
BcRot(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
SmartPetscObj< DM > dmElasticSchurBubble
Sub problem of dmElastic Schur.
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
SmartPetscObj< DM > dmElasticSchurStreach
Sub problem of dmElastic Schur.
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
MoFEMErrorCode setGenericFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_lhs)
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
MoFEMErrorCode setElasticElementOps(const int tag)
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode 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)
SmartPetscObj< DM > dmElasticSchurOmega
Sub problem of dmElastic Schur.
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
EshelbianCore(MoFEM::Interface &m_field)
SmartPetscObj< DM > dM
Coupled problem all fields.
MoFEMErrorCode solveElastic(TS ts, Vec x)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe)
SmartPetscObj< DM > dmElastic
Elastic problem.
boost::shared_ptr< EpFEMethod > schurAssembly
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
boost::shared_ptr< PhysicalEquations > physicalEquations
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
MoFEMErrorCode addFields(const EntityHandle meshset=0)
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
int operator()(int p_row, int p_col, int p_data) const
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
int operator()(int p_row, int p_col, int p_data) const
Set integration rule to boundary elements.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Deprecated interface functions.
Class used to pass element data to calculate base functions on tet,triangle,edge.
Data on single entity (This is passed as argument to DataOperator::doWork)
Field evaluator interface.
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.
Provide data structure for (tensor) field approximation.
@ OPROW
operator doWork function is executed on FE rows
structure to get information form mofem into EntitiesFieldData
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
transform Hdiv base fluxes from reference element to physical triangle
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
Problem manager is used to build and partition problems.
intrusive_ptr for managing petsc objects
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.
Set integration rule.
VolEle::UserDataOperator VolOp
VolumeElementForcesAndSourcesCore VolEle