v0.9.1
Classes | Functions | Variables
simple_elasticity.cpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  OpK
 assemble the stiffness matrix More...
 
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 409 of file simple_elasticity.cpp.

409  {
410 
411  const string default_options = "-ksp_type fgmres \n"
412  "-pc_type lu \n"
413  "-pc_factor_mat_solver_package mumps \n"
414  "-ksp_monitor\n";
415 
416  string param_file = "param_file.petsc";
417  if (!static_cast<bool>(ifstream(param_file))) {
418  std::ofstream file(param_file.c_str(), std::ios::ate);
419  if (file.is_open()) {
420  file << default_options;
421  file.close();
422  }
423  }
424 
425  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
426 
427  // Create mesh database
428  moab::Core mb_instance; // create database
429  moab::Interface &moab = mb_instance; // create interface to database
430 
431  try {
432  // Create MoFEM database and link it to MoAB
433  MoFEM::Core core(moab);
434  MoFEM::Interface &m_field = core;
435 
436  CHKERR DMRegister_MoFEM("DMMOFEM");
437 
438  // Get command line options
439  int order = 3; // default approximation order
440  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
441  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
442  "none");
443 
444  // Set approximation order
445  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
446  PETSC_NULL);
447 
448  // Set testing (used by CTest)
449  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
450  &flg_test, PETSC_NULL);
451 
452  ierr = PetscOptionsEnd();
453  CHKERRG(ierr);
454 
455  Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
456 
457  CHKERR simple_interface->getOptions();
458  CHKERR simple_interface->loadFile();
459 
460  Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
461 
462  enum MyBcTypes {
463  FIX_BRICK_FACES = 1,
464  FIX_NODES = 2,
465  BRICK_PRESSURE_FACES = 3,
466  FIX_NODES_Y_DIR = 4
467  };
468 
470  EntityHandle meshset = bit->getMeshset();
471  int id = bit->getMeshsetId();
472 
473  if (id == FIX_BRICK_FACES) { // brick-faces
474 
475  CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
476  fix_faces, true);
477 
478  Range adj_ents;
479  CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
480  moab::Interface::UNION);
481 
482  CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
483  moab::Interface::UNION);
484  fix_faces.merge(adj_ents);
485  } else if (id == FIX_NODES) { // node(s)
486 
487  CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
488  fix_nodes, true);
489 
490  } else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
491  CHKERR m_field.get_moab().get_entities_by_dimension(
492  meshset, 2, pressure_faces, true);
493 
494  } else if (id ==
495  FIX_NODES_Y_DIR) { // restrained second node in y direction
496  CHKERR m_field.get_moab().get_entities_by_dimension(
497  meshset, 0, fix_second_node, true);
498 
499  } else {
500  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
501  }
502  }
503 
504  CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
505  3);
506  CHKERR simple_interface->setFieldOrder("U", order);
507 
508  CHKERR simple_interface->defineFiniteElements();
509 
510  // Add pressure element
511  CHKERR m_field.add_finite_element("PRESSURE");
512  CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
513  CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
514  CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
515 
516  CHKERR simple_interface->defineProblem();
517 
518  DM dm;
519  CHKERR simple_interface->getDM(&dm);
520 
521  CHKERR DMMoFEMAddElement(dm, "PRESSURE");
522  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
523 
524  CHKERR simple_interface->buildFields();
525  CHKERR simple_interface->buildFiniteElements();
526 
527  CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
528  "PRESSURE");
529  CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
530 
531  CHKERR simple_interface->buildProblem();
532 
533  // Create elements instances
534  boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
535  new VolumeElementForcesAndSourcesCore(m_field));
536  boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
537  new FaceElementForcesAndSourcesCore(m_field));
538 
539  // Set integration rule to elements instances
540  elastic_fe->getRuleHook = VolRule();
541  pressure_fe->getRuleHook = FaceRule();
542 
543  // Add operators to element instances
544  // push operators to elastic_fe
545  elastic_fe->getOpPtrVector().push_back(new OpK());
546  // push operators to pressure_fe
547  pressure_fe->getOpPtrVector().push_back(new OpPressure());
548 
549  boost::shared_ptr<FEMethod> fix_dofs_fe(
550  new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
551 
552  boost::shared_ptr<FEMethod> null_fe;
553 
554  // Set operators for KSP solver
556  dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
557 
558  CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
559  null_fe);
560 
561  // initialise matrix A used as the global stiffness matrix
562  Mat A;
563 
564  // initialise left hand side vector x and right hand side vector f
565  Vec x, f;
566 
567  // allocate memory handled by MoFEM discrete manager for matrix A
568  CHKERR DMCreateMatrix(dm, &A);
569 
570  // allocate memory handled by MoFEM discrete manager for vector x
571  CHKERR DMCreateGlobalVector(dm, &x);
572 
573  // allocate memory handled by MoFEM discrete manager for vector f of the
574  // same size as x
575  CHKERR VecDuplicate(x, &f);
576 
577  // precondition matrix A according to fix_dofs_fe and elastic_fe finite
578  // elements
579  elastic_fe->ksp_B = A;
580  fix_dofs_fe->ksp_B = A;
581 
582  // precondition the right hand side vector f according to fix_dofs_fe and
583  // elastic_fe finite elements
584  fix_dofs_fe->ksp_f = f;
585  pressure_fe->ksp_f = f;
586 
587  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
588  elastic_fe);
589 
590  CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
591 
592  // This is done because only post processor is implemented in the
593  // ApplyDirichletBc struct
594  CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
595 
596  // make available a KSP solver
597  KSP solver;
598 
599  // make the solver available for parallel computing by determining its MPI
600  // communicator
601  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
602 
603  // making available running all options available for KSP solver in running
604  // command
605  CHKERR KSPSetFromOptions(solver);
606 
607  // set A matrix with preconditioner
608  CHKERR KSPSetOperators(solver, A, A);
609 
610  // set up the solver data strucure for the iterative solver
611  CHKERR KSPSetUp(solver);
612 
613  // solve the system of linear equations
614  CHKERR KSPSolve(solver, f, x);
615 
616  // destroy solver no needed any more
617  CHKERR KSPDestroy(&solver);
618 
619  // make vector x available for parallel computations for visualization
620  // context
621  VecView(x, PETSC_VIEWER_STDOUT_WORLD);
622 
623  // save solution in vector x on mesh
624  CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
625 
626  // Set up post-processor. It is some generic implementation of finite
627  // element
628  PostProcVolumeOnRefinedMesh post_proc(m_field);
629  // Add operators to the elements, starting with some generic
630  CHKERR post_proc.generateReferenceElementMesh();
631 
632  CHKERR post_proc.addFieldValuesPostProc("U");
633 
634  CHKERR post_proc.addFieldValuesGradientPostProc("U");
635 
637  dm, simple_interface->getDomainFEName().c_str(), &post_proc);
638 
639  // write output
640  CHKERR post_proc.writeFile("out.h5m");
641 
642  {
643  if (flg_test == PETSC_TRUE) {
644  const double x_vec_norm_const = 0.4;
645  // Check norm_1 value
646  double norm_check;
647  // Takes maximal element of the vector, that should be maximal
648  // displacement at the end of the bar
649  CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
650  if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const > 1.e-10) {
651  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
652  "test failed (nrm 2 %6.4e)", norm_check);
653  }
654  }
655  }
656 
657  // free memory handled by mofem discrete manager for A, x and f
658  CHKERR MatDestroy(&A);
659  CHKERR VecDestroy(&x);
660  CHKERR VecDestroy(&f);
661 
662  // free memory allocated for mofem discrete manager
663  CHKERR DMDestroy(&dm);
664 
665  // This is a good reference for the future
666  }
667  CATCH_ERRORS;
668 
669  // finish work cleaning memory, getting statistics, etc
671 
672  return 0;
673 }
Set integration rule.
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
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:633
Deprecated interface functions.
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
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:707
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:387
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: DMMMoFEM.cpp:556
virtual moab::Interface & get_moab()=0
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
Core (interface) class.
Definition: Core.hpp:50
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: DMMMoFEM.cpp:597
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
static char help[]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:659
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:491
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:457
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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:269
Post processing.
Set integration rule to boundary elements.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:478
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:602
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:281
constexpr int order
assemble the stiffness matrix
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
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 CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:483
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:460
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:77

Variable Documentation

◆ help

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

Definition at line 26 of file simple_elasticity.cpp.