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 876 of file incompressible_elasticity.cpp.

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

◆ 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 582 of file incompressible_elasticity.cpp.

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

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:488
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:468
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:127
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:169
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:139
BiLinearForm
FTensor::Index< 'i', SPACE_DIM >
Incompressible::tsSolve
MoFEMErrorCode tsSolve()
Definition: incompressible_elasticity.cpp:582
AT
constexpr AssemblyType AT
Definition: incompressible_elasticity.cpp:13
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
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:1974
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
approx_order
int approx_order
Definition: test_broken_space.cpp:54
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:64
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:802
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:876
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:1291
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:709
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