constexpr
double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1));
};
constexpr
double third = boost::math::constants::third<double>();
return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2));
};
};
};
GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
constexpr
double omega = 1.;
template <int DIM_0, int DIM_1>
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) {}
const size_t nb_gauss_pts = getGaussPts().size2();
defGradStabPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
defGradStabPtr->clear();
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);
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) {
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;
};
template <int DIM_0, int DIM_1>
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) {}
const size_t nb_gauss_pts = getGaussPts().size2();
firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
firstPiolaPtr->clear();
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_dk *
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>
boost::shared_ptr<MatrixDouble> reference_pos_ptr,
boost::shared_ptr<MatrixDouble> u_ptr)
xPtr(spatial_pos_ptr), XPtr(reference_pos_ptr), uPtr(u_ptr) {}
const size_t nb_gauss_pts = getGaussPts().size2();
uPtr->resize(DIM, nb_gauss_pts, false);
uPtr->clear();
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) {}
const size_t nb_gauss_pts = getGaussPts().size2();
firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
firstPiolaPtr->clear();
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) {
t_P(
i,
j) = bulk_mod * (t_det - 1.) * t_det * t_inv_F(
j,
i);
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) {}
const size_t nb_gauss_pts = getGaussPts().size2();
defGradPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
defGradPtr->clear();
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;
++t_H;
}
}
private:
boost::shared_ptr<MatrixDouble> gradTensorPtr;
boost::shared_ptr<MatrixDouble> defGradPtr;
};
PetscInt stageindex,
Vec *Y);
};
double getScale(const double time) {
};
};
};
private:
double getScale(const double time) { return 0.001 * sin(0.1 * time); };
};
double getScale(const double time) { return 0.001; };
};
};
}
}
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
break;
}
auto project_ho_geometry = [&]() {
CHKERR mField.loop_dofs(
"x_1", ent_method_x);
CHKERR mField.loop_dofs(
"x_2", ent_method_x_2);
return mField.loop_dofs("GEOMETRY", ent_method);
};
}
auto time_scale = boost::make_shared<TimeScale>();
PetscBool sin_time_function = PETSC_FALSE;
&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);
};
simple->getProblemName(),
"V");
auto get_pre_proc_hook = [&]() {
mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
};
pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
}
PetscInt stageindex,
Vec *Y) {
auto &m_field = ptr->fsRawPtr->mField;
double time;
PetscInt num_stages;
CHKERR TSGetStages(ts, &num_stages, &stage_solutions);
PetscPrintf(PETSC_COMM_WORLD, "Check timestep %d time %e dt %e\n",
const double inv_num_step = (
double)num_stages;
CHKERR fb->fieldCopy(1.,
"x_1",
"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");
}
}
auto &m_field = ptr->fsRawPtr->mField;
double time;
}
}
auto &m_field = ptr->fsRawPtr->mField;
double time;
int step_num;
CHKERR TSGetStepNumber(ts, &step_num);
}
}
t_source(0) = 0.1;
t_source(1) = 1.;
return t_source;
};
auto get_time_scale = [this](const double time) {
return sin(time *
omega * M_PI);
};
auto apply_rhs = [&](auto &pip) {
"GEOMETRY");
auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
"V", mat_v_grad_ptr));
auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
auto set_body_force = [&]() {
auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
double unit_weight = 0.;
PETSC_NULL);
t_force(1) = -unit_weight;
t_force(2) = unit_weight;
}
};
[](double, double, double) { return 1.; }));
auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
"F", mat_H_tensor_ptr));
double tau = 0.2;
double xi = 0.;
};
return -1.;
};
auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
"F_dot", mat_dot_F_tensor_ptr));
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;
&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(
mat_F_stab_ptr, inv_F, det_ptr));
}
};
CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
};
}
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){};
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
} else {
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
}
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);
<< "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()))) {
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)) {
<< "Velocities: " << ent_ptr->getEntFieldData()[0] << " "
<< ent_ptr->getEntFieldData()[1] << " "
<< ent_ptr->getEntFieldData()[2] << "\n";
}
};
print_vets, "V", &ents);
}
PetscBool print_volume = PETSC_FALSE;
PETSC_NULL);
PetscBool print_skin = PETSC_FALSE;
PETSC_NULL);
if (print_volume) {
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
if (print_skin) {
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;
};
auto calculate_stress_ops = [&](auto &pip) {
auto v_ptr = boost::make_shared<MatrixDouble>();
auto X_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
auto x_ptr = boost::make_shared<MatrixDouble>();
auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
"F", mat_H_tensor_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
auto mat_F_ptr = boost::make_shared<MatrixDouble>();
mat_F_ptr, mat_H_tensor_ptr));
PetscBool is_linear_elasticity = PETSC_TRUE;
&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>();
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 =
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(
boundary_post_proc_fe->getPostProcMesh(),
boundary_post_proc_fe->getMapGaussPts(),
{"GEOMETRY", boundary_X_ptr},
{"x", boundary_x_ptr},
{"U", boundary_u_ptr}},
{"F", boundary_mat_F_ptr}},
)
);
return boundary_post_proc_fe;
};
};
auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
boost::shared_ptr<DomainEle> vol_mass_ele(
new DomainEle(mField));
};
vol_mass_ele->getOpPtrVector().push_back(
new OpMassV(
"V",
"V", get_rho));
vol_mass_ele->getOpPtrVector().push_back(
new OpMassF(
"F",
"F"));
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatDiagonalSet(
M, lumpVec, INSERT_VALUES);
CHKERR KSPSetFromOptions(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;
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);
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;
field_eval_coords.data(), &dim,
&field_eval_flag);
if (field_eval_flag) {
field_eval_data =
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(
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(
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>(
post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
geometry_field_ptr, field_eval_coords, field_eval_data);
null_fe, monitor_ptr);
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
boost::shared_ptr<ForcesAndSourcesCore> null;
ptr->fsRawPtr = this;
}
}
PetscBool test_flg = PETSC_FALSE;
if (test_flg) {
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)
"Regression test failed; wrong norm value.");
}
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
}
}