v0.14.0
Public Member Functions | Private Member Functions | Private Attributes | List of all members
Incompressible Struct Reference
Collaboration diagram for Incompressible:
[legend]

Public Member Functions

 Incompressible (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run problem] More...
 

Private Member Functions

MoFEMErrorCode setupProblem ()
 [Run problem] More...
 
MoFEMErrorCode createCommonData ()
 [Set up problem] More...
 
MoFEMErrorCode bC ()
 [Create common data] More...
 
MoFEMErrorCode OPs ()
 [Boundary condition] More...
 
MoFEMErrorCode tsSolve ()
 
MoFEMErrorCode checkResults ()
 [Solve] More...
 

Private Attributes

MoFEM::Interface & mField
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 

Detailed Description

Definition at line 203 of file incompressible_elasticity.cpp.

Constructor & Destructor Documentation

◆ Incompressible()

Incompressible::Incompressible ( MoFEM::Interface &  m_field)
inline

Definition at line 205 of file incompressible_elasticity.cpp.

205 : mField(m_field) {}

Member Function Documentation

◆ bC()

MoFEMErrorCode Incompressible::bC ( )
private

[Create common data]

[Boundary condition]

[Natural boundary condition]

[Define gravity vector]

Definition at line 391 of file incompressible_elasticity.cpp.

391  {
392  MoFEMFunctionBegin;
393 
394  auto pipeline_mng = mField.getInterface<PipelineManager>();
395  auto simple = mField.getInterface<Simple>();
396  auto bc_mng = mField.getInterface<BcManager>();
397  auto time_scale = boost::make_shared<TimeScale>();
398 
399  auto integration_rule = [](int, int, int approx_order) {
400  return 2 * (approx_order - 1);
401  };
402 
403  using DomainNaturalBC =
404  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
405  using OpBodyForce =
406  DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
407  using BoundaryNaturalBC =
408  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
409  using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
410 
411  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
412  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
413  pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
414  //! [Natural boundary condition]
415  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
416  pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
417  "FORCE", Sev::inform);
418 
419  //! [Define gravity vector]
420  CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
421  pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
422  "BODY_FORCE", Sev::inform);
423 
424  // Essential BC
425  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
426  "U", 0, 0);
427  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
428  "U", 1, 1);
429  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
430  "U", 2, 2);
431  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
432  simple->getProblemName(), "U");
433 
434  MoFEMFunctionReturn(0);
435 }

◆ checkResults()

MoFEMErrorCode Incompressible::checkResults ( )
private

[Solve]

[Check]

Definition at line 880 of file incompressible_elasticity.cpp.

880 { return 0; }

◆ createCommonData()

MoFEMErrorCode Incompressible::createCommonData ( )
private

[Set up problem]

[Create common data]

Definition at line 361 of file incompressible_elasticity.cpp.

361  {
362  MoFEMFunctionBegin;
363 
364  auto get_options = [&]() {
365  MoFEMFunctionBegin;
366  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
367  &young_modulus, PETSC_NULL);
368  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
369  &poisson_ratio, PETSC_NULL);
370 
371  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
372  << "Young modulus " << young_modulus;
373  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
374  << "Poisson_ratio " << poisson_ratio;
375 
376  mu = young_modulus / (2. + 2. * poisson_ratio);
377  const double lambda_denom =
378  (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
379  lambda = young_modulus * poisson_ratio / lambda_denom;
380 
381  MoFEMFunctionReturn(0);
382  };
383 
384  CHKERR get_options();
385 
386  MoFEMFunctionReturn(0);
387 }

◆ OPs()

MoFEMErrorCode Incompressible::OPs ( )
private

[Boundary condition]

[Push operators to pip]

[Operators used for RHS incompressible elasticity]

[Operators used for incompressible elasticity]

[Operators used for incompressible elasticity]

[Operators used for RHS incompressible elasticity]

[Operators used for RHS incompressible elasticity]

Definition at line 439 of file incompressible_elasticity.cpp.

