390 {
391
392 const string default_options = "-ksp_type fgmres \n"
393 "-pc_type lu \n"
394 "-pc_factor_mat_solver_type mumps \n"
395 "-mat_mumps_icntl_20 0 \n"
396 "-ksp_monitor\n";
397
398 string param_file = "param_file.petsc";
399 if (!static_cast<bool>(ifstream(param_file))) {
400 std::ofstream file(param_file.c_str(), std::ios::ate);
401 if (file.is_open()) {
402 file << default_options;
403 file.close();
404 }
405 }
406
408
409
410 moab::Core mb_instance;
411 moab::Interface &moab = mb_instance;
412
413 try {
414
417
419
420
422 PetscBool flg_test = PETSC_FALSE;
423 PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
424 "none");
425
426
428 PETSC_NULLPTR);
429
430
431 CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
432 &flg_test, PETSC_NULLPTR);
433
434 PetscOptionsEnd();
435
437
440
441 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
442
443 enum MyBcTypes {
444 FIX_BRICK_FACES = 1,
445 FIX_NODES = 2,
446 BRICK_PRESSURE_FACES = 3,
447 FIX_NODES_Y_DIR = 4
448 };
449
452 int id =
bit->getMeshsetId();
453
454 if (id == FIX_BRICK_FACES) {
455
457 fix_faces, true);
458
460 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
461 moab::Interface::UNION);
462
463 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
464 moab::Interface::UNION);
465 fix_faces.merge(adj_ents);
466 } else if (id == FIX_NODES) {
467
469 fix_nodes, true);
470
471 } else if (id == BRICK_PRESSURE_FACES) {
473 meshset, 2, pressure_faces, true);
474
475 } else if (id ==
476 FIX_NODES_Y_DIR) {
478 meshset, 0, fix_second_node, true);
479
480 } else {
482 }
483 }
484
486 3);
488
490
491
496
498
499 DM dm;
501
504
507
509 "PRESSURE");
511
513
514
515 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
517 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
519
520
521 elastic_fe->getRuleHook =
VolRule();
522 pressure_fe->getRuleHook =
FaceRule();
523
524
525
526 elastic_fe->getOpPtrVector().push_back(
new OpK());
527
528 pressure_fe->getOpPtrVector().push_back(
new OpPressure());
529
530 boost::shared_ptr<FEMethod> fix_dofs_fe(
532
533 boost::shared_ptr<FEMethod> null_fe;
534
535
537 dm, simple_interface->
getDomainFEName(), elastic_fe, null_fe, null_fe);
538
540 null_fe);
541
542
544
545
547
548
549 CHKERR DMCreateMatrix(dm, &A);
550
551
552 CHKERR DMCreateGlobalVector(dm, &x);
553
554
555
556 CHKERR VecDuplicate(x, &f);
557
558
559
560 elastic_fe->ksp_B =
A;
561 fix_dofs_fe->ksp_B =
A;
562
563
564
565 fix_dofs_fe->ksp_f =
f;
566 pressure_fe->ksp_f =
f;
567
569 elastic_fe);
570
572
573
574
576
577
578 KSP solver;
579
580
581
582 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
583
584
585
586 CHKERR KSPSetFromOptions(solver);
587
588
589 CHKERR KSPSetOperators(solver, A, A);
590
591
593
594
595 CHKERR KSPSolve(solver, f, x);
596
597
598 CHKERR KSPDestroy(&solver);
599
600
601
602 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
603
604
606
607
608
609 auto get_post_proc_ele = [&]() {
610 auto jac_ptr = boost::make_shared<MatrixDouble>();
611 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
612 auto det_ptr = boost::make_shared<VectorDouble>();
613
614 auto post_proc_ele = boost::make_shared<
616
618 post_proc_ele->getOpPtrVector().push_back(
620 post_proc_ele->getOpPtrVector().push_back(
622 post_proc_ele->getOpPtrVector().push_back(
624
625 auto u_ptr = boost::make_shared<MatrixDouble>();
626 auto grad_ptr = boost::make_shared<MatrixDouble>();
627 post_proc_ele->getOpPtrVector().push_back(
629 post_proc_ele->getOpPtrVector().push_back(
631 grad_ptr));
633
634 post_proc_ele->getOpPtrVector().push_back(
635
637
638 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
639 {},
640
641 {{"U", u_ptr}},
642
643 {{"GRAD", grad_ptr}},
644
645 {})
646
647 );
648
649 return post_proc_ele;
650 };
651
652 if (auto post_proc = get_post_proc_ele()) {
654 post_proc);
655
656 CHKERR post_proc->writeFile(
"out.h5m");
657 }
658
659 {
660 if (flg_test == PETSC_TRUE) {
661 const double x_vec_norm_const = 0.4;
662
663 double norm_check;
664
665
666 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
667 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
668 1.e-10) {
670 "test failed (nrm 2 %6.4e)", norm_check);
671 }
672 }
673 }
674
675
679
680
682
683
684 }
686
687
689
690 return 0;
691}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set compute operator for KSP solver via sub-matrix and IS.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set integration rule to boundary elements.
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Simple interface for fast problem set-up.
MoFEMErrorCode buildProblem()
Build problem.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Volume finite element base.
Set integration rule to volume elements.