v0.14.0
Loading...
Searching...
No Matches
Classes | Functions | Variables
simple_elasticity.cpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  OpK
 
struct  OpPressure
 
struct  ApplyDirichletBc
 
struct  VolRule
 Set integration rule. More...
 
struct  FaceRule
 Set integration rule to boundary elements. More...
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help []
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 390 of file simple_elasticity.cpp.

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
407 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
408
409 // Create mesh database
410 moab::Core mb_instance; // create database
411 moab::Interface &moab = mb_instance; // create interface to database
412
413 try {
414 // Create MoFEM database and link it to MoAB
415 MoFEM::Core core(moab);
416 MoFEM::Interface &m_field = core;
417
418 CHKERR DMRegister_MoFEM("DMMOFEM");
419
420 // Get command line options
421 int order = 3; // default approximation order
422 PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
423 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
424 "none");
425
426 // Set approximation order
427 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
428 PETSC_NULL);
429
430 // Set testing (used by CTest)
431 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
432 &flg_test, PETSC_NULL);
433
434 ierr = PetscOptionsEnd();
435 CHKERRG(ierr);
436
437 Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
438
439 CHKERR simple_interface->getOptions();
440 CHKERR simple_interface->loadFile();
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
452 EntityHandle meshset = bit->getMeshset();
453 int id = bit->getMeshsetId();
454
455 if (id == FIX_BRICK_FACES) { // brick-faces
456
457 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
458 fix_faces, true);
459
460 Range adj_ents;
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) { // node(s)
468
469 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
470 fix_nodes, true);
471
472 } else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
473 CHKERR m_field.get_moab().get_entities_by_dimension(
474 meshset, 2, pressure_faces, true);
475
476 } else if (id ==
477 FIX_NODES_Y_DIR) { // restrained second node in y direction
478 CHKERR m_field.get_moab().get_entities_by_dimension(
479 meshset, 0, fix_second_node, true);
480
481 } else {
482 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
483 }
484 }
485
486 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
487 3);
488 CHKERR simple_interface->setFieldOrder("U", order);
489
490 CHKERR simple_interface->defineFiniteElements();
491
492 // Add pressure element
493 CHKERR m_field.add_finite_element("PRESSURE");
494 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
495 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
496 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
497
498 CHKERR simple_interface->defineProblem();
499
500 DM dm;
501 CHKERR simple_interface->getDM(&dm);
502
503 CHKERR DMMoFEMAddElement(dm, "PRESSURE");
504 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
505
506 CHKERR simple_interface->buildFields();
507 CHKERR simple_interface->buildFiniteElements();
508
509 CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
510 "PRESSURE");
511 CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
512
513 CHKERR simple_interface->buildProblem();
514
515 // Create elements instances
516 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
518 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
520
521 // Set integration rule to elements instances
522 elastic_fe->getRuleHook = VolRule();
523 pressure_fe->getRuleHook = FaceRule();
524
525 // Add operators to element instances
526 // push operators to elastic_fe
527 elastic_fe->getOpPtrVector().push_back(new OpK());
528 // push operators to pressure_fe
529 pressure_fe->getOpPtrVector().push_back(new OpPressure());
530
531 boost::shared_ptr<FEMethod> fix_dofs_fe(
532 new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
533
534 boost::shared_ptr<FEMethod> null_fe;
535
536 // Set operators for KSP solver
538 dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
539
540 CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
541 null_fe);
542
543 // initialise matrix A used as the global stiffness matrix
544 Mat A;
545
546 // initialise left hand side vector x and right hand side vector f
547 Vec x, f;
548
549 // allocate memory handled by MoFEM discrete manager for matrix A
550 CHKERR DMCreateMatrix(dm, &A);
551
552 // allocate memory handled by MoFEM discrete manager for vector x
553 CHKERR DMCreateGlobalVector(dm, &x);
554
555 // allocate memory handled by MoFEM discrete manager for vector f of the
556 // same size as x
557 CHKERR VecDuplicate(x, &f);
558
559 // precondition matrix A according to fix_dofs_fe and elastic_fe finite
560 // elements
561 elastic_fe->ksp_B = A;
562 fix_dofs_fe->ksp_B = A;
563
564 // precondition the right hand side vector f according to fix_dofs_fe and
565 // elastic_fe finite elements
566 fix_dofs_fe->ksp_f = f;
567 pressure_fe->ksp_f = f;
568
569 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
570 elastic_fe);
571
572 CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
573
574 // This is done because only post processor is implemented in the
575 // ApplyDirichletBc struct
576 CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
577
578 // make available a KSP solver
579 KSP solver;
580
581 // make the solver available for parallel computing by determining its MPI
582 // communicator
583 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
584
585 // making available running all options available for KSP solver in running
586 // command
587 CHKERR KSPSetFromOptions(solver);
588
589 // set A matrix with preconditioner
590 CHKERR KSPSetOperators(solver, A, A);
591
592 // set up the solver data strucure for the iterative solver
593 CHKERR KSPSetUp(solver);
594
595 // solve the system of linear equations
596 CHKERR KSPSolve(solver, f, x);
597
598 // destroy solver no needed any more
599 CHKERR KSPDestroy(&solver);
600
601 // make vector x available for parallel computations for visualization
602 // context
603 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
604
605 // save solution in vector x on mesh
606 CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
607
608 // Set up post-processor. It is some generic implementation of finite
609 // element
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 // Add operators to the elements, starting with some generic
618 constexpr int SPACE_DIM = 3;
619 post_proc_ele->getOpPtrVector().push_back(
620 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
621 post_proc_ele->getOpPtrVector().push_back(
622 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
623 post_proc_ele->getOpPtrVector().push_back(
624 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
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
637 new OpPPMap(
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()) {
654 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
655 post_proc);
656 // write output
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 // Check norm_1 value
664 double norm_check;
665 // Takes maximal element of the vector, that should be maximal
666 // displacement at the end of the bar
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) {
670 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
671 "test failed (nrm 2 %6.4e)", norm_check);
672 }
673 }
674 }
675
676 // free memory handled by mofem discrete manager for A, x and f
677 CHKERR MatDestroy(&A);
678 CHKERR VecDestroy(&x);
679 CHKERR VecDestroy(&f);
680
681 // free memory allocated for mofem discrete manager
682 CHKERR DMDestroy(&dm);
683
684 // This is a good reference for the future
685 }
687
688 // finish work cleaning memory, getting statistics, etc
690
691 return 0;
692}
const std::string default_options
std::string param_file
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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
Definition: DMMoFEM.cpp:623
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
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.
Definition: DMMoFEM.cpp:664
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.
auto bit
set bit
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static char help[]
Set integration rule to boundary elements.
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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.
Definition: Simple.hpp:27
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:597
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:380
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:571
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.
Definition: Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:482
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:452
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Set integration rule.

Variable Documentation

◆ help

char help[]
static
Initial value:
= "-order approximation order\n"
"\n"

Definition at line 13 of file simple_elasticity.cpp.