439  {
440  MoFEMFunctionBegin;
441  auto simple = mField.getInterface<Simple>();
442  auto pip_mng = mField.getInterface<PipelineManager>();
443  auto bc_mng = mField.getInterface<BcManager>();
444 
445  auto integration_rule_vol = [](int, int, int approx_order) {
446  return 2 * approx_order + geom_order - 1;
447  };
448 
449  auto add_domain_base_ops = [&](auto &pip) {
450  MoFEMFunctionBegin;
451  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, L2},
452  "GEOMETRY");
453  MoFEMFunctionReturn(0);
454  };
455 
456  auto add_domain_ops_lhs = [&](auto &pip) {
457  MoFEMFunctionBegin;
458 
459  // This assemble A-matrix
460  using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
461  AT>::BiLinearForm<GAUSS>::OpMass<1, 1>;
462  // This assemble B-matrix (preconditioned)
463  using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
464  AT>::BiLinearForm<GAUSS>::OpMass<1, 1>;
465  //! [Operators used for RHS incompressible elasticity]
466 
467  //! [Operators used for incompressible elasticity]
468  using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
469  AT>::BiLinearForm<IT>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
470  using OpMixScalarTimesDiv = FormsIntegrators<DomainEleOp>::Assembly<
471  AT>::BiLinearForm<IT>::OpMixScalarTimesDiv<SPACE_DIM, coord_type>;
472  //! [Operators used for incompressible elasticity]
473 
474  auto mat_D_ptr = boost::make_shared<MatrixDouble>();
475  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
476  mat_D_ptr->resize(size_symm * size_symm, 1);
477 
478  FTensor::Index<'i', SPACE_DIM> i;
479  FTensor::Index<'j', SPACE_DIM> j;
480  FTensor::Index<'k', SPACE_DIM> k;
481  FTensor::Index<'l', SPACE_DIM> l;
482 
483  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
484  auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
485  t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
486 
487  pip.push_back(new OpMixScalarTimesDiv(
488  "P", "U",
489  [](const double, const double, const double) constexpr { return -1.; },
490  true, false));
491  pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
492 
493  auto get_lambda_reciprocal = [](const double, const double, const double) {
494  return 1. / lambda;
495  };
496  if (poisson_ratio < 0.5)
497  pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
498 
499  double eps_stab = 0;
500  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
501  PETSC_NULL);
502  if (eps_stab > 0)
503  pip.push_back(new OpMassPressureStab(
504  "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
505 
506  MoFEMFunctionReturn(0);
507  };
508 
509  auto add_domain_ops_rhs = [&](auto &pip) {
510  MoFEMFunctionBegin;
511 
512  //! [Operators used for RHS incompressible elasticity]
513  using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
514  AT>::LinearForm<GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
515 
516  using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
517  AT>::LinearForm<GAUSS>::OpMixDivTimesU<1, SPACE_DIM, SPACE_DIM>;
518 
519  using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
520  AT>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
521 
522  //! [Operators used for RHS incompressible elasticity]
523 
524  auto pressure_ptr = boost::make_shared<VectorDouble>();
525  pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
526 
527  auto div_u_ptr = boost::make_shared<VectorDouble>();
528  pip.push_back(
529  new OpCalculateDivergenceVectorFieldValues<SPACE_DIM>("U", div_u_ptr));
530 
531  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
532  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
533  "U", grad_u_ptr));
534 
535  auto strain_ptr = boost::make_shared<MatrixDouble>();
536  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
537 
538  auto get_four_mu = [](const double, const double, const double) {
539  return -2. * mu;
540  };
541  auto minus_one = [](const double, const double, const double) constexpr {
542  return -1.;
543  };
544 
545  pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
546 
547  pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
548 
549  pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
550 
551  auto get_lambda_reciprocal = [](const double, const double, const double) {
552  return 1. / lambda;
553  };
554  if (poisson_ratio < 0.5) {
555  pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
556  get_lambda_reciprocal));
557  }
558 
559  MoFEMFunctionReturn(0);
560  };
561 
562  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
563  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
564  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
565  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
566 
567  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
568  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
569 
570  MoFEMFunctionReturn(0);
571 }

◆ runProblem()

MoFEMErrorCode Incompressible::runProblem ( )

[Run problem]

Definition at line 273 of file incompressible_elasticity.cpp.

273  {
274  MoFEMFunctionBegin;
275  CHKERR setupProblem();
276  CHKERR createCommonData();
277  CHKERR bC();
278  CHKERR OPs();
279  CHKERR tsSolve();
280  CHKERR checkResults();
281  MoFEMFunctionReturn(0);
282 }

◆ setupProblem()

MoFEMErrorCode Incompressible::setupProblem ( )
private

[Run problem]

[Set up problem]

Definition at line 286 of file incompressible_elasticity.cpp.

