v0.14.0
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[] 
)
Examples
simple_elasticity.cpp.

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(
517  new VolumeElementForcesAndSourcesCore(m_field));
518  boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
519  new FaceElementForcesAndSourcesCore(m_field));
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  }
686  CATCH_ERRORS;
687 
688  // finish work cleaning memory, getting statistics, etc
690 
691  return 0;
692 }

Variable Documentation

◆ help

char help[]
static
Initial value:
= "-order approximation order\n"
"\n"
Examples
simple_elasticity.cpp.

Definition at line 13 of file simple_elasticity.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
OpK
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradSymTensorGrad< BASE_DIM, SPACE_DIM, SPACE_DIM, 0 > OpK
[Define entities]
Definition: elastic.cpp:39
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMKSPSetComputeRHS
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:637
EntityHandle
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::Simple::buildProblem
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:615
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
VolRule
Set integration rule.
Definition: simple_interface.cpp:88
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DMMoFEMKSPSetComputeOperators
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:678
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateHOJac
Definition: HODataOperators.hpp:267
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:693
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::add_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
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::Simple::addDomainField
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
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
help
static char help[]
Definition: simple_elasticity.cpp:13
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
MoFEM::PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:670
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:491
MoFEM::Simple::buildFiniteElements
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:589
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
Range
MoFEM::CoreTmp< 0 >::Initialize
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Simple::defineFiniteElements
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:398
MoFEM::Simple::buildFields
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:500
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
ApplyDirichletBc
Definition: simple_elasticity.cpp:301
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
OpPressure
Definition: simple_elasticity.cpp:229
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::Simple::defineProblem
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:470
MoFEM::DMoFEMLoopFiniteElements
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:586
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698