v0.14.0
dynamic_first_order_con_law.cpp

Explicit first order conservation laws for solid dynamics

/**
* \file dynamic_first_order_con_law.cpp
* \example dynamic_first_order_con_law.cpp
*
* Explicit first order conservation laws for solid dynamics
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
template <typename T> inline double trace(FTensor::Tensor2<T, 2, 2> &t_stress) {
constexpr double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1));
};
template <typename T> inline double trace(FTensor::Tensor2<T, 3, 3> &t_stress) {
constexpr double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2));
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using BoundaryEle =
template <int DIM> struct PostProcEleByDim;
template <> struct PostProcEleByDim<2> {
};
template <> struct PostProcEleByDim<3> {
};
PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
IntegrationType::GAUSS>::OpBaseTimesVector<1, SPACE_DIM * SPACE_DIM, 1>;
/** \brief Save field DOFS on vertices/tags
*/
constexpr double omega = 1.;
constexpr double young_modulus = 1.;
constexpr double poisson_ratio = 0.;
double bulk_modulus_K = young_modulus / (3. * (1. - 2. * poisson_ratio));
double shear_modulus_G = young_modulus / (2. * (1. + poisson_ratio));
double mu = young_modulus / (2. * (1. + poisson_ratio));
((1. + poisson_ratio) * (1. - 2. * poisson_ratio));
// Operator to Calculate F
template <int DIM_0, int DIM_1>
OpCalculateFStab(boost::shared_ptr<MatrixDouble> def_grad_ptr,
boost::shared_ptr<MatrixDouble> def_grad_stab_ptr,
boost::shared_ptr<MatrixDouble> def_grad_dot_ptr,
double tau_F_ptr, double xi_F_ptr,
boost::shared_ptr<MatrixDouble> grad_x_ptr,
boost::shared_ptr<MatrixDouble> grad_vel_ptr)
defGradPtr(def_grad_ptr), defGradStabPtr(def_grad_stab_ptr),
defGradDotPtr(def_grad_dot_ptr), tauFPtr(tau_F_ptr), xiF(xi_F_ptr),
gradxPtr(grad_x_ptr), gradVelPtr(grad_vel_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
// Define Indicies
// Number of Gauss points
const size_t nb_gauss_pts = getGaussPts().size2();
defGradStabPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
defGradStabPtr->clear();
// Extract matrix from data matrix
auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
auto t_Fstab = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradStabPtr);
auto t_F_dot = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradDotPtr);
// tau_F = alpha deltaT
auto tau_F = tauFPtr;
double xi_F = xiF;
auto t_gradx = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradxPtr);
auto t_gradVel = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradVelPtr);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
// Stabilised Deformation Gradient
t_Fstab(i, j) = t_F(i, j) + tau_F * (t_gradVel(i, j) - t_F_dot(i, j)) +
xi_F * (t_gradx(i, j) - t_F(i, j));
++t_F;
++t_Fstab;
++t_gradVel;
++t_F_dot;
++t_gradx;
}
}
private:
double tauFPtr;
double xiF;
boost::shared_ptr<MatrixDouble> defGradPtr;
boost::shared_ptr<MatrixDouble> defGradStabPtr;
boost::shared_ptr<MatrixDouble> defGradDotPtr;
boost::shared_ptr<MatrixDouble> gradxPtr;
boost::shared_ptr<MatrixDouble> gradVelPtr;
};
// Operator to Calculate P
template <int DIM_0, int DIM_1>
OpCalculatePiola(double shear_modulus, double bulk_modulus, double m_u,
double lambda_lamme,
boost::shared_ptr<MatrixDouble> first_piola_ptr,
boost::shared_ptr<MatrixDouble> def_grad_ptr)
shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
defGradPtr(def_grad_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
// Define Indicies
// Define Kronecker Delta
// Number of Gauss points
const size_t nb_gauss_pts = getGaussPts().size2();
// Resize Piola
firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check
firstPiolaPtr->clear();
// Extract matrix from data matrix
auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
const double two_o_three = 2. / 3.;
const double trace_t_dk = DIM_0;
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_P(i, j) = shearModulus * (t_F(i, j) + t_F(j, i) - 2. * t_kd(i, j) -
two_o_three * trace(t_F) * t_kd(i, j) +
two_o_three * trace_t_dk * t_kd(i, j)) +
bulkModulus * trace(t_F) * t_kd(i, j) -
bulkModulus * trace_t_dk * t_kd(i, j);
++t_F;
++t_P;
}
}
private:
double shearModulus;
double bulkModulus;
double mU;
double lammeLambda;
boost::shared_ptr<MatrixDouble> firstPiolaPtr;
boost::shared_ptr<MatrixDouble> defGradPtr;
};
template <int DIM>
OpCalculateDisplacement(boost::shared_ptr<MatrixDouble> spatial_pos_ptr,
boost::shared_ptr<MatrixDouble> reference_pos_ptr,
boost::shared_ptr<MatrixDouble> u_ptr)
xPtr(spatial_pos_ptr), XPtr(reference_pos_ptr), uPtr(u_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
// Define Indicies
// Number of Gauss points
const size_t nb_gauss_pts = getGaussPts().size2();
uPtr->resize(DIM, nb_gauss_pts, false); // ignatios check
uPtr->clear();
// Extract matrix from data matrix
auto t_x = getFTensor1FromMat<DIM>(*xPtr);
auto t_X = getFTensor1FromMat<DIM>(*XPtr);
auto t_u = getFTensor1FromMat<DIM>(*uPtr);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_u(i) = t_x(i) - t_X(i);
++t_u;
++t_x;
++t_X;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> XPtr;
boost::shared_ptr<MatrixDouble> uPtr;
};
template <int DIM_0, int DIM_1>
double shear_modulus, double bulk_modulus, double m_u,
double lambda_lamme, boost::shared_ptr<MatrixDouble> first_piola_ptr,
boost::shared_ptr<MatrixDouble> def_grad_ptr,
boost::shared_ptr<MatrixDouble> inv_def_grad_ptr,
boost::shared_ptr<VectorDouble> det)
shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
defGradPtr(def_grad_ptr), invDefGradPtr(inv_def_grad_ptr), dEt(det) {}
MoFEMErrorCode doWork(int side, EntityType type,
// Define Indicies
// Define Kronecker Delta
// Number of Gauss points
const size_t nb_gauss_pts = getGaussPts().size2();
// Resize Piola
firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check
firstPiolaPtr->clear();
// Extract matrix from data matrix
auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
auto t_inv_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*invDefGradPtr);
auto t_det = getFTensor0FromVec<1>(*dEt);
const double two_o_three = 2. / 3.;
const double one_o_three = 1. / 3.;
const double bulk_mod = bulkModulus;
const double shear_mod = shearModulus;
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
// Nearly incompressible NH
// volumetric part
t_P(i, j) = bulk_mod * (t_det - 1.) * t_det * t_inv_F(j, i);
// deviatoric part
t_P(i, j) +=
shear_mod * pow(t_det, two_o_three) *
(t_F(i, j) - one_o_three * (t_F(l, k) * t_F(l, k)) * t_inv_F(j, i));
++t_F;
++t_P;
++t_inv_F;
++t_det;
}
}
private:
double shearModulus;
double bulkModulus;
double mU;
double lammeLambda;
boost::shared_ptr<MatrixDouble> firstPiolaPtr;
boost::shared_ptr<MatrixDouble> defGradPtr;
boost::shared_ptr<MatrixDouble> invDefGradPtr;
boost::shared_ptr<VectorDouble> dEt;
};
template <int DIM_0, int DIM_1>
boost::shared_ptr<MatrixDouble> def_grad_ptr,
boost::shared_ptr<MatrixDouble> grad_tensor_ptr)
defGradPtr(def_grad_ptr), gradTensorPtr(grad_tensor_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
// Define Indicies
// Define Kronecker Delta
// Number of Gauss points
const size_t nb_gauss_pts = getGaussPts().size2();
// Resize Piola
defGradPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
defGradPtr->clear();
// Extract matrix from data matrix
auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
auto t_H = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradTensorPtr);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_F(i, j) = t_H(i, j) + t_kd(i, j);
++t_F;
++t_H;
}
}
private:
boost::shared_ptr<MatrixDouble> gradTensorPtr;
boost::shared_ptr<MatrixDouble> defGradPtr;
};
struct Example;
struct TSPrePostProc {
TSPrePostProc() = default;
virtual ~TSPrePostProc() = default;
/**
* @brief Used to setup TS solver
*
* @param ts
* @return MoFEMErrorCode
*/
MoFEMErrorCode tsSetUp(TS ts);
// SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
Example *fsRawPtr;
static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime,
PetscInt stageindex, Vec *Y);
static MoFEMErrorCode tsPostStep(TS ts);
static MoFEMErrorCode tsPreStep(TS ts);
};
static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
double getScale(const double time) {
return sin(2. * M_PI * MoFEM::TimeScale::getScale(time));
};
};
struct CommonData {
SmartPetscObj<Mat> M; ///< Mass matrix
SmartPetscObj<KSP> ksp; ///< Linear solver
};
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
friend struct TSPrePostProc;
struct DynamicFirstOrderConsSinusTimeScale : public MoFEM::TimeScale {
double getScale(const double time) { return 0.001 * sin(0.1 * time); };
};
struct DynamicFirstOrderConsConstantTimeScale : public MoFEM::TimeScale {
double getScale(const double time) { return 0.001; };
};
};
//! [Run problem]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
}
//! [Run problem]
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
MOFEM_LOG("WORLD", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
MOFEM_LOG("WORLD", Sev::inform)
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
base = LASTBASE;
break;
}
// Add field
CHKERR simple->addDomainField("V", H1, base, SPACE_DIM);
CHKERR simple->addBoundaryField("V", H1, base, SPACE_DIM);
CHKERR simple->addDomainField("F", H1, base, SPACE_DIM * SPACE_DIM);
CHKERR simple->addDataField("x_1", H1, base, SPACE_DIM);
CHKERR simple->addDataField("x_2", H1, base, SPACE_DIM);
CHKERR simple->addDataField("F_0", H1, base, SPACE_DIM * SPACE_DIM);
CHKERR simple->addDataField("F_dot", H1, base, SPACE_DIM * SPACE_DIM);
CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
int order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("V", order);
CHKERR simple->setFieldOrder("F", order);
CHKERR simple->setFieldOrder("F_0", order);
CHKERR simple->setFieldOrder("F_dot", order);
CHKERR simple->setFieldOrder("x_1", order);
CHKERR simple->setFieldOrder("x_2", order);
CHKERR simple->setFieldOrder("GEOMETRY", order);
CHKERR simple->setUp();
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method_x(mField, "x_1");
CHKERR mField.loop_dofs("x_1", ent_method_x);
Projection10NodeCoordsOnField ent_method_x_2(mField, "x_2");
CHKERR mField.loop_dofs("x_2", ent_method_x_2);
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
}
//! [Set up problem]
//! [Boundary condition]
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto time_scale = boost::make_shared<TimeScale>();
PetscBool sin_time_function = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-sin_time_function",
&sin_time_function, PETSC_NULL);
if (sin_time_function)
time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();
else
time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();
pipeline_mng->getBoundaryExplicitRhsFE().reset();
pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE}, "GEOMETRY");
pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField, "V",
{time_scale}, "FORCE", Sev::inform);
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setBoundaryExplicitRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "V");
auto get_pre_proc_hook = [&]() {
mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
};
pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
}
//! [Boundary condition]
MoFEMErrorCode TSPrePostProc::tsPostStage(TS ts, PetscReal stagetime,
PetscInt stageindex, Vec *Y) {
// cerr << "tsPostStage " <<"\n";
if (auto ptr = tsPrePostProc.lock()) {
auto &m_field = ptr->fsRawPtr->mField;
auto *simple = m_field.getInterface<Simple>();
auto *pipeline_mng = m_field.getInterface<PipelineManager>();
auto fb = m_field.getInterface<FieldBlas>();
double dt;
CHKERR TSGetTimeStep(ts, &dt);
double time;
CHKERR TSGetTime(ts, &time);
PetscInt num_stages;
Vec *stage_solutions;
CHKERR TSGetStages(ts, &num_stages, &stage_solutions);
PetscPrintf(PETSC_COMM_WORLD, "Check timestep %d time %e dt %e\n",
num_stages, time, dt);
const double inv_num_step = (double)num_stages;
CHKERR fb->fieldCopy(1., "x_1", "x_2");
CHKERR fb->fieldAxpy(dt, "V", "x_2");
CHKERR fb->fieldCopy(1., "x_2", "x_1");
CHKERR fb->fieldCopy(-inv_num_step / dt, "F_0", "F_dot");
CHKERR fb->fieldAxpy(inv_num_step / dt, "F", "F_dot");
CHKERR fb->fieldCopy(1., "F", "F_0");
}
}
if (auto ptr = tsPrePostProc.lock()) {
auto &m_field = ptr->fsRawPtr->mField;
auto fb = m_field.getInterface<FieldBlas>();
double dt;
CHKERR TSGetTimeStep(ts, &dt);
double time;
CHKERR TSGetTime(ts, &time);
}
}
if (auto ptr = tsPrePostProc.lock()) {
auto &m_field = ptr->fsRawPtr->mField;
auto *simple = m_field.getInterface<Simple>();
auto *pipeline_mng = m_field.getInterface<PipelineManager>();
double dt;
CHKERR TSGetTimeStep(ts, &dt);
double time;
CHKERR TSGetTime(ts, &time);
int step_num;
CHKERR TSGetStepNumber(ts, &step_num);
}
}
//! [Push operators to pipeline]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto get_body_force = [this](const double, const double, const double) {
t_source(i) = 0.;
t_source(0) = 0.1;
t_source(1) = 1.;
return t_source;
};
// specific time scaling
auto get_time_scale = [this](const double time) {
return sin(time * omega * M_PI);
};
auto apply_rhs = [&](auto &pip) {
"GEOMETRY");
// Calculate Gradient of velocity
auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
"V", mat_v_grad_ptr));
auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
gravity_vector_ptr->resize(SPACE_DIM, 1);
auto set_body_force = [&]() {
auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
double unit_weight = 0.;
CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-unit_weight", &unit_weight,
PETSC_NULL);
t_force(i) = 0;
if (SPACE_DIM == 2) {
t_force(1) = -unit_weight;
} else if (SPACE_DIM == 3) {
t_force(2) = unit_weight;
}
};
CHKERR set_body_force();
pip.push_back(new OpBodyForce("V", gravity_vector_ptr,
[](double, double, double) { return 1.; }));
// Calculate unknown F
auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
"F", mat_H_tensor_ptr));
// // Calculate F
double tau = 0.2;
CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-tau", &tau, PETSC_NULL);
double xi = 0.;
CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-xi", &xi, PETSC_NULL);
// Calculate P stab
auto one = [&](const double, const double, const double) {
return 3. * bulk_modulus_K;
};
auto minus_one = [](const double, const double, const double) {
return -1.;
};
auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
"F_dot", mat_dot_F_tensor_ptr));
// Calculate Gradient of Spatial Positions
auto mat_x_grad_ptr = boost::make_shared<MatrixDouble>();
"x_2", mat_x_grad_ptr));
auto mat_F_tensor_ptr = boost::make_shared<MatrixDouble>();
mat_F_tensor_ptr, mat_H_tensor_ptr));
auto mat_F_stab_ptr = boost::make_shared<MatrixDouble>();
mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,
mat_x_grad_ptr, mat_v_grad_ptr));
PetscBool is_linear_elasticity = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",
&is_linear_elasticity, PETSC_NULL);
auto mat_P_stab_ptr = boost::make_shared<MatrixDouble>();
if (is_linear_elasticity) {
mat_F_stab_ptr));
} else {
auto inv_F = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pip.push_back(
new OpInvertMatrix<SPACE_DIM>(mat_F_stab_ptr, det_ptr, inv_F));
// OpCalculatePiolaIncompressibleNH
mat_F_stab_ptr, inv_F, det_ptr));
}
pip.push_back(new OpGradTimesTensor2("V", mat_P_stab_ptr, minus_one));
pip.push_back(new OpRhsTestPiola("F", mat_v_grad_ptr, one));
};
CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
}
//! [Push operators to pipeline]
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
boost::shared_ptr<PostProcEle> post_proc,
boost::shared_ptr<PostProcFaceEle> post_proc_bdry,
boost::shared_ptr<MatrixDouble> velocity_field_ptr,
boost::shared_ptr<MatrixDouble> x2_field_ptr,
boost::shared_ptr<MatrixDouble> geometry_field_ptr,
std::array<double, 3> pass_field_eval_coords,
boost::shared_ptr<SetPtsData> pass_field_eval_data)
: dM(dm), mField(m_field), postProc(post_proc),
postProcBdy(post_proc_bdry), velocityFieldPtr(velocity_field_ptr),
x2FieldPtr(x2_field_ptr), geometryFieldPtr(geometry_field_ptr),
fieldEvalCoords(pass_field_eval_coords),
fieldEvalData(pass_field_eval_data){};
MoFEMErrorCode postProcess() {
// cerr << "wagawaga\n";
auto *simple = mField.getInterface<Simple>();
if (SPACE_DIM == 3) {
CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);
} else {
CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);
}
if (velocityFieldPtr->size1()) {
auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velocityFieldPtr);
auto t_x2_field = getFTensor1FromMat<SPACE_DIM>(*x2FieldPtr);
auto t_geom = getFTensor1FromMat<SPACE_DIM>(*geometryFieldPtr);
double u_x = t_x2_field(0) - t_geom(0);
double u_y = t_x2_field(1) - t_geom(1);
double u_z = t_x2_field(2) - t_geom(2);
MOFEM_LOG("SYNC", Sev::inform)
<< "Velocities x: " << t_vel(0) << " y: " << t_vel(1)
<< " z: " << t_vel(2) << "\n";
MOFEM_LOG("SYNC", Sev::inform) << "Displacement x: " << u_x
<< " y: " << u_y << " z: " << u_z << "\n";
}
std::regex((boost::format("%s(.*)") % "Data_Vertex").str()))) {
Range ents;
mField.get_moab().get_entities_by_dimension(m->getMeshset(), 0, ents,
true);
auto print_vets = [](boost::shared_ptr<FieldEntity> ent_ptr) {
if (!(ent_ptr->getPStatus() & PSTATUS_NOT_OWNED)) {
MOFEM_LOG("SYNC", Sev::inform)
<< "Velocities: " << ent_ptr->getEntFieldData()[0] << " "
<< ent_ptr->getEntFieldData()[1] << " "
<< ent_ptr->getEntFieldData()[2] << "\n";
}
};
CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
print_vets, "V", &ents);
}
MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::inform);
PetscBool print_volume = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-print_volume", &print_volume,
PETSC_NULL);
PetscBool print_skin = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-print_skin", &print_skin,
PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step",
&save_every_nth_step, PETSC_NULL);
if (ts_step % save_every_nth_step == 0) {
if (print_volume) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
if (print_skin) {
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdy);
CHKERR postProcBdy->writeFile(
"out_boundary_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<PostProcFaceEle> postProcBdy;
boost::shared_ptr<MatrixDouble> velocityFieldPtr;
boost::shared_ptr<MatrixDouble> x2FieldPtr;
boost::shared_ptr<MatrixDouble> geometryFieldPtr;
std::array<double, 3> fieldEvalCoords;
boost::shared_ptr<SetPtsData> fieldEvalData;
};
//! [Solve]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
auto calculate_stress_ops = [&](auto &pip) {
auto v_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("V", v_ptr));
auto X_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("x_1", x_ptr));
// Calculate unknown F
auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
"F", mat_H_tensor_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateDisplacement<SPACE_DIM>(x_ptr, X_ptr, u_ptr));
// Calculate P
auto mat_F_ptr = boost::make_shared<MatrixDouble>();
mat_F_ptr, mat_H_tensor_ptr));
PetscBool is_linear_elasticity = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",
&is_linear_elasticity, PETSC_NULL);
auto mat_P_ptr = boost::make_shared<MatrixDouble>();
if (is_linear_elasticity) {
mat_F_ptr));
} else {
auto inv_F = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pip.push_back(new OpInvertMatrix<SPACE_DIM>(mat_F_ptr, det_ptr, inv_F));
mat_F_ptr, inv_F, det_ptr));
}
auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
"V", mat_v_grad_ptr));
return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);
};
auto post_proc_boundary = [&]() {
auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
boundary_post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
auto op_loop_side =
new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
// push ops to side element, through op_loop_side operator
auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
boundary_mat_F_ptr, boundary_u_ptr] =
calculate_stress_ops(op_loop_side->getOpPtrVector());
boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
boundary_post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
boundary_post_proc_fe->getPostProcMesh(),
boundary_post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapMat{{"V", boundary_v_ptr},
{"GEOMETRY", boundary_X_ptr},
{"x", boundary_x_ptr},
{"U", boundary_u_ptr}},
OpPPMap::DataMapMat{{"FIRST_PIOLA", boundary_mat_P_ptr},
{"F", boundary_mat_F_ptr}},
)
);
return boundary_post_proc_fe;
};
// Add monitor to time solver
double rho = 1.;
CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-density", &rho, PETSC_NULL);
auto get_rho = [rho](const double, const double, const double) {
return rho;
};
SmartPetscObj<Mat> M; ///< Mass matrix
SmartPetscObj<KSP> ksp; ///< Linear solver
auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
tsPrePostProc = ts_pre_post_proc;
CHKERR MatZeroEntries(M);
boost::shared_ptr<DomainEle> vol_mass_ele(new DomainEle(mField));
vol_mass_ele->B = M;
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
vol_mass_ele->getRuleHook = integration_rule;
vol_mass_ele->getOpPtrVector().push_back(new OpMassV("V", "V", get_rho));
vol_mass_ele->getOpPtrVector().push_back(new OpMassF("F", "F"));
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_mass_ele);
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
auto lumpVec = createDMVector(simple->getDM());
CHKERR MatGetRowSum(M, lumpVec);
CHKERR MatZeroEntries(M);
CHKERR MatDiagonalSet(M, lumpVec, INSERT_VALUES);
// Create and septup KSP (linear solver), we need this to calculate g(t,u) =
// M^-1G(t,u)
ksp = createKSP(mField.get_comm());
CHKERR KSPSetOperators(ksp, M, M);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSetUp(ksp);
auto solve_boundary_for_g = [&]() {
if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {
CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) = false;
auto D =
smartVectorDuplicate(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecCopy(D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
}
};
pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =
solve_boundary_for_g;
ts = pipeline_mng->createTSEX(dm);
// Field eval
PetscBool field_eval_flag = PETSC_TRUE;
boost::shared_ptr<MatrixDouble> velocity_field_ptr;
boost::shared_ptr<MatrixDouble> geometry_field_ptr;
boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;
boost::shared_ptr<SetPtsData> field_eval_data;
std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};
int dim = 3;
CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
field_eval_coords.data(), &dim,
&field_eval_flag);
if (field_eval_flag) {
field_eval_data =
mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
if (SPACE_DIM == 3) {
field_eval_data, simple->getDomainFEName());
} else {
field_eval_data, simple->getDomainFEName());
}
field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
auto no_rule = [](int, int, int) { return -1; };
auto fe_ptr = field_eval_data->feMethodPtr.lock();
fe_ptr->getRuleHook = no_rule;
velocity_field_ptr = boost::make_shared<MatrixDouble>();
geometry_field_ptr = boost::make_shared<MatrixDouble>();
spatial_position_field_ptr = boost::make_shared<MatrixDouble>();
fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("V", velocity_field_ptr));
fe_ptr->getOpPtrVector().push_back(
geometry_field_ptr));
fe_ptr->getOpPtrVector().push_back(
"x_2", spatial_position_field_ptr));
}
auto post_proc_domain = [&]() {
auto post_proc_fe_vol = boost::make_shared<PostProcEle>(mField);
auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
boundary_mat_F_ptr, boundary_u_ptr] =
calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());
post_proc_fe_vol->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe_vol->getPostProcMesh(),
post_proc_fe_vol->getMapGaussPts(),
{},
{{"V", boundary_v_ptr},
{"GEOMETRY", boundary_X_ptr},
{"x", boundary_x_ptr},
{"U", boundary_u_ptr}},
{{"FIRST_PIOLA", boundary_mat_P_ptr}, {"F", boundary_mat_F_ptr}},
{}
)
);
return post_proc_fe_vol;
};
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(
SmartPetscObj<DM>(dm, true), mField, post_proc_domain(),
post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
geometry_field_ptr, field_eval_coords, field_eval_data);
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
null_fe, monitor_ptr);
double ftime = 1;
// CHKERR TSSetMaxTime(ts, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR TSSetSolution(ts, T);
CHKERR TSSetFromOptions(ts);
auto fb = mField.getInterface<FieldBlas>();
CHKERR TSSetPostStage(ts, TSPrePostProc::tsPostStage);
CHKERR TSSetPostStep(ts, TSPrePostProc::tsPostStep);
CHKERR TSSetPreStep(ts, TSPrePostProc::tsPreStep);
boost::shared_ptr<ForcesAndSourcesCore> null;
if (auto ptr = tsPrePostProc.lock()) {
ptr->fsRawPtr = this;
CHKERR TSSetUp(ts);
CHKERR TSSolve(ts, NULL);
CHKERR TSGetTime(ts, &ftime);
}
}
//! [Solve]
//! [Postprocess results]
PetscBool test_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
if (test_flg) {
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
double nrm2;
CHKERR VecNorm(T, NORM_2, &nrm2);
MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
constexpr double regression_value = 0.0194561;
if (fabs(nrm2 - regression_value) > 1e-2)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Regression test failed; wrong norm value.");
}
}
//! [Postprocess results]
//! [Check]
}
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: adolc_plasticity.cpp:97
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBodyForce
Definition: dynamic_first_order_con_law.cpp:66
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
young_modulus
constexpr double young_modulus
Definition: dynamic_first_order_con_law.cpp:94
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:93
OpMassV
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassV
Definition: dynamic_first_order_con_law.cpp:58
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpGradTimesTensor2
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesTensor2
Definition: dynamic_first_order_con_law.cpp:76
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
main
int main(int argc, char *argv[])
Definition: dynamic_first_order_con_law.cpp:1213
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:887
OpMassF
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM *SPACE_DIM > OpMassF
Definition: dynamic_first_order_con_law.cpp:60
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
rho
double rho
Definition: plastic.cpp:191
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
OpRhsTestPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpBaseTimesVector< 1, SPACE_DIM *SPACE_DIM, 1 > OpRhsTestPiola
Definition: dynamic_first_order_con_law.cpp:80
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
OpGradTimesPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesPiola
Definition: dynamic_first_order_con_law.cpp:84
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
SPACE_DIM
constexpr int SPACE_DIM
Definition: dynamic_first_order_con_law.cpp:24
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
TSPrePostProc::tsPostStage
static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
[Boundary condition]
Definition: dynamic_first_order_con_law.cpp:581
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:104
OpCalculateDeformationGradient
Definition: dynamic_first_order_con_law.cpp:339
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: adolc_plasticity.cpp:98
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: dynamic_first_order_con_law.cpp:405
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpCalculatePiolaIncompressibleNH
Definition: dynamic_first_order_con_law.cpp:267
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EntData
EntitiesFieldData::EntData EntData
Definition: dynamic_first_order_con_law.cpp:27
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: adolc_plasticity.cpp:99
poisson_ratio
constexpr double poisson_ratio
Definition: dynamic_first_order_con_law.cpp:95
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
PostProcEleByDim
Definition: adolc_plasticity.cpp:82
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
trace
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
Definition: dynamic_first_order_con_law.cpp:14
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
TSPrePostProc
Set of functions called by PETSc solver used to refine and update mesh.
Definition: dynamic_first_order_con_law.cpp:384
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
TSPrePostProc::tsPostStep
static MoFEMErrorCode tsPostStep(TS ts)
Definition: dynamic_first_order_con_law.cpp:614
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
OpCalculateDisplacement
Definition: dynamic_first_order_con_law.cpp:226
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::FieldEvaluatorInterface::SetPtsData
Definition: FieldEvaluator.hpp:29
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:228
double
convert.type
type
Definition: convert.py:64
DomainEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Definition: dynamic_first_order_con_law.cpp:28
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
OpCalculateFStab
Definition: dynamic_first_order_con_law.cpp:104
help
static char help[]
[Check]
Definition: dynamic_first_order_con_law.cpp:1211
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:880
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
OpCalculatePiola
Definition: dynamic_first_order_con_law.cpp:168
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128
SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: dynamic_first_order_con_law.cpp:37
MatrixFunction.hpp
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:226
lamme_lambda
double lamme_lambda
Definition: dynamic_first_order_con_law.cpp:99
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
FTensor::Index< 'i', SPACE_DIM >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
Range
DomainEleOp
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:276
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
Monitor
[Push operators to pipeline]
Definition: adolc_plasticity.cpp:333
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1060
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:264
TSPrePostProc::tsPreStep
static MoFEMErrorCode tsPreStep(TS ts)
Definition: dynamic_first_order_con_law.cpp:628
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:235
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
LinMomTimeScale
Definition: dynamic_first_order_con_law.cpp:407
QUIET
@ QUIET
Definition: definitions.h:208
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:575
dt
double dt
Definition: heat_method.cpp:26
MoFEM::SmartPetscObj< Mat >
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182