12 using namespace MoFEM;
82 CHKERR boundaryCondition();
97 char meshFileName[255];
99 meshFileName, 255, PETSC_NULL);
111 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
112 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
113 PetscInt choice_base_value = AINSWORTH;
115 LASBASETOPT, &choice_base_value, PETSC_NULL);
118 switch (choice_base_value) {
122 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
127 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
151 auto time_scale = boost::make_shared<TimeScale>();
162 pipeline_mng->getOpBoundaryRhsPipeline(), mField,
"U", {time_scale},
163 "FORCE", Sev::inform);
167 pipeline_mng->getOpDomainRhsPipeline(), mField,
"U", {time_scale},
168 "BODY_FORCE", Sev::inform);
171 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
173 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
175 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
178 simple->getProblemName(),
"U");
181 CHKERR bc_mng->addBlockDOFsToMPCs(
simple->getProblemName(),
"U");
193 auto add_domain_ops_lhs = [&](
auto &pip) {
196 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
197 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
201 auto add_domain_ops_rhs = [&](
auto &pip) {
204 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
205 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
209 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
210 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
224 : dM(dm), postProc(post_proc){};
231 postProc->mField, postProc->getPostProcMesh(), {
"U"});
233 CHKERR postProc->writeFile(
234 "out_step_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
250 auto dm =
simple->getDM();
251 auto ts = pipeline_mng->createTSIM();
254 auto create_post_proc_fe = [dm,
this,
simple]() {
255 auto post_proc_ele_domain = [dm,
this](
auto &pip_domain) {
257 auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
258 mField, pip_domain,
"U",
"MAT_ELASTIC", Sev::inform);
262 auto post_proc_fe_bdy = boost::make_shared<PostProcEleBdy>(mField);
263 auto u_ptr = boost::make_shared<MatrixDouble>();
264 post_proc_fe_bdy->getOpPtrVector().push_back(
268 auto common_ptr = post_proc_ele_domain(op_loop_side->getOpPtrVector());
269 post_proc_fe_bdy->getOpPtrVector().push_back(op_loop_side);
273 post_proc_fe_bdy->getOpPtrVector().push_back(
277 post_proc_fe_bdy->getPostProcMesh(),
278 post_proc_fe_bdy->getMapGaussPts(),
284 {{
"GRAD", common_ptr->matGradPtr},
285 {
"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
292 return post_proc_fe_bdy;
295 auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
298 auto pre_proc_ptr = boost::make_shared<FEMethod>();
299 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
300 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
302 auto time_scale = boost::make_shared<TimeScale>();
304 auto get_bc_hook_rhs = [
this, pre_proc_ptr, time_scale]() {
307 {time_scale},
false)();
311 pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
313 auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
316 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
318 mField, post_proc_rhs_ptr, 1.)();
322 auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
325 mField, post_proc_lhs_ptr, 1.)();
329 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
330 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
334 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
335 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
336 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
337 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
343 CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
345 auto create_monitor_fe = [dm](
auto &&post_proc_fe) {
346 return boost::make_shared<Monitor>(dm, post_proc_fe);
350 boost::shared_ptr<FEMethod> null_fe;
351 auto monitor_ptr = create_monitor_fe(create_post_proc_fe());
353 null_fe, monitor_ptr);
357 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
358 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
363 CHKERR TSSetFromOptions(ts);
366 CHKERR TSGetTime(ts, &ftime);
368 PetscInt steps, snesfails, rejects, nonlinits, linits;
369 CHKERR TSGetStepNumber(ts, &steps);
370 CHKERR TSGetSNESFailures(ts, &snesfails);
371 CHKERR TSGetStepRejections(ts, &rejects);
372 CHKERR TSGetSNESIterations(ts, &nonlinits);
373 CHKERR TSGetKSPIterations(ts, &linits);
375 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
377 steps, rejects, snesfails, ftime, nonlinits, linits);
388 auto dm =
simple->getDM();
394 CHKERR VecNorm(T, NORM_2, &nrm2);
395 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Solution norm " << nrm2;
397 auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
399 auto post_proc_norm_rule_hook = [](
int,
int,
int p) ->
int {
return 2 * p; };
400 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
403 post_proc_norm_fe->getOpPtrVector(), {H1});
405 enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
408 (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
409 CHKERR VecZeroEntries(norms_vec);
411 auto u_ptr = boost::make_shared<MatrixDouble>();
412 post_proc_norm_fe->getOpPtrVector().push_back(
415 post_proc_norm_fe->getOpPtrVector().push_back(
418 auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
419 mField, post_proc_norm_fe->getOpPtrVector(),
"U",
"MAT_ELASTIC",
422 post_proc_norm_fe->getOpPtrVector().push_back(
424 common_ptr->getMatFirstPiolaStress(), norms_vec, PIOLA_NORM));
429 CHKERR VecAssemblyBegin(norms_vec);
430 CHKERR VecAssemblyEnd(norms_vec);
433 if (mField.get_comm_rank() == 0) {
435 CHKERR VecGetArrayRead(norms_vec, &norms);
437 <<
"norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
439 <<
"norm_piola: " << std::scientific << std::sqrt(norms[PIOLA_NORM]);
440 CHKERR VecRestoreArrayRead(norms_vec, &norms);
450 PetscInt test_nb = 0;
451 PetscBool test_flg = PETSC_FALSE;
460 CHKERR VecNorm(T, NORM_2, &nrm2);
461 MOFEM_LOG(
"EXAMPLE", Sev::verbose) <<
"Regression norm " << nrm2;
462 double regression_value = 0;
465 regression_value = 1.02789;
468 regression_value = 1.8841e+00;
471 regression_value = 1.8841e+00;
474 regression_value = 4.9625e+00;
477 regression_value = 6.6394e+00;
480 regression_value = 4.98764e+00;
483 regression_value = 4.9473e+00;
486 regression_value = 2.5749e-01;
493 if (fabs(nrm2 - regression_value) > 1e-2)
495 "Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
509 static char help[] =
"...\n\n";
511 int main(
int argc,
char *argv[]) {
514 const char param_file[] =
"param_file.petsc";
518 auto core_log = logging::core::get();
527 DMType dm_name =
"DMMOFEM";