286  {
287  MoFEMFunctionBegin;
288  Simple *simple = mField.getInterface<Simple>();
289 
290  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
291  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
292  PETSC_NULL);
293  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_discontinuous_pressure",
294  &isDiscontinuousPressure, PETSC_NULL);
295 
296  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
297  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
298 
299  // Select base
300  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
301  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
302  PetscInt choice_base_value = AINSWORTH;
303  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
304  LASBASETOPT, &choice_base_value, PETSC_NULL);
305 
306  FieldApproximationBase base;
307  switch (choice_base_value) {
308  case AINSWORTH:
309  base = AINSWORTH_LEGENDRE_BASE;
310  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
311  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
312  break;
313  case DEMKOWICZ:
314  base = DEMKOWICZ_JACOBI_BASE;
315  MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
316  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
317  break;
318  default:
319  base = LASTBASE;
320  break;
321  }
322 
323  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
324  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
325  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
326  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
327 
328  // Adding fields related to incompressible elasticity
329  // Add displacement domain and boundary fields
330  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
331  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
332  CHKERR simple->setFieldOrder("U", order);
333 
334  // Add pressure domain and boundary fields
335  // Choose either Crouzeix-Raviart element:
336  if (isDiscontinuousPressure) {
337  CHKERR simple->addDomainField("P", L2, base, 1);
338  CHKERR simple->setFieldOrder("P", order - 2);
339  } else {
340  // ... or Taylor-Hood element:
341  CHKERR simple->addDomainField("P", H1, base, 1);
342  CHKERR simple->setFieldOrder("P", order - 1);
343  }
344 
345  // Add geometry data field
346  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
347  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
348 
349  CHKERR simple->setUp();
350 
351  auto project_ho_geometry = [&]() {
352  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
353  return mField.loop_dofs("GEOMETRY", ent_method);
354  };
355  CHKERR project_ho_geometry();
356 
357  MoFEMFunctionReturn(0);
358 } //! [Set up problem]

◆ tsSolve()

MoFEMErrorCode Incompressible::tsSolve ( )
private

Definition at line 584 of file incompressible_elasticity.cpp.

