#include <MoFEM.hpp>
#include <Hdiv.hpp>

Classes
struct	SomeUserPolynomialBase
	Class used to calculate base functions at integration points. More...

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

Variables
static PetscBool	quiet = PETSC_FALSE

static PetscBool	base_cache = PETSC_FALSE

static char	help [] = "...\n\n"

Function Documentation

◆ main()

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

Simple user data operator which main purpose is to print values of base functions at integration points.

Examples: forces_and_sources_testing_users_base.cpp.

Definition at line 172 of file forces_and_sources_testing_users_base.cpp.

                                  {
  
   // Initialise MoFEM, MPI and petsc
   MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
  
   try {
  
     // create moab
     moab::Core mb_instance;
     // get interface to moab databse
     moab::Interface &moab = mb_instance;
  
     // get file
     PetscBool flg = PETSC_TRUE;
     char mesh_file_name[255];
 #if PETSC_VERSION_GE(3, 6, 4)
     CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
                                  255, &flg);
 #else
     CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
                                  mesh_file_name, 255, &flg);
 #endif
     if (flg != PETSC_TRUE) {
       SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "*** ERROR -my_file (MESH FILE NEEDED)");
     }
  
     // create MoFEM database
     MoFEM::Core core(moab);
     // get interface to moab database
     MoFEM::Interface &m_field = core;
  
     // load mesh file
     const char *option;
     option = "";
     CHKERR moab.load_file(mesh_file_name, 0, option);
  
     // set bit refinement level
     CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
         0, 3, BitRefLevel().set(0));
  
     // Create fields, field "FIELD_CGG" has user base, it means that recipe how
     // to construct approximation is provided by user. Is set that user provided
     // base is in h-div space.
     CHKERR m_field.add_field("FILED_CGG", HDIV, USER_BASE, 1);
     CHKERR m_field.add_field("FILED_RT", HDIV, DEMKOWICZ_JACOBI_BASE, 1);
  
     // get access to "FIELD_CGG" data structure
     auto field_ptr = m_field.get_field_structure("FILED_CGG");
     // get table associating number of dofs to entities depending on
     // approximation order set on those entities.
     auto field_order_table =
         const_cast<Field *>(field_ptr)->getFieldOrderTable();
  
     // function set zero number of dofs
     auto get_cgg_bubble_order_zero = [](int p) { return 0; };
     // function set non-zero number of dofs on tetrahedrons
     auto get_cgg_bubble_order_face = [](int p) {
       return NBFACETRI_DEMKOWICZ_HDIV(p);
     };
     auto get_cgg_bubble_order_tet = [](int p) {
       return NBVOLUMETET_DEMKOWICZ_HDIV(p);
     };
     field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
     field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
     field_order_table[MBTRI] = get_cgg_bubble_order_face;
     field_order_table[MBTET] = get_cgg_bubble_order_tet;
     const_cast<Field *>(field_ptr)->rebuildDofsOrderMap();
  
     auto &dof_order_map = field_ptr->getDofOrderMap(MBTET);
     for(auto d = 0; d!=10; ++d) {
       MOFEM_LOG("WORLD", Sev::noisy) << "dof " << dof_order_map[d];
     }
  
     CHKERR m_field.add_finite_element("FE");
  
     // define rows/cols and element data
     CHKERR m_field.modify_finite_element_add_field_row("FE", "FILED_CGG");
     CHKERR m_field.modify_finite_element_add_field_col("FE", "FILED_CGG");
     CHKERR m_field.modify_finite_element_add_field_data("FE", "FILED_CGG");
     CHKERR m_field.modify_finite_element_add_field_row("FE", "FILED_RT");
     CHKERR m_field.modify_finite_element_add_field_col("FE", "FILED_RT");
     CHKERR m_field.modify_finite_element_add_field_data("FE", "FILED_RT");
  
     // add problem
     CHKERR m_field.add_problem("PROBLEM");
  
     // set finite elements for problem
     CHKERR m_field.modify_problem_add_finite_element("PROBLEM", "FE");
     // set refinement level for problem
     CHKERR m_field.modify_problem_ref_level_add_bit("PROBLEM",
                                                     BitRefLevel().set(0));
  
     // meshset consisting all entities in mesh
     EntityHandle root_set = moab.get_root_set();
     // add entities to field
     CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FILED_CGG");
     CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FILED_RT");
     // add entities to finite element
     CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "FE");
  
     // set app. order
     int order = 3;
     CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
     CHKERR m_field.set_field_order(root_set, MBTRI, "FILED_CGG", order);
     CHKERR m_field.set_field_order(root_set, MBTET, "FILED_CGG", order);
     CHKERR m_field.set_field_order(root_set, MBTRI, "FILED_RT", order);
     CHKERR m_field.set_field_order(root_set, MBTET, "FILED_RT", order);
  
     /****/
     // build database
     // build field
     CHKERR m_field.build_fields();
     // build finite elemnts
     CHKERR m_field.build_finite_elements();
     // build adjacencies
     CHKERR m_field.build_adjacencies(BitRefLevel().set(0));
  
     // build problem
     CHKERR m_field.getInterface<ProblemsManager>()->buildProblem("PROBLEM",
                                                                  true);
     // dofs partitioning
     CHKERR m_field.getInterface<ProblemsManager>()->partitionSimpleProblem(
         "PROBLEM");
     CHKERR m_field.getInterface<ProblemsManager>()->partitionFiniteElements(
         "PROBLEM");
     // what are ghost nodes, see Petsc Manual
     CHKERR m_field.getInterface<ProblemsManager>()->partitionGhostDofs(
         "PROBLEM");
  
     typedef tee_device<std::ostream, std::ofstream> TeeDevice;
     typedef stream<TeeDevice> TeeStream;
  
     std::ofstream ofs("forces_and_sources_testing_users_base.txt");
     TeeDevice my_tee(std::cout, ofs);
     TeeStream my_split(my_tee);
  
     /**
      * Simple user data operator which main purpose is to print values
      * of base functions at integration points.
      *
      */
     struct MyOp1 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
  
       TeeStream &my_split;
       MyOp1(const std::string &row_field, const std::string &col_field,
             TeeStream &_my_split, char type)
           : VolumeElementForcesAndSourcesCore::UserDataOperator(
                 row_field, col_field, type),
             my_split(_my_split) {
         sYmm = false;
       }
  
       MoFEMErrorCode doWork(int side, EntityType type,
                             EntitiesFieldData::EntData &data) {
         MoFEMFunctionBeginHot;
         if (data.getIndices().empty()) {
           MoFEMFunctionReturnHot(0);
         }
         if (!quiet) {
           my_split << rowFieldName << endl;
           my_split << "side: " << side << " type: " << type << std::endl;
           my_split << data << endl;
           my_split << data.getN() << endl;
           my_split << endl;
         }
         MoFEMFunctionReturnHot(0);
       }
  
       MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                             EntityType col_type,
                             EntitiesFieldData::EntData &row_data,
                             EntitiesFieldData::EntData &col_data) {
         MoFEMFunctionBeginHot;
         if (row_data.getIndices().empty())
           MoFEMFunctionReturnHot(0);
         if (col_data.getIndices().empty())
           MoFEMFunctionReturnHot(0);
         if (!quiet) {
           my_split << rowFieldName << " : " << colFieldName << endl;
           my_split << "row side: " << row_side << " row_type: " << row_type
                    << std::endl;
           my_split << "col side: " << col_side << " col_type: " << col_type
                    << std::endl;
           my_split << row_data.getIndices().size() << " : "
                    << col_data.getIndices().size() << endl;
           my_split << endl;
         }
         MoFEMFunctionReturnHot(0);
       }
     };
  
     // create finite element instance
     VolumeElementForcesAndSourcesCore fe1(m_field);
     // set class needed to construct user approximation base
     fe1.getUserPolynomialBase() =
         boost::shared_ptr<BaseFunction>(new SomeUserPolynomialBase());
  
     CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quiet", &quiet, PETSC_NULL);
  
     // push user data operators
     fe1.getOpPtrVector().push_back(
         new MyOp1("FILED_CGG", "FILED_CGG", my_split,
                   ForcesAndSourcesCore::UserDataOperator::OPROW));
     fe1.getOpPtrVector().push_back(
         new MyOp1("FILED_CGG", "FILED_RT", my_split,
                   ForcesAndSourcesCore::UserDataOperator::OPROWCOL));
  
     CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-base_cache", &base_cache,
                                PETSC_NULL);
  
     if (base_cache) {
       if (!TetPolynomialBase::switchCacheBaseFace<HDIV>(DEMKOWICZ_JACOBI_BASE,
                                                        &fe1)) {
         TetPolynomialBase::switchCacheBaseFace<HDIV>(DEMKOWICZ_JACOBI_BASE,
                                                     &fe1);
       }
       if (!TetPolynomialBase::switchCacheBaseInterior<HDIV>(
               DEMKOWICZ_JACOBI_BASE, &fe1)) {
         TetPolynomialBase::switchCacheBaseInterior<HDIV>(DEMKOWICZ_JACOBI_BASE,
                                                       &fe1);
       }
     }
  
     // iterate over finite elements, and execute user data operators on each
     // of them
     CHKERR m_field.loop_finite_elements("PROBLEM", "FE", fe1);
  
     if (base_cache) {
       if (TetPolynomialBase::switchCacheBaseFace<HDIV>(DEMKOWICZ_JACOBI_BASE,
                                                     &fe1)) {
       }
       if (TetPolynomialBase::switchCacheBaseInterior<HDIV>(DEMKOWICZ_JACOBI_BASE,
                                                           &fe1)) {
       };
     }
   }
   CATCH_ERRORS;
  
   CHKERR MoFEM::Core::Finalize();
  
   return 0;
 }