v0.15.0
Loading...
Searching...
No Matches
Incompressible Struct Reference
Collaboration diagram for Incompressible:
[legend]

Public Member Functions

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

Private Member Functions

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

Private Attributes

MoFEM::InterfacemField
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
std::array< double, 3 > fieldEvalCoords
 

Detailed Description

Definition at line 185 of file incompressible_elasticity.cpp.

Constructor & Destructor Documentation

◆ Incompressible()

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

Definition at line 187 of file incompressible_elasticity.cpp.

187: 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 378 of file incompressible_elasticity.cpp.

378 {
380
381 auto pipeline_mng = mField.getInterface<PipelineManager>();
383 auto bc_mng = mField.getInterface<BcManager>();
384 auto time_scale = boost::make_shared<TimeScale>();
385
386 auto integration_rule = [](int, int, int approx_order) {
387 return 2 * (approx_order - 1);
388 };
389
390 using DomainNaturalBC =
392 using OpBodyForce =
394 using BoundaryNaturalBC =
397
398 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
400 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
401 //! [Natural boundary condition]
403 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
404 "FORCE", Sev::inform);
405
406 //! [Define gravity vector]
408 pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
409 "BODY_FORCE", Sev::inform);
410
411 // Essential BC
412 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
413 "U", 0, 0);
414 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
415 "U", 1, 1);
416 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
417 "U", 2, 2);
418 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
419 simple->getProblemName(), "U");
420
422}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto integration_rule
constexpr int SPACE_DIM
static constexpr int approx_order
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Assembly methods.
Definition Natural.hpp:65
PipelineManager interface.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

◆ checkResults()

MoFEMErrorCode Incompressible::checkResults ( )
private

[Solve]

[Check]

no elements to evaluate

Definition at line 853 of file incompressible_elasticity.cpp.

853 {
855
856 int atom_test = 0;
857 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
858 PETSC_NULLPTR);
859
860 if (atom_test) {
861 double eps = 1e-6;
862 double Ux_ref, Uy_ref, Uz_ref;
863 string test_name;
864
865 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
866
867 auto field_eval_data =
868 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
869 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
870 field_eval_data, mField.getInterface<Simple>()->getDomainFEName());
871
872 field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
873 auto no_rule = [](int, int, int) { return -1; };
874 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
875 field_eval_fe_ptr->getRuleHook = no_rule;
876 field_eval_fe_ptr->getOpPtrVector().push_back(
877 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
878
880 ->evalFEAtThePoint<SPACE_DIM>(
881 fieldEvalCoords.data(), 1e-12,
883 mField.getInterface<Simple>()->getDomainFEName(), field_eval_data,
885 QUIET);
886
887 auto eval_num_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
888 CHKERR VecZeroEntries(eval_num_vec);
889 if (vector_field_ptr->size1()) {
890 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
891 }
892 CHKERR VecAssemblyBegin(eval_num_vec);
893 CHKERR VecAssemblyEnd(eval_num_vec);
894
895 double eval_num;
896 CHKERR VecSum(eval_num_vec, &eval_num);
897 if (!(int)eval_num) { /// no elements to evaluate
898 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
899 "atom test %d failed: did not find elements to evaluate "
900 "the field, check the coordinates",
901 atom_test);
902 }
903
904 if (vector_field_ptr->size1()) {
905 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vector_field_ptr);
906 switch (atom_test) {
907 case 1: // Cooke's Membrane
908 Ux_ref = -5.89757;
909 Uy_ref = 8.15644;
910 Uz_ref = 0.0;
911 test_name = "Cooke's Membrane";
912 break;
913 default:
914 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
915 "atom test %d does not exist", atom_test);
916 }
917
918 auto calculate_error = [](double computed, double reference,
919 double tolerance) {
920 return fabs(computed - reference) > tolerance * fabs(reference);
921 };
922
923 if (calculate_error(t_disp(0), Ux_ref, eps) ||
924 calculate_error(t_disp(1), Uy_ref, eps) ||
925 (SPACE_DIM == 3 && calculate_error(t_disp(2), Uz_ref, eps))) {
926 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
927 "atom test %d:%s failed: Displacement exceeds tolerance: Ux = "
928 "%f, Uy = %f",
929 atom_test, test_name.c_str(), t_disp(0), t_disp(1));
930 }
931 }
932 }
933
935}
static const double eps
@ QUIET
@ MF_EXIST
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_INVALID_DATA
Definition definitions.h:36
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
int atom_test
Atom test.
Definition plastic.cpp:121
std::array< double, 3 > fieldEvalCoords
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Field evaluator interface.
Get values at integration pts for tensor field rank 1, i.e. vector field.
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:410
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:389

◆ createCommonData()

