v0.15.0
Loading...
Searching...
No Matches
simple_elasticity.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

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

Variable Documentation

◆ help

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

Definition at line 13 of file simple_elasticity.cpp.