12 using namespace MoFEM;
15 constexpr
double third = boost::math::constants::third<double>();
16 return (t_stress(0, 0) + t_stress(1, 1));
20 return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2));
65 GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
102 template <
int DIM_0,
int DIM_1>
105 boost::shared_ptr<MatrixDouble> def_grad_stab_ptr,
106 boost::shared_ptr<MatrixDouble> def_grad_dot_ptr,
107 double tau_F_ptr,
double xi_F_ptr,
108 boost::shared_ptr<MatrixDouble> grad_x_ptr,
109 boost::shared_ptr<MatrixDouble> grad_vel_ptr)
111 defGradPtr(def_grad_ptr), defGradStabPtr(def_grad_stab_ptr),
112 defGradDotPtr(def_grad_dot_ptr), tauFPtr(tau_F_ptr), xiF(xi_F_ptr),
113 gradxPtr(grad_x_ptr), gradVelPtr(grad_vel_ptr) {}
123 const size_t nb_gauss_pts = getGaussPts().size2();
125 defGradStabPtr->resize(DIM_0 * DIM_1, nb_gauss_pts,
false);
126 defGradStabPtr->clear();
129 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
130 auto t_Fstab = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradStabPtr);
131 auto t_F_dot = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradDotPtr);
134 auto tau_F = tauFPtr;
136 auto t_gradx = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradxPtr);
137 auto t_gradVel = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradVelPtr);
139 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
141 t_Fstab(
i,
j) = t_F(
i,
j) + tau_F * (t_gradVel(
i,
j) - t_F_dot(
i,
j)) +
142 xi_F * (t_gradx(
i,
j) - t_F(
i,
j));
166 template <
int DIM_0,
int DIM_1>
170 boost::shared_ptr<MatrixDouble> first_piola_ptr,
171 boost::shared_ptr<MatrixDouble> def_grad_ptr)
173 shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
174 lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
175 defGradPtr(def_grad_ptr) {}
189 const size_t nb_gauss_pts = getGaussPts().size2();
192 firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts,
false);
193 firstPiolaPtr->clear();
196 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
197 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
198 const double two_o_three = 2. / 3.;
199 const double trace_t_dk = DIM_0;
200 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
202 t_P(
i,
j) = shearModulus * (t_F(
i,
j) + t_F(
j,
i) - 2. *
t_kd(
i,
j) -
204 two_o_three * trace_t_dk *
t_kd(
i,
j)) +
206 bulkModulus * trace_t_dk *
t_kd(
i,
j);
227 boost::shared_ptr<MatrixDouble> reference_pos_ptr,
228 boost::shared_ptr<MatrixDouble> u_ptr)
230 xPtr(spatial_pos_ptr), XPtr(reference_pos_ptr), uPtr(u_ptr) {}
239 const size_t nb_gauss_pts = getGaussPts().size2();
241 uPtr->resize(DIM, nb_gauss_pts,
false);
245 auto t_x = getFTensor1FromMat<DIM>(*xPtr);
246 auto t_X = getFTensor1FromMat<DIM>(*XPtr);
247 auto t_u = getFTensor1FromMat<DIM>(*uPtr);
248 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
250 t_u(
i) = t_x(
i) - t_X(
i);
260 boost::shared_ptr<MatrixDouble>
xPtr;
261 boost::shared_ptr<MatrixDouble>
XPtr;
262 boost::shared_ptr<MatrixDouble>
uPtr;
265 template <
int DIM_0,
int DIM_1>
269 double shear_modulus,
double bulk_modulus,
double m_u,
270 double lambda_lamme, boost::shared_ptr<MatrixDouble> first_piola_ptr,
271 boost::shared_ptr<MatrixDouble> def_grad_ptr,
272 boost::shared_ptr<MatrixDouble> inv_def_grad_ptr,
273 boost::shared_ptr<VectorDouble> det)
275 shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
276 lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
277 defGradPtr(def_grad_ptr), invDefGradPtr(inv_def_grad_ptr), dEt(det) {}
292 const size_t nb_gauss_pts = getGaussPts().size2();
295 firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts,
false);
296 firstPiolaPtr->clear();
299 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
300 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
301 auto t_inv_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*invDefGradPtr);
302 auto t_det = getFTensor0FromVec<1>(*dEt);
303 const double two_o_three = 2. / 3.;
304 const double one_o_three = 1. / 3.;
305 const double bulk_mod = bulkModulus;
306 const double shear_mod = shearModulus;
307 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
311 t_P(
i,
j) = bulk_mod * (t_det - 1.) * t_det * t_inv_F(
j,
i);
314 shear_mod * pow(t_det, two_o_three) *
315 (t_F(
i,
j) - one_o_three * (t_F(
l,
k) * t_F(
l,
k)) * t_inv_F(
j,
i));
334 boost::shared_ptr<VectorDouble>
dEt;
337 template <
int DIM_0,
int DIM_1>
341 boost::shared_ptr<MatrixDouble> def_grad_ptr,
342 boost::shared_ptr<MatrixDouble> grad_tensor_ptr)
344 defGradPtr(def_grad_ptr), gradTensorPtr(grad_tensor_ptr) {}
357 const size_t nb_gauss_pts = getGaussPts().size2();
360 defGradPtr->resize(DIM_0 * DIM_1, nb_gauss_pts,
false);
364 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
365 auto t_H = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradTensorPtr);
366 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
399 PetscInt stageindex,
Vec *Y);
438 double getScale(
const double time) {
return 0.001 * sin(0.1 * time); };
443 double getScale(
const double time) {
return 0.001; };
452 CHKERR boundaryCondition();
476 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
477 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
478 PetscInt choice_base_value = AINSWORTH;
480 LASBASETOPT, &choice_base_value, PETSC_NULL);
483 switch (choice_base_value) {
487 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
492 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
519 auto project_ho_geometry = [&]() {
521 CHKERR mField.loop_dofs(
"x_1", ent_method_x);
523 CHKERR mField.loop_dofs(
"x_2", ent_method_x_2);
526 return mField.loop_dofs(
"GEOMETRY", ent_method);
528 CHKERR project_ho_geometry();
541 auto time_scale = boost::make_shared<TimeScale>();
543 PetscBool sin_time_function = PETSC_FALSE;
545 &sin_time_function, PETSC_NULL);
547 if (sin_time_function)
548 time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();
550 time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();
552 pipeline_mng->getBoundaryExplicitRhsFE().reset();
554 pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE},
"GEOMETRY");
557 pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField,
"V",
558 {time_scale},
"FORCE", Sev::inform);
568 simple->getProblemName(),
"V");
570 auto get_pre_proc_hook = [&]() {
572 mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
574 pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
581 PetscInt stageindex,
Vec *Y) {
585 auto &m_field = ptr->fsRawPtr->mField;
587 auto fb = m_field.getInterface<
FieldBlas>();
591 CHKERR TSGetTime(ts, &time);
593 Vec *stage_solutions;
595 CHKERR TSGetStages(ts, &num_stages, &stage_solutions);
596 PetscPrintf(PETSC_COMM_WORLD,
"Check timestep %d time %e dt %e\n",
597 num_stages, time,
dt);
599 const double inv_num_step = (
double)num_stages;
600 CHKERR fb->fieldCopy(1.,
"x_1",
"x_2");
601 CHKERR fb->fieldAxpy(
dt,
"V",
"x_2");
602 CHKERR fb->fieldCopy(1.,
"x_2",
"x_1");
604 CHKERR fb->fieldCopy(-inv_num_step /
dt,
"F_0",
"F_dot");
605 CHKERR fb->fieldAxpy(inv_num_step /
dt,
"F",
"F_dot");
606 CHKERR fb->fieldCopy(1.,
"F",
"F_0");
615 auto &m_field = ptr->fsRawPtr->mField;
619 CHKERR TSGetTime(ts, &time);
628 auto &m_field = ptr->fsRawPtr->mField;
632 CHKERR TSGetTime(ts, &time);
634 CHKERR TSGetStepNumber(ts, &step_num);
652 auto get_time_scale = [
this](
const double time) {
653 return sin(time *
omega * M_PI);
656 auto apply_rhs = [&](
auto &pip) {
663 auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
665 "V", mat_v_grad_ptr));
667 auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
668 gravity_vector_ptr->resize(
SPACE_DIM, 1);
669 auto set_body_force = [&]() {
672 auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
673 double unit_weight = 0.;
678 t_force(1) = -unit_weight;
680 t_force(2) = unit_weight;
686 pip.push_back(
new OpBodyForce(
"V", gravity_vector_ptr,
687 [](
double,
double,
double) {
return 1.; }));
690 auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
692 "F", mat_H_tensor_ptr));
709 auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
711 "F_dot", mat_dot_F_tensor_ptr));
714 auto mat_x_grad_ptr = boost::make_shared<MatrixDouble>();
716 "x_2", mat_x_grad_ptr));
718 auto mat_F_tensor_ptr = boost::make_shared<MatrixDouble>();
720 mat_F_tensor_ptr, mat_H_tensor_ptr));
722 auto mat_F_stab_ptr = boost::make_shared<MatrixDouble>();
724 mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,
725 mat_x_grad_ptr, mat_v_grad_ptr));
727 PetscBool is_linear_elasticity = PETSC_TRUE;
729 &is_linear_elasticity, PETSC_NULL);
731 auto mat_P_stab_ptr = boost::make_shared<MatrixDouble>();
732 if (is_linear_elasticity) {
737 auto inv_F = boost::make_shared<MatrixDouble>();
738 auto det_ptr = boost::make_shared<VectorDouble>();
746 mat_F_stab_ptr, inv_F, det_ptr));
756 CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
777 boost::shared_ptr<PostProcEle> post_proc,
778 boost::shared_ptr<PostProcFaceEle> post_proc_bdry,
779 boost::shared_ptr<MatrixDouble> velocity_field_ptr,
780 boost::shared_ptr<MatrixDouble> x2_field_ptr,
781 boost::shared_ptr<MatrixDouble> geometry_field_ptr,
782 std::array<double, 3> pass_field_eval_coords,
783 boost::shared_ptr<SetPtsData> pass_field_eval_data)
784 : dM(dm), mField(m_field), postProc(post_proc),
785 postProcBdy(post_proc_bdry), velocityFieldPtr(velocity_field_ptr),
786 x2FieldPtr(x2_field_ptr), geometryFieldPtr(geometry_field_ptr),
787 fieldEvalCoords(pass_field_eval_coords),
788 fieldEvalData(pass_field_eval_data){};
797 fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
802 fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
807 if (velocityFieldPtr->size1()) {
808 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velocityFieldPtr);
809 auto t_x2_field = getFTensor1FromMat<SPACE_DIM>(*x2FieldPtr);
810 auto t_geom = getFTensor1FromMat<SPACE_DIM>(*geometryFieldPtr);
812 double u_x = t_x2_field(0) - t_geom(0);
813 double u_y = t_x2_field(1) - t_geom(1);
814 double u_z = t_x2_field(2) - t_geom(2);
817 <<
"Velocities x: " << t_vel(0) <<
" y: " << t_vel(1)
818 <<
" z: " << t_vel(2) <<
"\n";
819 MOFEM_LOG(
"SYNC", Sev::inform) <<
"Displacement x: " << u_x
820 <<
" y: " << u_y <<
" z: " << u_z <<
"\n";
824 std::regex((boost::format(
"%s(.*)") %
"Data_Vertex").str()))) {
826 mField.
get_moab().get_entities_by_dimension(
m->getMeshset(), 0, ents,
828 auto print_vets = [](boost::shared_ptr<FieldEntity> ent_ptr) {
830 if (!(ent_ptr->getPStatus() & PSTATUS_NOT_OWNED)) {
832 <<
"Velocities: " << ent_ptr->getEntFieldData()[0] <<
" "
833 << ent_ptr->getEntFieldData()[1] <<
" "
834 << ent_ptr->getEntFieldData()[2] <<
"\n";
839 print_vets,
"V", &ents);
843 PetscBool print_volume = PETSC_FALSE;
847 PetscBool print_skin = PETSC_FALSE;
858 CHKERR postProc->writeFile(
859 "out_step_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
864 CHKERR postProcBdy->writeFile(
865 "out_boundary_" + boost::lexical_cast<std::string>(ts_step) +
889 auto dm =
simple->getDM();
891 auto calculate_stress_ops = [&](
auto &pip) {
894 auto v_ptr = boost::make_shared<MatrixDouble>();
896 auto X_ptr = boost::make_shared<MatrixDouble>();
900 auto x_ptr = boost::make_shared<MatrixDouble>();
904 auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
906 "F", mat_H_tensor_ptr));
908 auto u_ptr = boost::make_shared<MatrixDouble>();
912 auto mat_F_ptr = boost::make_shared<MatrixDouble>();
914 mat_F_ptr, mat_H_tensor_ptr));
916 PetscBool is_linear_elasticity = PETSC_TRUE;
918 &is_linear_elasticity, PETSC_NULL);
920 auto mat_P_ptr = boost::make_shared<MatrixDouble>();
921 if (is_linear_elasticity) {
926 auto inv_F = boost::make_shared<MatrixDouble>();
927 auto det_ptr = boost::make_shared<VectorDouble>();
933 mat_F_ptr, inv_F, det_ptr));
936 auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
938 "V", mat_v_grad_ptr));
940 return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);
943 auto post_proc_boundary = [&]() {
944 auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
947 boundary_post_proc_fe->getOpPtrVector(), {},
"GEOMETRY");
951 auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
952 boundary_mat_F_ptr, boundary_u_ptr] =
953 calculate_stress_ops(op_loop_side->getOpPtrVector());
954 boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
958 boundary_post_proc_fe->getOpPtrVector().push_back(
962 boundary_post_proc_fe->getPostProcMesh(),
963 boundary_post_proc_fe->getMapGaussPts(),
968 {
"GEOMETRY", boundary_X_ptr},
969 {
"x", boundary_x_ptr},
970 {
"U", boundary_u_ptr}},
973 {
"F", boundary_mat_F_ptr}},
980 return boundary_post_proc_fe;
994 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
1000 boost::shared_ptr<DomainEle> vol_mass_ele(
new DomainEle(mField));
1002 vol_mass_ele->B =
M;
1010 vol_mass_ele->getOpPtrVector().push_back(
new OpMassV(
"V",
"V", get_rho));
1011 vol_mass_ele->getOpPtrVector().push_back(
new OpMassF(
"F",
"F"));
1014 CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
1015 CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
1018 CHKERR MatGetRowSum(
M, lumpVec);
1021 CHKERR MatDiagonalSet(
M, lumpVec, INSERT_VALUES);
1027 CHKERR KSPSetFromOptions(ksp);
1030 auto solve_boundary_for_g = [&]() {
1032 if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {
1034 CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1035 ADD_VALUES, SCATTER_REVERSE);
1036 CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1037 ADD_VALUES, SCATTER_REVERSE);
1038 CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1039 CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1040 *(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) =
false;
1044 CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
D);
1045 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1046 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1047 CHKERR VecCopy(
D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1053 pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =
1054 solve_boundary_for_g;
1057 ts = pipeline_mng->createTSEX(dm);
1060 PetscBool field_eval_flag = PETSC_TRUE;
1061 boost::shared_ptr<MatrixDouble> velocity_field_ptr;
1062 boost::shared_ptr<MatrixDouble> geometry_field_ptr;
1063 boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;
1064 boost::shared_ptr<SetPtsData> field_eval_data;
1066 std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};
1069 field_eval_coords.data(), &dim,
1072 if (field_eval_flag) {
1077 field_eval_data,
simple->getDomainFEName());
1080 field_eval_data,
simple->getDomainFEName());
1083 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1085 auto no_rule = [](
int,
int,
int) {
return -1; };
1087 auto fe_ptr = field_eval_data->feMethodPtr.lock();
1088 fe_ptr->getRuleHook = no_rule;
1089 velocity_field_ptr = boost::make_shared<MatrixDouble>();
1090 geometry_field_ptr = boost::make_shared<MatrixDouble>();
1091 spatial_position_field_ptr = boost::make_shared<MatrixDouble>();
1092 fe_ptr->getOpPtrVector().push_back(
1094 fe_ptr->getOpPtrVector().push_back(
1096 geometry_field_ptr));
1097 fe_ptr->getOpPtrVector().push_back(
1099 "x_2", spatial_position_field_ptr));
1102 auto post_proc_domain = [&]() {
1103 auto post_proc_fe_vol = boost::make_shared<PostProcEle>(mField);
1107 auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
1108 boundary_mat_F_ptr, boundary_u_ptr] =
1109 calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());
1111 post_proc_fe_vol->getOpPtrVector().push_back(
1115 post_proc_fe_vol->getPostProcMesh(),
1116 post_proc_fe_vol->getMapGaussPts(),
1120 {{
"V", boundary_v_ptr},
1121 {
"GEOMETRY", boundary_X_ptr},
1122 {
"x", boundary_x_ptr},
1123 {
"U", boundary_u_ptr}},
1125 {{
"FIRST_PIOLA", boundary_mat_P_ptr}, {
"F", boundary_mat_F_ptr}},
1132 return post_proc_fe_vol;
1135 boost::shared_ptr<FEMethod> null_fe;
1136 auto monitor_ptr = boost::make_shared<Monitor>(
1138 post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
1139 geometry_field_ptr, field_eval_coords, field_eval_data);
1142 null_fe, monitor_ptr);
1146 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1151 CHKERR TSSetSolution(ts, T);
1152 CHKERR TSSetFromOptions(ts);
1158 boost::shared_ptr<ForcesAndSourcesCore>
null;
1161 ptr->fsRawPtr =
this;
1163 CHKERR TSSolve(ts, NULL);
1164 CHKERR TSGetTime(ts, &ftime);
1174 PetscBool test_flg = PETSC_FALSE;
1182 CHKERR VecNorm(T, NORM_2, &nrm2);
1183 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Regression norm " << nrm2;
1184 constexpr
double regression_value = 0.0194561;
1185 if (fabs(nrm2 - regression_value) > 1e-2)
1187 "Regression test failed; wrong norm value.");
1205 const char param_file[] =
"param_file.petsc";
1209 auto core_log = logging::core::get();
1211 LogManager::createSink(LogManager::getStrmWorld(),
"EXAMPLE"));
1212 LogManager::setLog(
"EXAMPLE");
1218 DMType dm_name =
"DMMOFEM";