Eshelbian plasticity implementation.
#include <boost/math/constants/constants.hpp>
#ifdef PYTHON_SFD
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
#pragma message "With PYTHON_SFD"
#else
#pragma message "Without PYTHON_SFD"
#endif
*iface = const_cast<EshelbianCore *>(this);
return 0;
}
}
: mField(m_field), piolaStress("P"), eshelbyStress("S"),
spatialL2Disp("wL2"), spatialH1Disp("wH1"), materialL2Disp("W"),
streachTensor("u"), rotAxis("omega"), materialGradient("G"),
tauField("TAU"), lambdaField("LAMBDA"), bubbleField("BUBBLE"),
elementVolumeName("EP"), naturalBcElement("NATURAL_BC"),
essentialBcElement("ESSENTIAL_BC"), skinElement("SKIN_ELEMENT"),
contactElement("CONTACT_ELEMENT") {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
const char *list_rots[] = {"small", "moderate", "large"};
const char *list_stretches[] = {"linear", "log"};
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"",
alphaU,
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"",
alphaW,
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"",
alphaRho,
CHKERR PetscOptionsScalar(
"-preconditioner_eps",
"preconditioner_eps",
"",
CHKERR PetscOptionsScalar(
"-preconditioner_eps_omega",
"preconditioner_eps",
CHKERR PetscOptionsScalar(
"-preconditioner_eps_w",
"preconditioner_eps",
"",
CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
PETSC_NULL);
list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
ierr = PetscOptionsEnd();
switch (choice_stretch) {
break;
break;
default:
break;
};
else
MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
MOFEM_LOG(
"EP", Sev::inform) <<
"Stretch " << list_stretches[choice_stretch];
}
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
tets_skin = subtract(tets_skin,
v.faces);
}
tets_skin = subtract(tets_skin,
v.faces);
}
tets_skin = subtract(tets_skin,
v.faces);
}
auto subtract_faces_where_displacements_are_applied =
[&](const std::string block_name) {
auto reg_exp = std::regex((boost::format("%s(.*)") % block_name).str());
for (
auto m : msh_mng->getCubitMeshsetPtr(reg_exp)) {
tets_skin = subtract(tets_skin, faces);
<<
"Subtracting " <<
m->getName() <<
" " << faces.size();
}
};
CHKERR subtract_faces_where_displacements_are_applied(
"CONTACT");
moab::Interface::UNION);
Range faces_not_on_the_skin = subtract(faces, tets_skin);
};
auto add_l2_field = [
this, meshset](
const std::string
field_name,
};
auto add_h1_field = [
this, meshset](
const std::string
field_name,
};
auto add_bubble_field = [
this, meshset](
const std::string
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;
};
}
auto add_field_to_fe = [this](const std::string fe,
};
}
}
auto bc_elements_add_to_range = [&](
const std::string block_name,
Range &r) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
true);
r.merge(faces);
}
};
auto add_field_to_fe = [this](const std::string fe,
};
Range natural_bc_elements;
natural_bc_elements.merge(
v.faces);
}
}
natural_bc_elements.merge(
v.faces);
}
}
Range essentail_bc_elements;
essentail_bc_elements.merge(
v.faces);
}
}
auto get_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto &&skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto skin = filter_true_skin(get_skin());
CHKERR bc_elements_add_to_range(
"CONTACT", contact_range);
MOFEM_LOG(
"EP", Sev::inform) <<
"Contact elements " << contact_range.size();
}
CHKERR DMMoFEMSetDestroyProblem(
dM, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(
dM, PETSC_TRUE);
auto remove_dofs_on_essential_spatial_stress_boundary =
[&](const std::string prb_name) {
for (int d : {0, 1, 2})
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary(
"ESHELBY_PLASTICITY");
{
PetscSection section;
§ion);
CHKERR PetscSectionDestroy(§ion);
}
}
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
<< "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " << faces.size();
}
: blockName(name), faces(faces) {
vals.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
}
theta = attr[3];
}
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
<< "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " << faces.size();
}
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
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;
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
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 p_data + p_data + (p_data - 1);
}
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data; }
};
return 0;
}
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");
}
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);
nb_gauss_pts);
}
};
const int tag, const bool do_rhs, const bool do_lhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
}
fe->getOpPtrVector().push_back(
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
}
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(
dataAtPts));
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto time_scale = boost::make_shared<TimeScale>();
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<GAUSS>;
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
fe_rhs->getOpPtrVector(),
mField,
"w", {time_scale},
"BODY_FORCE",
Sev::inform);
}
if (add_elastic) {
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(
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(
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
if constexpr (
A == AssemblyType::SCHUR) {
fe_lhs->getOpPtrVector().push_back(
}));
std::numeric_limits<double>::epsilon()) {
fe_lhs->getOpPtrVector().push_back(new OpMassStab(
[
this](
double,
double,
double) {
return precEpsW; }));
}
}
if (add_material) {
}
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
{HDIV});
{HDIV});
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
fe_rhs->getOpPtrVector().push_back(
}
}
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<BoundaryEle>(
mField);
fe_lhs = boost::make_shared<BoundaryEle>(
mField);
{HDIV});
{HDIV});
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_side = [&](auto disp_ptr, auto col_ent_data_ptr) {
op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
op_loop_side->getOpPtrVector().push_back(
if (col_ent_data_ptr != nullptr) {
op_loop_side->getOpPtrVector().push_back(
}
return op_loop_side;
};
auto add_ops_lhs = [&](auto &pip) {
auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(
get_op_side(common_data_ptr->contactDispPtr(), col_ent_data_ptr));
};
auto add_ops_rhs = [&](auto &pip) {
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(get_op_side(common_data_ptr->contactDispPtr(), nullptr));
};
CHKERR add_ops_lhs(fe_lhs->getOpPtrVector());
CHKERR add_ops_rhs(fe_rhs->getOpPtrVector());
}
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_contact_bc = [&]() {
return op_loop_side;
};
auto op_contact_bc = get_op_contact_bc();
}
boost::shared_ptr<FEMethod> null;
boost::shared_ptr<FeTractionBc> spatial_traction_bc(
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);
}
}
#ifdef PYTHON_SFD
boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
if (file_exists("sdf.py")) {
MOFEM_LOG(
"EP", Sev::inform) <<
"sdf.py file found";
sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdf_python_ptr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
} else {
MOFEM_LOG(
"EP", Sev::warning) <<
"sdf.py file NOT found";
}
#else
#endif
boost::shared_ptr<TsCtx>
ts_ctx;
CHKERR TSMonitorSet(ts, TsMonitorSet,
ts_ctx.get(), PETSC_NULL);
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(), {HDIV, H1, L2});
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));
eP.spatialL2Disp, eP.dataAtPts->getSmallWAtPts(), MBTET));
vol_ele->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
};
auto set_element_for_post_process = [&]() {
volPostProcEnergy->getRuleHook =
VolRule();
volPostProcEnergy->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
volPostProcEnergy->getOpPtrVector(), {HDIV, H1, L2});
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.spatialL2Disp, eP.dataAtPts->getSmallWAtPts(), MBTET));
volPostProcEnergy->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
volPostProcEnergy->getOpPtrVector().push_back(
new OpCalculateStrainEnergy(eP.spatialL2Disp, 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 = 0;
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 (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;
return s.str();
};
std::ostringstream print;
<< add() << "comm rank " << eP.mField.get_comm_rank();
<< add() << "point " << getVectorAdaptor(point.data(), 3);
<< 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., 0.};
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));
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
SNES snes;
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
PetscSection section;
int num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
}
CHKERR DMoFEMMeshToLocalVector(
dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
boost::shared_ptr<SetUpSchur> schur_ptr;
if constexpr (
A == AssemblyType::SCHUR) {
auto p = matDuplicate(
m, MAT_DO_NOT_COPY_VALUES);
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR TSSetI2Function(ts,
f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetI2Jacobian(ts,
m,
p, PETSC_NULL, PETSC_NULL);
} else {
CHKERR TSSetIFunction(ts,
f, PETSC_NULL, PETSC_NULL);
CHKERR TSSetIJacobian(ts,
m,
p, PETSC_NULL, PETSC_NULL);
}
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
};
};
};
};
}
if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon()) {
Vec xx;
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
auto create_post_step_ksp = [this]() {
auto set_up = [&]() {
GAUSS>::OpBaseTimesVector<1, 3, 1>;
auto fe_lhs = boost::make_shared<DomainEle>(
mField);
auto fe_rhs = boost::make_shared<DomainEle>(
mField);
fe_lhs->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
fe_rhs->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
fe_lhs->getOpPtrVector().push_back(
auto w_ptr = boost::make_shared<MatrixDouble>();
fe_rhs->getOpPtrVector().push_back(
fe_lhs, nullptr, nullptr);
nullptr, nullptr);
CHKERR KSPAppendOptionsPrefix(ksp,
"prjspatial_");
CHKERR KSPSetFromOptions(ksp);
};
return ksp;
};
prjKsp = create_post_step_ksp();
auto post_step_fun = [](TS ts) {
CHKERR VecGhostUpdateBegin(
prjD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
prjD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(
prjDM,
prjD, INSERT_VALUES, SCATTER_REVERSE);
};
CHKERR TSSetPostStep(ts, post_step_fun);
CHKERR TSSolve(ts, PETSC_NULL);
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL,
"",
"-test_cook", &test_cook_flg,
PETSC_NULL);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
}
}
const std::string file) {
}
auto domain_ops = [&](auto &fe) {
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
fe.getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(
dataAtPts));
};
CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
post_proc.getOpPtrVector().push_back(op_loop_side);
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
spatialL2Disp,
CHKERR post_proc.writeFile(file.c_str());
}
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(
mField);
auto post_proc_norm_rule_hook = [](int, int,
int p) ->
int {
};
post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
post_proc_norm_fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
auto norms_vec =
CHKERR VecZeroEntries(norms_vec);
auto u_l2_ptr = boost::make_shared<MatrixDouble>();
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
u_h1_ptr));
auto piola_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
*post_proc_norm_fe);
CHKERR VecAssemblyBegin(norms_vec);
CHKERR VecAssemblyEnd(norms_vec);
const double *norms;
CHKERR VecGetArrayRead(norms_vec, &norms);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
<< "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
<< "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
}
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
}
private:
};
auto get_ents_by_dim = [&](
int dim) {
return ents;
};
std::vector<std::string> schur_field_list{
CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(
SPACE_DIM - 1));
std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr};
if (field_list.size() != dm_range_list.size()) {
"Both ranges should have the same size");
}
int r_idx = 0;
for (
auto f : field_list) {
MOFEM_LOG(
"EP", Sev::inform) <<
"Add schur field: " <<
f;
CHKERR DMMoFEMAddSubFieldRow(dm_sub,
f.c_str(), dm_range_list[r_idx]);
CHKERR DMMoFEMAddSubFieldCol(dm_sub,
f.c_str(), dm_range_list[r_idx]);
++r_idx;
}
};
std::vector<SmartPetscObj<IS>> is_vec;
std::vector<Range *> range_list_ptr(schur_field_list.size(), nullptr);
range_list_ptr.back() = &vols;
int r_idx = 0;
for (
auto f : schur_field_list) {
range_list_ptr[r_idx]);
is_vec.push_back(is);
++r_idx;
}
if (is_vec.size()) {
is_a00 = is_vec[0];
for (
auto i = 1;
i < is_vec.size(); ++
i) {
IS is_union_raw;
CHKERR ISExpand(is_a00, is_vec[
i], &is_union_raw);
}
}
};
PC pc;
CHKERR KSPSetFromOptions(solver);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
CHKERR create_schur_dm(schur_dm);
"Is expected that schur matrix is not allocated. This is "
"possible only is PC is set up twice");
}
S = createDMMatrix(schur_dm);
auto set_ops = [&]() {
std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
nullptr);
ranges_list.back() =
boost::make_shared<Range>(get_ents_by_dim(
SPACE_DIM));
std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
std::vector<bool> symm_list(schur_field_list.size(), false);
auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
aoUp = createAOMappingIS(dm_is, PETSC_NULL);
ao_list, mat_list, symm_list,
false));
ao_list, mat_list, symm_list,
false));
};
auto set_pc = [&]() {
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
};
} else {
}
}
MOFEM_LOG(
"EP", Sev::verbose) <<
"Zero Schur";
}
}
MOFEM_LOG(
"EP", Sev::verbose) <<
"Assemble Schur";
CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
}
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyBegin(
P, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
P, MAT_FINAL_ASSEMBLY);
CHKERR MatAXPY(
P, 1,
M, SAME_NONZERO_PATTERN);
}
boost::shared_ptr<EshelbianCore::SetUpSchur>
EshelbianCore *ep_core_ptr) {
return boost::shared_ptr<SetUpSchur>(
}
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto disp_bc = bc.second->dispBcPtr) {
std::vector<double> block_attributes(6, 0.);
if (disp_bc->data.flag1 == 1) {
block_attributes[0] = disp_bc->data.value1;
block_attributes[3] = 1;
}
if (disp_bc->data.flag2 == 1) {
block_attributes[1] = disp_bc->data.value2;
block_attributes[4] = 1;
}
if (disp_bc->data.flag3 == 1) {
block_attributes[2] = disp_bc->data.value3;
block_attributes[5] = 1;
}
auto faces = bc.second->bcEnts.subset_by_dimension(2);
}
}
}
false, false);
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto force_bc = bc.second->forceBcPtr) {
std::vector<double> block_attributes(6, 0.);
block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
block_attributes[3] = 1;
block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
block_attributes[4] = 1;
block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
block_attributes[5] = 1;
auto faces = bc.second->bcEnts.subset_by_dimension(2);
}
}
}
}
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Eshelbian plasticity interface.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
static PetscErrorCode ierr
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FieldEvaluatorInterface::SetPtsData SetPtsData
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
FTensor::Index< 'm', SPACE_DIM > m
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
static double exponentBase
double dd_f_linear(const double v)
double f_log(const double v)
SmartPetscObj< DM > prjDM
static boost::function< double(const double)> d_f
double d_f_linear(const double v)
static boost::function< double(const double)> f
double dd_f_log(const double v)
double f_linear(const double v)
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.
static boost::function< double(const double)> dd_f
SmartPetscObj< Vec > prjD
Unknown vector for projection spatial displacement.
SmartPetscObj< Vec > prjF
double d_f_log(const double v)
static enum RotSelector rotSelector
SmartPetscObj< KSP > prjKsp
KSP for projection spatial displacement on H1 space.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
constexpr double t
plate stiffness
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)
EntPolynomialBaseCtx * cTx
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
MoFEMErrorCode getSpatialDispBc()
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
const std::string piolaStress
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
const std::string essentialBcElement
const std::string bubbleField
MoFEMErrorCode solveElastic(TS ts, Mat m, Vec f, Vec x)
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe)
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
MoFEMErrorCode setContactElementOps(boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
EshelbianCore(MoFEM::Interface &m_field)
SmartPetscObj< DM > dM
Coupled problem all fields.
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
SmartPetscObj< DM > dmElastic
Elastic problem.
const std::string skinElement
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
boost::shared_ptr< PhysicalEquations > physicalEquations
const std::string spatialL2Disp
const std::string spatialH1Disp
MoFEMErrorCode setVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
boost::shared_ptr< FaceElementForcesAndSourcesCore > contactRhs
MoFEMErrorCode gettingNorms()
[Getting norms]
MoFEMErrorCode getSpatialTractionBc()
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
MoFEMErrorCode addFields(const EntityHandle meshset=0)
const std::string naturalBcElement
const std::string streachTensor
MoFEMErrorCode getOptions()
const std::string rotAxis
MoFEMErrorCode setElasticElementToTs(DM dm)
const std::string elementVolumeName
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
Remove all, but entities where kinematic constrains are applied.
const std::string contactElement
MoFEM::Interface & mField
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)
SmartPetscObj< IS > isA00
MoFEMErrorCode setUp(KSP solver)
MoFEM::Interface & mField
EshelbianCore * epCorePtr
virtual ~SetUpSchurImpl()
SmartPetscObj< IS > isShift
MoFEMErrorCode postProc()
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.
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
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.
Definition of the displacement bc data structure.
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)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
default operator for TRI element
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.
Definition of the force bc data structure.
@ 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.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.
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.
Calculate trace of vector (Hdiv/Hcurl) space.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Clear Schur complement internal data.
Assemble Schur complement.
Problem manager is used to build and partition problems.
intrusive_ptr for managing petsc objects
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
Base volume element used to integrate on skeleton.
[Push operators to pipeline]
VolEle::UserDataOperator VolOp
VolumeElementForcesAndSourcesCore VolEle