#include <BasicFiniteElements.hpp>
#include <boost/math/constants/constants.hpp>
*iface = NULL;
*iface = const_cast<EshelbianCore *>(this);
}
}
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") {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
EshelbianCore::~EshelbianCore() {}
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
"none");
spaceOrder = 1;
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
spaceOrder, &spaceOrder, PETSC_NULL);
materialOrder = 1;
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
"", materialOrder, &materialOrder, PETSC_NULL);
alphaU = 0;
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"", alphaU,
&alphaU, PETSC_NULL);
alphaW = 0;
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"", alphaW,
&alphaW, PETSC_NULL);
alphaRho = 0;
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"", alphaRho,
&alphaRho, PETSC_NULL);
precEps = 1e-3;
CHKERR PetscOptionsScalar(
"-preconditioner_eps",
"preconditioner_eps",
"",
precEps, &precEps, PETSC_NULL);
ierr = PetscOptionsEnd();
}
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 =
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,
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,
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,
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,
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,
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);
};
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);
CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
CHKERR add_l2_field(materialGradient +
"0", materialOrder - 1, 9);
}
EshelbianCore::addVolumeFiniteElement(
const EntityHandle meshset) {
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_ents_to_finite_element_by_type(meshset, MBTET,
elementVolumeName);
CHKERR add_field_to_fe(elementVolumeName, piolaStress);
CHKERR add_field_to_fe(elementVolumeName, bubbleField);
CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
CHKERR add_field_to_fe(elementVolumeName, streachTensor);
CHKERR add_field_to_fe(elementVolumeName, rotAxis);
CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
CHKERR add_field_to_fe(elementVolumeName, streachTensor);
CHKERR add_field_to_fe(elementVolumeName, materialGradient);
CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
materialGradient + "0");
}
CHKERR mField.build_finite_elements(elementVolumeName);
}
EshelbianCore::addBoundaryFiniteElement(
const EntityHandle meshset) {
auto bc_elements_add_to_range = [&](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);
}
}
};
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_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_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);
}
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
auto remove_dofs_on_essential_spatial_stress_boundary =
[&](const std::string prb_name) {
prb_name, piolaStress, (*bcSpatialFreeTraction)[
d],
d,
d,
NOISY,
true);
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary(
"ESHELBY_PLASTICITY");
dmElasticSchurStreach =
createSmartDM(mField.get_comm(),
"DMMOFEM");
"ELASTIC_PROBLEM_STREACH_SCHUR");
CHKERR DMSetUp(dmElasticSchurStreach);
dmElasticSchurBubble =
createSmartDM(mField.get_comm(),
"DMMOFEM");
"ELASTIC_PROBLEM_BUBBLE_SCHUR");
CHKERR DMSetUp(dmElasticSchurBubble);
dmElasticSchurOmega =
createSmartDM(mField.get_comm(),
"DMMOFEM");
"ELASTIC_PROBLEM_OMEGA_SCHUR");
CHKERR DMSetUp(dmElasticSchurOmega);
dmElasticSchurSpatialDisp =
createSmartDM(mField.get_comm(),
"DMMOFEM");
"ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
elementVolumeName.c_str());
essentialBcElement.c_str());
CHKERR DMSetUp(dmElasticSchurSpatialDisp);
{
PetscSection section;
§ion);
CHKERR DMSetDefaultSection(dmElastic, section);
CHKERR DMSetDefaultGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(§ion);
}
}
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 =
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) {
"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);
}
}
}
int operator()(int p_row, int p_col, int p_data) const {
return 2 * (p_data + 1);
}
};
int operator()(int p_row, int p_col, int p_data) const {
return 2 * (p_data);
}
};
CGGUserPolynomialBase() {}
~CGGUserPolynomialBase() {}
*iface = NULL;
*iface = const_cast<CGGUserPolynomialBase *>(this);
} else {
}
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
CHKERR getValueHdivForCGGBubble(pts);
break;
default:
}
}
private:
"Wrong base, should be USER_BASE");
}
const int nb_gauss_pts = pts.size2();
shapeFun.resize(nb_gauss_pts, 4, false);
&pts(2, 0), nb_gauss_pts);
double diff_shape_fun[12];
phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
&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));
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());
if (!dataAtPts) {
dataAtPts =
dataAtPts->physicsPtr = physicalEquations;
}
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
piolaStress, dataAtPts->getDivPAtPts()));
eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
eshelbyStress, dataAtPts->getDivSigmaAtPts()));
streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
materialGradient, dataAtPts->getBigGAtPts(), MBTET));
materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
fe->getOpPtrVector().push_back(
streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
}
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
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) {
CHKERR setBaseVolumeElementOps(tag,
true,
false, fe_rhs);
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));
}
CHKERR setBaseVolumeElementOps(tag,
true,
true, fe_lhs);
if (add_elastic) {
fe_lhs->getOpPtrVector().push_back(
new OpSpatialSchurBegin(spatialDisp, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
streachTensor, streachTensor, dataAtPts, alphaU));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
streachTensor, bubbleField, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
spatialDisp, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
spatialDisp, spatialDisp, dataAtPts, alphaW, alphaRho));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
bubbleField, rotAxis, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
streachTensor, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
streachTensor, rotAxis, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
rotAxis, piolaStress, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
rotAxis, bubbleField, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
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);
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr));
fe_rhs->getOpPtrVector().push_back(
new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
}
}
CHKERR setGenericVolumeElementOps(tag,
true,
false, elasticFeRhs,
elasticFeLhs);
CHKERR setGenericFaceElementOps(
true,
false, elasticBcRhs, elasticBcLhs);
}
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()) {
null);
null);
null);
spatial_traction_bc);
null);
null);
} else {
null);
null);
null);
spatial_traction_bc);
null);
null);
}
}
boost::shared_ptr<TsCtx> 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);
}
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addStreachSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addStreachSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addBubbleSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addBubbleSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
fe_cast->addOmegaSchurMatrix(S, aoS);
else
};
for (auto &fe : list)
if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
fe_cast->addOmegaSchurMatrix(S, ao);
else
};
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 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 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);
&schur_spatial_disp_prb_ptr);
if (auto spatial_disp_data =
CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
aoSw);
CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
aoSw);
} else
"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");
EshelbianCore &eP;
boost::shared_ptr<SetPtsData> dataFieldEval;
boost::shared_ptr<VolEle> volPostProcEnergy;
boost::shared_ptr<double> gEnergy;
: eP(ep),
volPostProcEnergy(
new VolEle(ep.mField)), gEnergy(
new double) {
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));
eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
vol_ele->getOpPtrVector().push_back(
eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
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();
}
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;
PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
FILE_MODE_WRITE, &viewer);
CHKERR PetscViewerDestroy(&viewer);
CHKERR eP.postProcessResults(1,
"out_sol_elastic_" + get_step(ts_step) +
".h5m");
*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)
str(str) {}
if (getGaussPts().size2()) {
auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
auto t_approx_P =
getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
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;
<< add() << "comm rank " << eP.mField.get_comm_rank();
<< add() << "coords at gauss pts " << getCoordsAtGaussPts();
<< add() << "w " << eP.dataAtPts->wAtPts;
<< add() << "Piola " << eP.dataAtPts->approxPAtPts;
<< add() << "Cauchy " << print_tensor(t_cauchy);
}
}
}
};
if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
fe_ptr->getOpPtrVector().push_back(
new OpPrint(eP,
point, str));
->evalFEAtThePoint3D(
point.data(), 1e-12, problemPtr->getName(),
"EP", dataFieldEval,
eP.mField.get_comm_rank(),
fe_ptr->getOpPtrVector().pop_back();
}
};
std::array<double, 3> pointA = {48., 60., 5.};
CHKERR post_proc_at_points(pointA,
"Point A");
std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
CHKERR post_proc_at_points(pointB,
"Point B");
std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
CHKERR post_proc_at_points(pointC,
"Point C");
}
};
boost::shared_ptr<FEMethod> monitor_ptr(
new Monitor(*
this));
ts_ctx->get_loops_to_do_Monitor().push_back(
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
}
CHKERR TSSetDM(ts, dmElastic);
SNES snes;
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
PetscSection section;
CHKERR DMGetDefaultSection(dmElastic, §ion);
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 VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_uu_field_split;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
if (is_uu_field_split) {
map<std::string, IS> is_map;
for (int ff = 0; ff != num_fields; ff++) {
const char *field_name;
CHKERR PetscSectionGetFieldName(section, ff, &field_name);
&is_map[field_name]);
}
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);
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) {
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);
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) {
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);
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) {
&schur_w_ptr);
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);
}
}
}
}
}
}
}
}
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
CHKERR TSSolve(ts, PETSC_NULL);
}
const std::string file) {
if (!dataAtPts) {
dataAtPts =
}
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
materialGradient, dataAtPts->getBigGAtPts(), MBTET));
spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
post_proc.
getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
spatialDisp, tag, true, false, dataAtPts, physicalEquations));
}
}