390 {
391
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
399 if (!
static_cast<bool>(ifstream(
param_file))) {
400 std::ofstream file(
param_file.c_str(), std::ios::ate);
401 if (file.is_open()) {
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 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"SimpleElasticProblem",
424 "none");
425
426
428 PETSC_NULL);
429
430
431 CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
432 &flg_test, PETSC_NULL);
433
434 ierr = PetscOptionsEnd();
436
438
441
442 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
443
444 enum MyBcTypes {
445 FIX_BRICK_FACES = 1,
446 FIX_NODES = 2,
447 BRICK_PRESSURE_FACES = 3,
448 FIX_NODES_Y_DIR = 4
449 };
450
453 int id =
bit->getMeshsetId();
454
455 if (id == FIX_BRICK_FACES) {
456
458 fix_faces, true);
459
461 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
462 moab::Interface::UNION);
463
464 CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
465 moab::Interface::UNION);
466 fix_faces.merge(adj_ents);
467 } else if (id == FIX_NODES) {
468
470 fix_nodes, true);
471
472 } else if (id == BRICK_PRESSURE_FACES) {
474 meshset, 2, pressure_faces, true);
475
476 } else if (id ==
477 FIX_NODES_Y_DIR) {
479 meshset, 0, fix_second_node, true);
480
481 } else {
483 }
484 }
485
487 3);
489
491
492
497
499
500 DM dm;
502
505
508
510 "PRESSURE");
512
514
515
516 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
518 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
520
521
522 elastic_fe->getRuleHook =
VolRule();
523 pressure_fe->getRuleHook =
FaceRule();
524
525
526
527 elastic_fe->getOpPtrVector().push_back(
new OpK());
528
529 pressure_fe->getOpPtrVector().push_back(
new OpPressure());
530
531 boost::shared_ptr<FEMethod> fix_dofs_fe(
533
534 boost::shared_ptr<FEMethod> null_fe;
535
536
538 dm, simple_interface->
getDomainFEName(), elastic_fe, null_fe, null_fe);
539
541 null_fe);
542
543
545
546
548
549
550 CHKERR DMCreateMatrix(dm, &A);
551
552
553 CHKERR DMCreateGlobalVector(dm, &x);
554
555
556
557 CHKERR VecDuplicate(x, &f);
558
559
560
561 elastic_fe->ksp_B =
A;
562 fix_dofs_fe->ksp_B =
A;
563
564
565
566 fix_dofs_fe->ksp_f =
f;
567 pressure_fe->ksp_f =
f;
568
570 elastic_fe);
571
573
574
575
577
578
579 KSP solver;
580
581
582
583 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
584
585
586
587 CHKERR KSPSetFromOptions(solver);
588
589
590 CHKERR KSPSetOperators(solver, A, A);
591
592
594
595
596 CHKERR KSPSolve(solver, f, x);
597
598
599 CHKERR KSPDestroy(&solver);
600
601
602
603 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
604
605
607
608
609
610 auto get_post_proc_ele = [&]() {
611 auto jac_ptr = boost::make_shared<MatrixDouble>();
612 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
613 auto det_ptr = boost::make_shared<VectorDouble>();
614
615 auto post_proc_ele = boost::make_shared<
617
619 post_proc_ele->getOpPtrVector().push_back(
621 post_proc_ele->getOpPtrVector().push_back(
623 post_proc_ele->getOpPtrVector().push_back(
625
626 auto u_ptr = boost::make_shared<MatrixDouble>();
627 auto grad_ptr = boost::make_shared<MatrixDouble>();
628 post_proc_ele->getOpPtrVector().push_back(
630 post_proc_ele->getOpPtrVector().push_back(
632 grad_ptr));
634
635 post_proc_ele->getOpPtrVector().push_back(
636
638
639 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
640 {},
641
642 {{"U", u_ptr}},
643
644 {{"GRAD", grad_ptr}},
645
646 {})
647
648 );
649
650 return post_proc_ele;
651 };
652
653 if (auto post_proc = get_post_proc_ele()) {
655 post_proc);
656
657 CHKERR post_proc->writeFile(
"out.h5m");
658 }
659
660 {
661 if (flg_test == PETSC_TRUE) {
662 const double x_vec_norm_const = 0.4;
663
664 double norm_check;
665
666
667 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
668 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
669 1.e-10) {
671 "test failed (nrm 2 %6.4e)", norm_check);
672 }
673 }
674 }
675
676
680
681
683
684
685 }
687
688
690
691 return 0;
692}
const std::string default_options
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
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 KSP right hand side evaluation function
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_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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 filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed 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 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 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 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 refernce to pointer of interface.
Volume finite element base.