584  {
585  MoFEMFunctionBegin;
586 
587  auto simple = mField.getInterface<Simple>();
588  auto pip_mng = mField.getInterface<PipelineManager>();
589  auto is_manager = mField.getInterface<ISManager>();
590 
591  auto set_section_monitor = [&](auto solver) {
592  MoFEMFunctionBegin;
593  SNES snes;
594  CHKERR TSGetSNES(solver, &snes);
595  PetscViewerAndFormat *vf;
596  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
597  PETSC_VIEWER_DEFAULT, &vf);
598  CHKERR SNESMonitorSet(
599  snes,
600  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
601  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
602  MoFEMFunctionReturn(0);
603  };
604 
605  auto scatter_create = [&](auto D, auto coeff) {
606  SmartPetscObj<IS> is;
607  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
608  ROW, "U", coeff, coeff, is);
609  int loc_size;
610  CHKERR ISGetLocalSize(is, &loc_size);
611  Vec v;
612  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
613  VecScatter scatter;
614  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
615  return std::make_tuple(SmartPetscObj<Vec>(v),
616  SmartPetscObj<VecScatter>(scatter));
617  };
618 
619  auto create_post_process_elements = [&]() {
620  auto pp_fe = boost::make_shared<PostProcEle>(mField);
621  auto &pip = pp_fe->getOpPtrVector();
622 
623  auto push_vol_ops = [this](auto &pip) {
624  MoFEMFunctionBegin;
625  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
626  "GEOMETRY");
627 
628  MoFEMFunctionReturn(0);
629  };
630 
631  auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
632  MoFEMFunctionBegin;
633 
634  auto &pip = pp_fe->getOpPtrVector();
635 
636  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
637 
638  auto x_ptr = boost::make_shared<MatrixDouble>();
639  pip.push_back(
640  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
641  auto u_ptr = boost::make_shared<MatrixDouble>();
642  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
643 
644  auto pressure_ptr = boost::make_shared<VectorDouble>();
645  pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
646 
647  auto div_u_ptr = boost::make_shared<VectorDouble>();
648  pip.push_back(new OpCalculateDivergenceVectorFieldValues<SPACE_DIM>(
649  "U", div_u_ptr));
650 
651  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
652  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
653  "U", grad_u_ptr));
654 
655  auto strain_ptr = boost::make_shared<MatrixDouble>();
656  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
657 
658  auto stress_ptr = boost::make_shared<MatrixDouble>();
659  pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
660  mu, stress_ptr, strain_ptr, pressure_ptr));
661 
662  pip.push_back(
663 
664  new OpPPMap(
665 
666  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
667 
668  {{"P", pressure_ptr}},
669 
670  {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
671 
672  {},
673 
674  {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
675 
676  )
677 
678  );
679 
680  MoFEMFunctionReturn(0);
681  };
682 
683  auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
684  PetscBool post_proc_vol = PETSC_FALSE;
685  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
686  &post_proc_vol, PETSC_NULL);
687  if (post_proc_vol == PETSC_FALSE)
688  return boost::shared_ptr<PostProcEle>();
689  auto pp_fe = boost::make_shared<PostProcEle>(mField);
690  CHK_MOAB_THROW(
691  push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
692  "push_vol_post_proc_ops");
693  return pp_fe;
694  };
695 
696  auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
697  PetscBool post_proc_skin = PETSC_TRUE;
698  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
699  &post_proc_skin, PETSC_NULL);
700  if (post_proc_skin == PETSC_FALSE)
701  return boost::shared_ptr<SkinPostProcEle>();
702 
703  auto simple = mField.getInterface<Simple>();
704  auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
705  auto op_side = new OpLoopSide<DomainEle>(
706  mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
707  pp_fe->getOpPtrVector().push_back(op_side);
708  CHK_MOAB_THROW(push_vol_post_proc_ops(
709  pp_fe, push_vol_ops(op_side->getOpPtrVector())),
710  "push_vol_post_proc_ops");
711  return pp_fe;
712  };
713 
714  return std::make_pair(vol_post_proc(), skin_post_proc());
715  };
716 
717  boost::shared_ptr<SetPtsData> field_eval_data;
718  boost::shared_ptr<MatrixDouble> vector_field_ptr;
719 
720  std::array<double, SPACE_DIM> field_eval_coords;
721  int coords_dim = SPACE_DIM;
722  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
723  field_eval_coords.data(), &coords_dim,
724  &doEvalField);
725 
726  if (doEvalField) {
727  vector_field_ptr = boost::make_shared<MatrixDouble>();
728  field_eval_data =
729  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
730  if constexpr (SPACE_DIM == 3) {
731  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
732  field_eval_data, simple->getDomainFEName());
733  } else {
734  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
735  field_eval_data, simple->getDomainFEName());
736  }
737 
738  field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
739  auto no_rule = [](int, int, int) { return -1; };
740  auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
741  field_eval_fe_ptr->getRuleHook = no_rule;
742  field_eval_fe_ptr->getOpPtrVector().push_back(
743  new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
744  }
745 
746  auto set_time_monitor = [&](auto dm, auto solver) {
747  MoFEMFunctionBegin;
748  boost::shared_ptr<MonitorIncompressible> monitor_ptr(
749  new MonitorIncompressible(SmartPetscObj<DM>(dm, true),
750  create_post_process_elements(), uXScatter,
751  uYScatter, uZScatter, field_eval_coords,
752  field_eval_data, vector_field_ptr));
753  boost::shared_ptr<ForcesAndSourcesCore> null;
754  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
755  monitor_ptr, null, null);
756  MoFEMFunctionReturn(0);
757  };
758 
759  auto set_essential_bc = [&]() {
760  MoFEMFunctionBegin;
761  // This is low level pushing finite elements (pipelines) to solver
762  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
763  auto pre_proc_ptr = boost::make_shared<FEMethod>();
764  auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
765 
766  // Add boundary condition scaling
767  auto time_scale = boost::make_shared<TimeScale>();
768 
769  pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
770  mField, pre_proc_ptr, {time_scale}, false);
771  post_proc_rhs_ptr->postProcessHook =
772  EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs_ptr,
773  1.);
774  ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
775  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
776  ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
777  MoFEMFunctionReturn(0);
778  };
779 
780  auto set_schur_pc = [&](auto solver) {
781  SNES snes;
782  CHKERR TSGetSNES(solver, &snes);
783  KSP ksp;
784  CHKERR SNESGetKSP(snes, &ksp);
785  PC pc;
786  CHKERR KSPGetPC(ksp, &pc);
787  PetscBool is_pcfs = PETSC_FALSE;
788  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
789  boost::shared_ptr<SetUpSchur> schur_ptr;
790  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
791 
792  auto A = createDMMatrix(simple->getDM());
793  auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
794 
795  CHK_MOAB_THROW(
796  TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
797  "set operators");
798  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
799  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
800  pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
801  MoFEMFunctionBegin;
802  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
803  CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
804  CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
805  MoFEMFunctionReturn(0);
806  };
807  post_proc_schur_lhs_ptr->postProcessHook = [this,
808  post_proc_schur_lhs_ptr]() {
809  MoFEMFunctionBegin;
810  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
811  *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
812  CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
813  CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
814  CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
815  mField, post_proc_schur_lhs_ptr, 1.,
816  SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
817 
818  CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
819  CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
820  CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
821  SAME_NONZERO_PATTERN);
822  MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
823  MoFEMFunctionReturn(0);
824  };
825  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
826  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
827 
828  if (is_pcfs == PETSC_TRUE) {
829  if (AT == AssemblyType::SCHUR) {
830  schur_ptr = SetUpSchur::createSetUpSchur(mField);
831  CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
832  } else {
833  auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
834  MoFEMFunctionBegin;
835  auto bc_mng = mField.getInterface<BcManager>();
836  auto name_prb = simple->getProblemName();
837  SmartPetscObj<IS> is_p;
838  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
839  name_prb, ROW, "P", 0, 1, is_p);
840  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_p);
841  MoFEMFunctionReturn(0);
842  };
843  CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
844  "set pcfieldsplit preconditioned");
845  }
846  return boost::make_tuple(schur_ptr, A, B);
847  }
848 
849  return boost::make_tuple(schur_ptr, A, B);
850  };
851 
852  auto dm = simple->getDM();
853  auto D = createDMVector(dm);
854 
855  uXScatter = scatter_create(D, 0);
856  uYScatter = scatter_create(D, 1);
857  if (SPACE_DIM == 3)
858  uZScatter = scatter_create(D, 2);
859 
860  // Add extra finite elements to SNES solver pipelines to resolve essential
861  // boundary conditions
862  CHKERR set_essential_bc();
863 
864  auto solver = pip_mng->createTSIM();
865  CHKERR TSSetFromOptions(solver);
866 
867  CHKERR set_section_monitor(solver);
868  CHKERR set_time_monitor(dm, solver);
869  auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
870 
871  CHKERR TSSetSolution(solver, D);
872  CHKERR TSSetUp(solver);
873  CHKERR TSSolve(solver, NULL);
874 
875  MoFEMFunctionReturn(0);
876 }