MoFEMErrorCode Incompressible::createCommonData ( )
private

[Set up problem]

[Create common data]

Definition at line 348 of file incompressible_elasticity.cpp.

348 {
350
351 auto get_options = [&]() {
353 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
354 &young_modulus, PETSC_NULLPTR);
355 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
356 &poisson_ratio, PETSC_NULLPTR);
357
358 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
359 << "Young modulus " << young_modulus;
360 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
361 << "Poisson_ratio " << poisson_ratio;
362
363 mu = young_modulus / (2. + 2. * poisson_ratio);
364 const double lambda_denom =
365 (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
366 lambda = young_modulus * poisson_ratio / lambda_denom;
367
369 };
370
371 CHKERR get_options();
372
374}
#define MOFEM_LOG(channel, severity)
Log.
static double young_modulus
static double lambda
static double poisson_ratio
static double mu
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)

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

426 {
428 auto pip_mng = mField.getInterface<PipelineManager>();
429
430 auto integration_rule_vol = [](int, int, int approx_order) {
431 return 2 * approx_order + geom_order - 1;
432 };
433
434 auto add_domain_base_ops = [&](auto &pip) {
437 "GEOMETRY");
439 };
440
441 auto add_domain_ops_lhs = [&](auto &pip) {
443
444 // This assemble A-matrix
445 using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
447 // This assemble B-matrix (preconditioned)
448 using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
450 //! [Operators used for RHS incompressible elasticity]
451
452 //! [Operators used for incompressible elasticity]
453 using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
457 //! [Operators used for incompressible elasticity]
458
459 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
460 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
461 mat_D_ptr->resize(size_symm * size_symm, 1);
462
467
469 auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
470 t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
471
472 pip.push_back(new OpMixScalarTimesDiv(
473 "P", "U",
474 [](const double, const double, const double) constexpr { return -1.; },
475 true, false));
476 pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
477
478 auto get_lambda_reciprocal = [](const double, const double, const double) {
479 return 1. / lambda;
480 };
481 if (poisson_ratio < 0.5)
482 pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
483
484 double eps_stab = 0;
485 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
486 PETSC_NULLPTR);
487 if (eps_stab > 0)
488 pip.push_back(new OpMassPressureStab(
489 "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
490
492 };
493
494 auto add_domain_ops_rhs = [&](auto &pip) {
496
497 //! [Operators used for RHS incompressible elasticity]
498 using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
500
501 using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
503
504 using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
506
507 //! [Operators used for RHS incompressible elasticity]
508
509 auto pressure_ptr = boost::make_shared<VectorDouble>();
510 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
511
512 auto div_u_ptr = boost::make_shared<VectorDouble>();
513 pip.push_back(
515
516 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
518 "U", grad_u_ptr));
519
520 auto strain_ptr = boost::make_shared<MatrixDouble>();
521 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
522
523 auto get_four_mu = [](const double, const double, const double) {
524 return -2. * mu;
525 };
526 auto minus_one = [](const double, const double, const double) constexpr {
527 return -1.;
528 };
529
530 pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
531
532 pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
533
534 pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
535
536 auto get_lambda_reciprocal = [](const double, const double, const double) {
537 return 1. / lambda;
538 };
539 if (poisson_ratio < 0.5) {
540 pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
541 get_lambda_reciprocal));
542 }
543
545 };
546
547 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
548 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
549 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
550 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
551
552 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
553 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
554
556}
Kronecker Delta class symmetric.
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
FTensor::Index< 'i', SPACE_DIM > i
constexpr AssemblyType AT
constexpr CoordinateTypes coord_type
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
constexpr auto size_symm
Definition plastic.cpp:42
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.

◆ runProblem()

MoFEMErrorCode Incompressible::runProblem ( )

[Run problem]

Definition at line 256 of file incompressible_elasticity.cpp.

256 {
260 CHKERR bC();
261 CHKERR OPs();
262 CHKERR tsSolve();
265}
MoFEMErrorCode checkResults()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode bC()
[Create common data]

◆ setupProblem()

MoFEMErrorCode Incompressible::setupProblem ( )
private

[Run problem]

[Set up problem]

Definition at line 269 of file incompressible_elasticity.cpp.

