12 using namespace MoFEM;
81 CHKERR boundaryCondition();
96 char meshFileName[255];
98 meshFileName, 255, PETSC_NULL);
110 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
111 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
112 PetscInt choice_base_value = AINSWORTH;
114 LASBASETOPT, &choice_base_value, PETSC_NULL);
117 switch (choice_base_value) {
121 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
126 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
150 auto time_scale = boost::make_shared<TimeScale>();
161 pipeline_mng->getOpBoundaryRhsPipeline(), mField,
"U", {time_scale},
162 "FORCE", Sev::inform);
166 pipeline_mng->getOpDomainRhsPipeline(), mField,
"U", {time_scale},
167 "BODY_FORCE", Sev::inform);
170 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
172 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
174 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
177 simple->getProblemName(),
"U");
180 CHKERR bc_mng->addBlockDOFsToMPCs(
simple->getProblemName(),
"U");
192 auto add_domain_ops_lhs = [&](
auto &pip) {
195 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
196 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
200 auto add_domain_ops_rhs = [&](
auto &pip) {
203 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
204 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
208 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
209 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
223 : dM(dm), postProc(post_proc){};
230 postProc->mField, postProc->getPostProcMesh(), {
"U"});
232 CHKERR postProc->writeFile(
233 "out_step_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
249 auto dm =
simple->getDM();
250 auto ts = pipeline_mng->createTSIM();
253 auto create_post_proc_fe = [dm,
this,
simple]() {
254 auto post_proc_ele_domain = [dm,
this](
auto &pip_domain) {
256 auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
257 mField, pip_domain,
"U",
"MAT_ELASTIC", Sev::inform);
261 auto post_proc_fe_bdy = boost::make_shared<PostProcEleBdy>(mField);
262 auto u_ptr = boost::make_shared<MatrixDouble>();
263 post_proc_fe_bdy->getOpPtrVector().push_back(
267 auto common_ptr = post_proc_ele_domain(op_loop_side->getOpPtrVector());
268 post_proc_fe_bdy->getOpPtrVector().push_back(op_loop_side);
272 post_proc_fe_bdy->getOpPtrVector().push_back(
276 post_proc_fe_bdy->getPostProcMesh(),
277 post_proc_fe_bdy->getMapGaussPts(),
283 {{
"GRAD", common_ptr->matGradPtr},
284 {
"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
291 return post_proc_fe_bdy;
294 auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
297 auto pre_proc_ptr = boost::make_shared<FEMethod>();
298 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
299 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
301 auto time_scale = boost::make_shared<TimeScale>();
303 auto get_bc_hook_rhs = [
this, pre_proc_ptr, time_scale]() {
306 {time_scale},
false)();
310 pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
312 auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
315 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
317 mField, post_proc_rhs_ptr, 1.)();
321 auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
324 mField, post_proc_lhs_ptr, 1.)();
328 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
329 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
333 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
334 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
335 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
336 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
342 CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
344 auto create_monitor_fe = [dm](
auto &&post_proc_fe) {
345 return boost::make_shared<Monitor>(dm, post_proc_fe);
349 boost::shared_ptr<FEMethod> null_fe;
350 auto monitor_ptr = create_monitor_fe(create_post_proc_fe());
352 null_fe, monitor_ptr);
356 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
357 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
362 CHKERR TSSetFromOptions(ts);
365 CHKERR TSGetTime(ts, &ftime);
367 PetscInt steps, snesfails, rejects, nonlinits, linits;
368 CHKERR TSGetStepNumber(ts, &steps);
369 CHKERR TSGetSNESFailures(ts, &snesfails);
370 CHKERR TSGetStepRejections(ts, &rejects);
371 CHKERR TSGetSNESIterations(ts, &nonlinits);
372 CHKERR TSGetKSPIterations(ts, &linits);
374 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
376 steps, rejects, snesfails, ftime, nonlinits, linits);
387 auto dm =
simple->getDM();
393 CHKERR VecNorm(T, NORM_2, &nrm2);
394 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Solution norm " << nrm2;
396 auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
398 auto post_proc_norm_rule_hook = [](
int,
int,
int p) ->
int {
return 2 * p; };
399 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
402 post_proc_norm_fe->getOpPtrVector(), {H1});
404 enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
407 (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
408 CHKERR VecZeroEntries(norms_vec);
410 auto u_ptr = boost::make_shared<MatrixDouble>();
411 post_proc_norm_fe->getOpPtrVector().push_back(
414 post_proc_norm_fe->getOpPtrVector().push_back(
417 auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
418 mField, post_proc_norm_fe->getOpPtrVector(),
"U",
"MAT_ELASTIC",
421 post_proc_norm_fe->getOpPtrVector().push_back(
423 common_ptr->getMatFirstPiolaStress(), norms_vec, PIOLA_NORM));
428 CHKERR VecAssemblyBegin(norms_vec);
429 CHKERR VecAssemblyEnd(norms_vec);
432 if (mField.get_comm_rank() == 0) {
434 CHKERR VecGetArrayRead(norms_vec, &norms);
436 <<
"norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
438 <<
"norm_piola: " << std::scientific << std::sqrt(norms[PIOLA_NORM]);
439 CHKERR VecRestoreArrayRead(norms_vec, &norms);
449 PetscInt test_nb = 0;
450 PetscBool test_flg = PETSC_FALSE;
459 CHKERR VecNorm(T, NORM_2, &nrm2);
460 MOFEM_LOG(
"EXAMPLE", Sev::verbose) <<
"Regression norm " << nrm2;
461 double regression_value = 0;
464 regression_value = 1.02789;
467 regression_value = 1.8841e+00;
470 regression_value = 1.8841e+00;
473 regression_value = 4.9625e+00;
476 regression_value = 6.6394e+00;
479 regression_value = 4.98764e+00;
482 regression_value = 4.9473e+00;
485 regression_value = 2.5749e-01;
492 if (fabs(nrm2 - regression_value) > 1e-2)
494 "Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
508 static char help[] =
"...\n\n";
510 int main(
int argc,
char *argv[]) {
513 const char param_file[] =
"param_file.petsc";
517 auto core_log = logging::core::get();
526 DMType dm_name =
"DMMOFEM";