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 
  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    PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
  424                             "none");
  425 
  426    
  428                           PETSC_NULLPTR);
  429 
  430    
  431    CHKERR PetscOptionsBool(
"-test", 
"if true is ctest", 
"", flg_test,
 
  432                            &flg_test, PETSC_NULLPTR);
  433 
  434    PetscOptionsEnd();
  435 
  437 
  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 
  452      int id = 
bit->getMeshsetId();
 
  453 
  454      if (id == FIX_BRICK_FACES) { 
  455 
  457                                                            fix_faces, true);
  458 
  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) { 
  467 
  469                                                            fix_nodes, true);
  470 
  471      } else if (id == BRICK_PRESSURE_FACES) { 
  473            meshset, 2, pressure_faces, true);
  474 
  475      } else if (id ==
  476                 FIX_NODES_Y_DIR) { 
  478            meshset, 0, fix_second_node, true);
  479 
  480      } else {
  482      }
  483    }
  484 
  486                                            3);
  488 
  490 
  491    
  496 
  498 
  499    DM dm;
  501 
  504 
  507 
  509                                                     "PRESSURE");
  511 
  513 
  514    
  515    boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
  517    boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
  519 
  520    
  521    elastic_fe->getRuleHook = 
VolRule();
 
  522    pressure_fe->getRuleHook = 
FaceRule();
 
  523 
  524    
  525    
  526    elastic_fe->getOpPtrVector().push_back(
new OpK());
 
  527    
  528    pressure_fe->getOpPtrVector().push_back(
new OpPressure());
 
  529 
  530    boost::shared_ptr<FEMethod> fix_dofs_fe(
  532 
  533    boost::shared_ptr<FEMethod> null_fe;
  534 
  535    
  537        dm, simple_interface->
getDomainFEName(), elastic_fe, null_fe, null_fe);
 
  538 
  540                                   null_fe);
  541 
  542    
  544 
  545    
  547 
  548    
  549    CHKERR DMCreateMatrix(dm, &A);
 
  550 
  551    
  552    CHKERR DMCreateGlobalVector(dm, &x);
 
  553 
  554    
  555    
  556    CHKERR VecDuplicate(x, &f);
 
  557 
  558    
  559    
  560    elastic_fe->ksp_B = 
A;
 
  561    fix_dofs_fe->ksp_B = 
A;
 
  562 
  563    
  564    
  565    fix_dofs_fe->ksp_f = 
f;
 
  566    pressure_fe->ksp_f = 
f;
 
  567 
  569                                    elastic_fe);
  570 
  572 
  573    
  574    
  576 
  577    
  578    KSP solver;
  579 
  580    
  581    
  582    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
 
  583 
  584    
  585    
  586    CHKERR KSPSetFromOptions(solver);
 
  587 
  588    
  589    CHKERR KSPSetOperators(solver, A, A);
 
  590 
  591    
  593 
  594    
  595    CHKERR KSPSolve(solver, f, x);
 
  596 
  597    
  598    CHKERR KSPDestroy(&solver);
 
  599 
  600    
  601    
  602    VecView(x, PETSC_VIEWER_STDOUT_WORLD);
  603 
  604    
  606 
  607    
  608    
  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      
  618      post_proc_ele->getOpPtrVector().push_back(
  620      post_proc_ele->getOpPtrVector().push_back(
  622      post_proc_ele->getOpPtrVector().push_back(
  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 
  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()) {
  654                                      post_proc);
  655      
  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        
  663        double norm_check;
  664        
  665        
  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) {
  670                   "test failed (nrm 2 %6.4e)", norm_check);
  671        }
  672      }
  673    }
  674 
  675    
  679 
  680    
  682 
  683    
  684  }
  686 
  687  
  689 
  690  return 0;
  691}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ 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_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.
const FTensor::Tensor2< T, Dim, Dim > Vec
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 field rank 0, i.e. vector field.
Specialization for double precision vector field values calculation.
Operator for inverting matrices at integration points.
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 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 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 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 reference to pointer of interface.
Volume finite element base.
Set integration rule to volume elements.