269 {
272
273 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
274 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
275 PETSC_NULLPTR);
276 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_discontinuous_pressure",
277 &isDiscontinuousPressure, PETSC_NULLPTR);
278
279 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
280 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
281
282 // Select base
283 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
284 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
285 PetscInt choice_base_value = AINSWORTH;
286 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, PETSC_NULLPTR, "-base", list_bases,
287 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
288 int coords_dim = 3;
289 CHKERR PetscOptionsGetRealArray(PETSC_NULLPTR, PETSC_NULLPTR,
290 "-field_eval_coords", fieldEvalCoords.data(),
291 &coords_dim, &doEvalField);
292
294 switch (choice_base_value) {
295 case AINSWORTH:
297 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
298 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
299 break;
300 case DEMKOWICZ:
302 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
303 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
304 break;
305 default:
306 base = LASTBASE;
307 break;
308 }
309
310 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
311 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
312 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
313 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
314
315 // Adding fields related to incompressible elasticity
316 // Add displacement domain and boundary fields
317 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
318 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
319 CHKERR simple->setFieldOrder("U", order);
320
321 // Add pressure domain and boundary fields
322 // Choose either Crouzeix-Raviart element:
324 CHKERR simple->addDomainField("P", L2, base, 1);
325 CHKERR simple->setFieldOrder("P", order - 2);
326 } else {
327 // ... or Taylor-Hood element:
328 CHKERR simple->addDomainField("P", H1, base, 1);
329 CHKERR simple->setFieldOrder("P", order - 1);
330 }
331
332 // Add geometry data field
333 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
334 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
335
336 CHKERR simple->setUp();
337
338 auto project_ho_geometry = [&]() {
339 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
340 return mField.loop_dofs("GEOMETRY", ent_method);
341 };
342 CHKERR project_ho_geometry();
343
345} //! [Set up problem]
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
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.
PetscBool isDiscontinuousPressure
int order
[Specialisation for assembly]
PetscBool doEvalField
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Projection of edge entities with one mid-node on hierarchical basis.

◆ tsSolve()

MoFEMErrorCode Incompressible::tsSolve ( )
private

Definition at line 569 of file incompressible_elasticity.cpp.