Member Data Documentation

◆ mField

MoFEM::Interface& Incompressible::mField
private

Definition at line 210 of file incompressible_elasticity.cpp.

◆ uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Incompressible::uXScatter
private

Definition at line 219 of file incompressible_elasticity.cpp.

◆ uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Incompressible::uYScatter
private

Definition at line 220 of file incompressible_elasticity.cpp.

◆ uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Incompressible::uZScatter
private

Definition at line 221 of file incompressible_elasticity.cpp.

The documentation for this struct was generated from the following file:
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
Incompressible::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: incompressible_elasticity.cpp:219
young_modulus
static double young_modulus
Definition: incompressible_elasticity.cpp:196
MonitorIncompressible
Definition: incompressible_elasticity.cpp:42
poisson_ratio
static double poisson_ratio
Definition: incompressible_elasticity.cpp:197
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
H1
@ H1
continuous field
Definition: definitions.h:85
LASTBASE
@ LASTBASE
Definition: definitions.h:69
Incompressible::bC
MoFEMErrorCode bC()
[Create common data]
Definition: incompressible_elasticity.cpp:391
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:486
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
isDiscontinuousPressure
PetscBool isDiscontinuousPressure
Definition: incompressible_elasticity.cpp:201
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
mu
static double mu
Definition: incompressible_elasticity.cpp:198
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
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::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:125
ROW
@ ROW
Definition: definitions.h:136
Incompressible::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: incompressible_elasticity.cpp:439
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
OpCalculateLameStress
Definition: incompressible_elasticity.cpp:225
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
geom_order
int geom_order
Definition: incompressible_elasticity.cpp:195
Incompressible::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: incompressible_elasticity.cpp:221
Incompressible::mField
MoFEM::Interface & mField
Definition: incompressible_elasticity.cpp:210
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
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
double
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
SPACE_DIM
constexpr int SPACE_DIM
Definition: incompressible_elasticity.cpp:10
Incompressible::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: incompressible_elasticity.cpp:361
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
Incompressible::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: incompressible_elasticity.cpp:286
Incompressible::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: incompressible_elasticity.cpp:220
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:234
order
int order
[Specialisation for assembly]
Definition: incompressible_elasticity.cpp:194
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
FTensor::Index< 'i', SPACE_DIM >
Incompressible::tsSolve
MoFEMErrorCode tsSolve()
Definition: incompressible_elasticity.cpp:584
AT
constexpr AssemblyType AT
Definition: incompressible_elasticity.cpp:13
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::OpSymmetrizeTensor
Definition: UserDataOperators.hpp:1948
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
approx_order
int approx_order
Definition: test_broken_space.cpp:50
coord_type
constexpr CoordinateTypes coord_type
Definition: incompressible_elasticity.cpp:19
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:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::SmartPetscObj< IS >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:768
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
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:429
Incompressible::checkResults
MoFEMErrorCode checkResults()
[Solve]
Definition: incompressible_elasticity.cpp:880
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow