#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.

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;
}