569 {
571
573 auto pip_mng = mField.getInterface<PipelineManager>();
574 auto is_manager = mField.getInterface<ISManager>();
575
576 auto set_section_monitor = [&](auto solver) {
578 SNES snes;
579 CHKERR TSGetSNES(solver, &snes);
580 PetscViewerAndFormat *vf;
581 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
582 PETSC_VIEWER_DEFAULT, &vf);
583 CHKERR SNESMonitorSet(
584 snes,
585 (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal,
586 void *))SNESMonitorFields,
587 vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
589 };
590
591 auto scatter_create = [&](auto D, auto coeff) {
593 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
594 ROW, "U", coeff, coeff, is);
595 int loc_size;
596 CHKERR ISGetLocalSize(is, &loc_size);
597 Vec v;
598 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
599 VecScatter scatter;
600 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
601 return std::make_tuple(SmartPetscObj<Vec>(v),
603 };
604
605 auto create_post_process_elements = [&]() {
606 auto pp_fe = boost::make_shared<PostProcEle>(mField);
607
608 auto push_vol_ops = [this](auto &pip) {
611 "GEOMETRY");
612
614 };
615
616 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
618
619 auto &pip = pp_fe->getOpPtrVector();
620
622
623 auto x_ptr = boost::make_shared<MatrixDouble>();
624 pip.push_back(
625 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
626 auto u_ptr = boost::make_shared<MatrixDouble>();
627 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
628
629 auto pressure_ptr = boost::make_shared<VectorDouble>();
630 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
631
632 auto div_u_ptr = boost::make_shared<VectorDouble>();
634 "U", div_u_ptr));
635
636 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
638 "U", grad_u_ptr));
639
640 auto strain_ptr = boost::make_shared<MatrixDouble>();
641 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
642
643 auto stress_ptr = boost::make_shared<MatrixDouble>();
644 pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
645 mu, stress_ptr, strain_ptr, pressure_ptr));
646
647 pip.push_back(
648
649 new OpPPMap(
650
651 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
652
653 {{"P", pressure_ptr}},
654
655 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
656
657 {},
658
659 {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
660
661 )
662
663 );
664
666 };
667
668 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
669 PetscBool post_proc_vol = PETSC_FALSE;
670 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
671 &post_proc_vol, PETSC_NULLPTR);
672 if (post_proc_vol == PETSC_FALSE)
673 return boost::shared_ptr<PostProcEle>();
674 auto pp_fe = boost::make_shared<PostProcEle>(mField);
676 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
677 "push_vol_post_proc_ops");
678 return pp_fe;
679 };
680
681 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
682 PetscBool post_proc_skin = PETSC_TRUE;
683 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
684 &post_proc_skin, PETSC_NULLPTR);
685 if (post_proc_skin == PETSC_FALSE)
686 return boost::shared_ptr<SkinPostProcEle>();
687
689 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
690 auto op_side = new OpLoopSide<DomainEle>(
691 mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
692 pp_fe->getOpPtrVector().push_back(op_side);
693 CHK_MOAB_THROW(push_vol_post_proc_ops(
694 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
695 "push_vol_post_proc_ops");
696 return pp_fe;
697 };
698
699 return std::make_pair(vol_post_proc(), skin_post_proc());
700 };
701
702 boost::shared_ptr<SetPtsData> field_eval_data;
703 boost::shared_ptr<MatrixDouble> vector_field_ptr;
704
705 if (doEvalField) {
706 vector_field_ptr = boost::make_shared<MatrixDouble>();
707 field_eval_data =
708 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
709 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
710 field_eval_data, simple->getDomainFEName());
711
712 field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
713 auto no_rule = [](int, int, int) { return -1; };
714 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
715 field_eval_fe_ptr->getRuleHook = no_rule;
716 field_eval_fe_ptr->getOpPtrVector().push_back(
717 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
718 }
719
720 auto set_time_monitor = [&](auto dm, auto solver) {
722 boost::shared_ptr<MonitorIncompressible> monitor_ptr(
724 create_post_process_elements(), uXScatter,
726 field_eval_data, vector_field_ptr));
727 boost::shared_ptr<ForcesAndSourcesCore> null;
728 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
729 monitor_ptr, null, null);
731 };
732
733 auto set_essential_bc = [&]() {
735 // This is low level pushing finite elements (pipelines) to solver
736 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
737 auto pre_proc_ptr = boost::make_shared<FEMethod>();
738 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
739
740 // Add boundary condition scaling
741 auto time_scale = boost::make_shared<TimeScale>();
742
743 pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
744 mField, pre_proc_ptr, {time_scale}, false);
745 post_proc_rhs_ptr->postProcessHook =
747 1.);
748 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
749 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
750 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
752 };
753
754 auto set_schur_pc = [&](auto solver) {
755 SNES snes;
756 CHKERR TSGetSNES(solver, &snes);
757 KSP ksp;
758 CHKERR SNESGetKSP(snes, &ksp);
759 PC pc;
760 CHKERR KSPGetPC(ksp, &pc);
761 PetscBool is_pcfs = PETSC_FALSE;
762 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
763 boost::shared_ptr<SetUpSchur> schur_ptr;
764 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
765
766 auto A = createDMMatrix(simple->getDM());
767 auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
768
770 TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
771 "set operators");
772 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
773 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
774 pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
776 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
777 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
778 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
780 };
781 post_proc_schur_lhs_ptr->postProcessHook = [this,
782 post_proc_schur_lhs_ptr]() {
784 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
785 *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
786 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
787 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
789 mField, post_proc_schur_lhs_ptr, 1.,
790 SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
791
792 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
793 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
794 CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
795 SAME_NONZERO_PATTERN);
796 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
798 };
799 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
800 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
801
802 if (is_pcfs == PETSC_TRUE) {
803 if (AT == AssemblyType::SCHUR) {
805 CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
806 } else {
807 auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
809 auto name_prb = simple->getProblemName();
811 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
812 name_prb, ROW, "P", 0, 1, is_p);
813 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_p);
815 };
816 CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
817 "set pcfieldsplit preconditioned");
818 }
819 return boost::make_tuple(schur_ptr, A, B);
820 }
821
822 return boost::make_tuple(schur_ptr, A, B);
823 };
824
825 auto dm = simple->getDM();
826 auto D = createDMVector(dm);
827
828 uXScatter = scatter_create(D, 0);
829 uYScatter = scatter_create(D, 1);
830 if (SPACE_DIM == 3)
831 uZScatter = scatter_create(D, 2);
832
833 // Add extra finite elements to SNES solver pipelines to resolve essential
834 // boundary conditions
835 CHKERR set_essential_bc();
836
837 auto solver = pip_mng->createTSIM();
838 CHKERR TSSetFromOptions(solver);
839
840 CHKERR set_section_monitor(solver);
841 CHKERR set_time_monitor(dm, solver);
842 auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
843
844 CHKERR TSSetSolution(solver, D);
845 CHKERR TSSetUp(solver);
846 CHKERR TSSolve(solver, NULL);
847
849}
@ ROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
double D
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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:1046
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
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
intrusive_ptr for managing petsc objects
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)

Member Data Documentation

◆ fieldEvalCoords

std::array<double, 3> Incompressible::fieldEvalCoords
private

Definition at line 204 of file incompressible_elasticity.cpp.

◆ mField

MoFEM::Interface& Incompressible::mField
private

Definition at line 192 of file incompressible_elasticity.cpp.

◆ uXScatter

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

Definition at line 201 of file incompressible_elasticity.cpp.

◆ uYScatter

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

Definition at line 202 of file incompressible_elasticity.cpp.

◆ uZScatter

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

Definition at line 203 of file incompressible_elasticity.cpp.


The documentation for this struct was generated from the following file: