Simple.hpp

Go to the documentation of this file.

/** \file Simple.hpp
 * \brief Header file for simple interface
 * \ingroup mofem_simple_interface
 *
 * Make simplified interface, to speed up problem setup and analysis.
 * See discussion here
 * <a
 * href=https://groups.google.com/d/msg/mofem-group/Vkc00aia4dU/o9RF3ZmPAAAJ>link
 * to google groups</a>
 *
 */
 
#ifndef __SIMPLE_HPP__
#define __SIMPLE_HPP__
 
#include "UnknownInterface.hpp"
 
namespace MoFEM {
 
template <int DIM> struct ParentFiniteElementAdjacencyFunction;
template <int DIM> struct ParentFiniteElementAdjacencyFunctionSkeleton;
 
/**
 * \brief Simple interface for fast problem set-up
 * \ingroup mofem_simple_interface
 */
struct Simple : public UnknownInterface {
 
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  Simple(const MoFEM::Core &core);
  virtual ~Simple() = default;
 
  /**
   * \brief get options
   * 
   * Retrieves command-line options and configuration settings for the Simple interface
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode getOptions();
 
  typedef boost::function<MoFEMErrorCode(Interface &, const char *,
                                         const char *)>
      LoadFileFunc;
  /**
   * \brief Set function for loading mesh file
   * 
   * Default function for loading mesh files using MOAB interface
   * 
   * @param  m_field MoFEM interface object
   * @param  file_name Name of the mesh file to load
   * @param  options Loader options string for MOAB
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  static MoFEMErrorCode defaultLoadFileFunc(Interface &m_field,
                                            const char *file_name,
                                            const char *options) {
    // Default behavior using the member moab.load_file
    return m_field.get_moab().load_file(file_name, 0, options);
  }
 
  /**
   * \brief Load mesh file
   * 
   * Loads mesh file with specified options and custom loading function
   * 
   * @param  options File load options string passed to the loader
   * @param  mesh_file_name Name of mesh file to load (if empty, uses default or command line)
   * @param  loadFunc Function for loading mesh file (default: defaultLoadFileFunc)
   *
   * \note If bitRefLevel is set to any, bit ref level of loaded entities is not
   * changed. After mesh is load, bit ref level should be set to create problem.
   * Default setting of bit ref level is on first bit, and if is set all mesh
   * entities on load are set to set level.
   *
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode loadFile(const std::string options,
                          const std::string mesh_file_name,
                          LoadFileFunc loadFunc = defaultLoadFileFunc);
  /**
   * \brief Load mesh file with parallel options if number of cores > 1
   * 
   * Automatically selects appropriate loading options for parallel execution
   * 
   * @param  mesh_file_name Name of mesh file to load (if empty, uses default or command line value)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode loadFile(const std::string mesh_file_name = "");
 
  /**
   * \brief Add field on domain
   * 
   * Adds a field to the domain entities (highest dimension entities in the mesh)
   * 
   * @param  name              Name of the field
   * @param  space             Function space (L2, H1, Hdiv, Hcurl)
   * @param  base              Approximation base (see FieldApproximationBase)
   * @param  nb_of_coefficients Number of field coefficients per DOF
   * @param  tag_type          Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param  bh                Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param  verb              Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  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);
 
  /**
   * \brief Add broken field on domain
   * 
   * Adds a discontinuous field to the domain entities with element-wise DOF mapping
   * 
   * @param  name              Name of the field
   * @param  space             Function space (L2, H1, Hdiv, Hcurl)
   * @param  base              Approximation base (see FieldApproximationBase)
   * @param  nb_of_coefficients Number of field coefficients per DOF
   * @param  tag_type          Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param  bh                Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param  verb              Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode
  addDomainBrokenField(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);
 
  /**
   * \brief Add field on boundary
   * 
   * Adds a field to the boundary entities (entities on the mesh boundary)
   * 
   * @param  name              Name of the field
   * @param  space             Function space (L2, H1, Hdiv, Hcurl)
   * @param  base              Approximation base (see FieldApproximationBase)
   * @param  nb_of_coefficients Number of field coefficients per DOF
   * @param  tag_type          Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param  bh                Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param  verb              Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode
  addBoundaryField(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);
 
  /**
   * \brief Add field on skeleton
   * 
   * Adds a field to the skeleton entities (internal interfaces between domain elements)
   * 
   * @param  name              Name of the field
   * @param  space             Function space (L2, H1, Hdiv, Hcurl)
   * @param  base              Approximation base (see FieldApproximationBase)
   * @param  nb_of_coefficients Number of field coefficients per DOF
   * @param  tag_type          Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param  bh                Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param  verb              Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode
  addSkeletonField(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);
 
  /**
   * \brief Add data field
   * 
   * Adds a field for storing data values (typically for material properties or constants)
   * 
   * @param  name              Name of the field
   * @param  space             Function space (L2, H1, Hdiv, Hcurl)
   * @param  base              Approximation base (see FieldApproximationBase)
   * @param  nb_of_coefficients Number of field coefficients per DOF
   * @param  tag_type          Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param  bh                Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param  verb              Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode addDataField(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);
 
  /**
   * @brief Add meshset field
   * 
   * Add fields managed by meshset finite element for global DOFs
   *
   * @param name Name of the field
   * @param space Function space (L2, H1, Hdiv, Hcurl)
   * @param base Approximation base (see FieldApproximationBase)
   * @param nb_of_coefficients Number of field coefficients per DOF
   * @param tag_type Tag storage type (MB_TAG_DENSE is faster, MB_TAG_SPARSE more flexible)
   * @param bh Behavior if field exists (MF_EXCL throws error, MF_ZERO ignores)
   * @param verb Verbosity level (-1 for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode
  addMeshsetField(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);
 
  /**
   * @brief Remove field from domain
   *
   * @param name Name of the field to remove from domain
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode removeDomainField(const std::string name);
 
  /**
   * @brief Remove field from boundary
   *
   * @param name Name of the field to remove from boundary
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode removeBoundaryField(const std::string name);
 
  /**
   * @brief Remove field from skeleton
   *
   * @param name Name of the field to remove from skeleton
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode removeSkeletonField(const std::string name);
 
  /**
   * \brief Define finite elements
   * 
   * Creates and defines domain, boundary, and skeleton finite elements based on added fields
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode defineFiniteElements();
 
  /**
   * \brief define problem
   * 
   * Creates the problem and associates finite elements with it
   * 
   * @param is_partitioned Whether the mesh should be partitioned for parallel execution
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode defineProblem(const PetscBool is_partitioned = PETSC_TRUE);
 
  /**
   * \brief Set field order
   * 
   * Sets the approximation order for a field on specified entities
   * 
   * @param  field_name Name of the field
   * @param  order Approximation order to set
   * @param  ents Range of entities to apply order to (NULL for all entities)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode setFieldOrder(const std::string field_name, const int order,
                               const Range *ents = NULL);
 
  /**
   * \brief Build fields
   * 
   * Builds the field data structures and creates DOFs
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode buildFields();
 
  /**
   * \brief Build finite elements
   * 
   * Builds the finite element data structures and adjacencies
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode buildFiniteElements();
 
  /**
   * @brief Set the skeleton adjacency object
   *
   * Sets up adjacency relationships for skeleton finite elements
   * 
   * @param dim Dimension for skeleton adjacency (-1 for automatic detection)
   * @param fe_name Name of finite element to set adjacency for (empty for default)
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode setSkeletonAdjacency(int dim = -1, std::string fe_name = "");
 
  /**
   * \brief Build problem
   * 
   * Builds the problem data structures and sets up DOF numbering
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode buildProblem();
 
  /**
   * \brief Setup problem
   * 
   * Complete setup of the problem including fields, finite elements, and problem definition
   * 
   * @param is_partitioned Whether the mesh should be partitioned for parallel execution
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode setUp(const PetscBool is_partitioned = PETSC_TRUE);
 
  /**
   * @brief Rebuild internal MoFEM data structures
   *
   * Call this function after you add field or remove it.
   *
   * \note If you add field, or remove it, finite element and problem needs to
   * be rebuild. However DM can remain the same.
   *
   * @param only_dm If true, only rebuild DM without rebuilding fields and finite elements
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode reSetUp(bool only_dm = false);
 
  /**
   * \brief Get DM
   * 
   * Retrieves the PETSc Discrete Manager (DM) object
   * 
   * @param  dm Pointer to store the DM object
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode getDM(DM *dm);
 
  /**
   * @brief Return smart DM object
   *
   * Returns a smart pointer to the DM object that handles automatic cleanup
   * 
   * \code
   * {
   *    auto dm  = simple_interface->getDM();
   *
   *    // ...
   *
   *    // dm is automatically destroyed when out of the scope
   * }
   * \endcode
   *
   * @return SmartPetscObj<DM> Smart pointer to DM object
   */
  inline SmartPetscObj<DM> getDM() { return dM; }
 
  /**
   * @brief Get the problem dimension
   *
   * Problem dimension is determined by highest dimension of entities on the
   * mesh.
   *
   * @return int Problem dimension (1D, 2D, or 3D)
   */
  inline int getDim() const { return dIm; }
 
  /**
   * @brief Set the problem dimension
   *
   * @param dim Problem dimension to set (1, 2, or 3)
   */
  void setDim(int dim) { dIm = dim; };
 
  /**
   * @deprecated Use getMeshset
   *
   * @return EntityHandle& Reference to domain meshset handle
   */
  DEPRECATED inline EntityHandle &getMeshSet() { return meshSet; }
 
  /**
   * @brief Get the MeshSet object
   *
   * @return EntityHandle& Reference to domain meshset handle
   */
  inline EntityHandle &getMeshset() { return meshSet; }
 
  /**
   * @brief Get the BoundaryMeshSet object
   *
   * @return EntityHandle& Reference to boundary meshset handle
   */
  inline EntityHandle &getBoundaryMeshSet() { return boundaryMeshset; }
 
  /**
   * @brief Get the SkeletonMeshSet object
   *
   * @return EntityHandle& Reference to skeleton meshset handle
   */
  inline EntityHandle &getSkeletonMeshSet() { return skeletonMeshset; }
 
  /**
   * @brief Get the BitRefLevel
   *
   * @return BitRefLevel& Reference to bit refinement level
   */
  inline BitRefLevel &getBitRefLevel() { return bitLevel; }
 
  /**
   * @brief Get the BitRefLevelMask
   *
   * @return BitRefLevel& Reference to bit refinement level mask
   */
  inline BitRefLevel &getBitRefLevelMask() { return bitLevelMask; }
 
  /**
   * @brief Get the Domain FE Name
   *
   * @return const std::string
   */
  inline const std::string getDomainFEName() const { return domainFE; }
 
  /**
   * @brief Get the Boundary FE Name
   *
   * @return const std::string
   */
  inline const std::string getBoundaryFEName() const { return boundaryFE; }
 
  /**
   * @brief Get the Skeleton FE Name
   *
   * @return const std::string
   */
  inline const std::string getSkeletonFEName() const { return skeletonFE; }
 
  /**
   * @brief Get the Problem Name
   *
   * @return const std::string
   */
  inline const std::string getProblemName() const { return nameOfProblem; }
 
  /**
   * @brief Get the Domain FE Name
   *
   * @return std::string&
   */
  inline std::string &getDomainFEName() { return domainFE; }
 
  /**
   * @brief Get the Boundary FE Name
   *
   * @return std::string&
   */
  inline std::string &getBoundaryFEName() { return boundaryFE; }
 
  /**
   * @brief Get the Skeleton FE Name
   *
   * @return std::string&
   */
  inline std::string &getSkeletonFEName() { return skeletonFE; }
 
  /**
   * @brief Get the Skeleton FE Name
   *
   * @return std::string&
   */
  inline std::string &getMeshsetFEName() { return meshsetFE; }
 
  /**
   * @brief Range of entities added to meshset finite element
   *
   * Meshset element have not defined integration rule, if you like to integrate
   * you have to add mesh entity to meshset finite element. In principle
   * meshset element become standard FE.
   *
   * @return std::vector<std::string>&
   */
  inline auto &getMeshsetFiniteElementEntities() {
    return meshsetFiniteElementEntities;
  }
 
  /**
   * @brief Get the Problem Name
   *
   * @return std::string&
   */
  inline std::string &getProblemName() { return nameOfProblem; }
 
  /**
   * @brief Get the Other Finite Elements
   *
   * User can create finite elements using directly core interface and
   * and them to simple problem by this function
   *
   * @return std::vector<std::string>&
   */
  inline std::vector<std::string> &getOtherFiniteElements() { return otherFEs; }
 
  /**
   * @brief Delete dm
   *
   * Destroys the PETSc DM object and cleans up resources
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode deleteDM();
 
  /**
   * @brief Delete finite elements
   *
   * Removes all finite elements created by the Simple interface
   * 
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode deleteFiniteElements();
 
  /**
   * @brief Get the addSkeletonFE flag
   *
   * If variable set to true, skeleton element is created regardless of whether
   * a field on skeleton is added or not.
   *
   * @return bool& Reference to flag indicating if skeleton FE should be created
   */
  bool &getAddSkeletonFE() { return addSkeletonFE; }
 
  /**
   * @brief Get the addBoundaryFE flag
   *
   * If variable set to true, boundary element is created regardless of whether
   * a field on boundary is added or not.
   *
   * @return bool& Reference to flag indicating if boundary FE should be created
   */
  bool &getAddBoundaryFE() { return addBoundaryFE; }
 
  /**
   * @brief Get the addParentAdjacencies flag
   *
   * If set true, parent adjacencies are added to finite elements
   *
   * @return bool& Reference to flag indicating if parent adjacencies should be added
   */
  bool &getParentAdjacencies() { return addParentAdjacencies; }
 
  /**
   * @brief Get bit ref level for parent
   *
   * @return auto& Reference to bit refinement level for parent adjacencies
   */
  auto &getBitAdjParent() { return bitAdjParent; }
 
  /**
   * @brief Get bit ref level mask for parent
   *
   * @return auto& Reference to bit refinement level mask for parent adjacencies
   */
  auto &getBitAdjParentMask() { return bitAdjParentMask; }
 
  /**
   * @brief Get bit ref level for entities
   *
   * @return auto& Reference to bit refinement level for entity adjacencies
   */
  auto &getBitAdjEnt() { return bitAdjEnt; }
 
  /**
   * @brief Get bit ref level mask for entities
   *
   * @return auto& Reference to bit refinement level mask for entity adjacencies
   */
  auto &getBitAdjEntMask() { return bitAdjEntMask; }
 
  /**
   * @brief Add empty block to problem
   *
   * MatrixManager assumes that all blocks, i.e. all field combinations are non-zero.
   * This is not always the case, to optimize code and reduce memory usage
   * user can specify which blocks are empty.
   *
   * @param row_field Name of the row field
   * @param col_field Name of the column field  
   * @return MoFEMErrorCode Error code indicating success or failure
   */
  MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string row_field,
                                            const std::string col_field) const;
 
private:
  MoFEM::Core &cOre;
 
  BitRefLevel bitLevel;     ///< BitRefLevel of the problem
  BitRefLevel bitLevelMask; ///< BitRefLevel of the problem
 
  PetscLogEvent MOFEM_EVENT_SimpleSetUP;
  PetscLogEvent MOFEM_EVENT_SimpleLoadMesh;
  PetscLogEvent MOFEM_EVENT_SimpleBuildFields;
  PetscLogEvent MOFEM_EVENT_SimpleBuildFiniteElements;
  PetscLogEvent MOFEM_EVENT_SimpleBuildProblem;
  PetscLogEvent MOFEM_EVENT_SimpleKSPSolve;
 
  EntityHandle meshSet;         ///< domain meshset
  EntityHandle boundaryMeshset; ///< meshset with boundary
  EntityHandle skeletonMeshset; ///< skeleton meshset with boundary
 
  bool addSkeletonFE;        ///< Add skeleton FE
  bool addBoundaryFE;        ///< Add boundary FE
  bool addParentAdjacencies; ///< If set to true parent adjacencies are build
 
  BitRefLevel bitAdjParent;     ///< bit ref level for parent
  BitRefLevel bitAdjParentMask; ///< bit ref level for parent parent
  BitRefLevel bitAdjEnt;        ///< bit ref level for parent
  BitRefLevel bitAdjEntMask;    ///< bit ref level for parent parent
 
  std::vector<std::string> domainFields;      ///< domain fields
  std::vector<std::string> boundaryFields;    ///< boundary fields
  std::vector<std::string> skeletonFields;    ///< fields on the skeleton
  std::vector<std::string> dataFields;        ///< Data fields
  std::vector<std::string> meshsetFields;     ///< meshset fields
  std::vector<std::string> noFieldFields;     ///< NOFIELD field name
  std::vector<std::string> noFieldDataFields; ///< NOFIELD field name
 
  std::list<std::tuple<std::string, int, Range, bool>>
      fieldsOrder; ///< fields order. 1: field name, order, range, set by range
                   ///< if true
 
  std::string nameOfProblem = "SimpleProblem"; ///< problem name
  std::string domainFE = "dFE";                ///< domain finite element
  std::string boundaryFE = "bFE";              ///< boundary finite element
  std::string skeletonFE = "sFE";              ///< skeleton finite element
  std::string meshsetFE = "mFE";               ///< meshset finite element
 
  std::vector<std::string> otherFEs; ///< Other finite elements
 
  std::vector<Range> meshsetFiniteElementEntities; ///< Meshset element entities
 
  char meshFileName[255]; ///< mesh file name
  int dIm;                ///< dimension of problem
 
  SmartPetscObj<DM>
      dM; ///< Discrete manager (interface to PETSc/MoFEM functions)
 
  MoFEMErrorCode createBoundaryMeshset();
  MoFEMErrorCode createSkeletonMeshset();
  MoFEMErrorCode exchangeGhostCells();
 
  template <int DIM = -1>
  MoFEMErrorCode setSkeletonAdjacency(std::string fe_name);
 
  template <int DIM = -1> MoFEMErrorCode setParentAdjacency();
 
  boost::shared_ptr<ParentFiniteElementAdjacencyFunction<3>>
      parentAdjFunctionDim3;
  boost::shared_ptr<ParentFiniteElementAdjacencyFunction<2>>
      parentAdjFunctionDim2;
  boost::shared_ptr<ParentFiniteElementAdjacencyFunction<1>>
      parentAdjFunctionDim1;
 
  boost::shared_ptr<ParentFiniteElementAdjacencyFunctionSkeleton<2>>
      parentAdjSkeletonFunctionDim2;
  boost::shared_ptr<ParentFiniteElementAdjacencyFunctionSkeleton<1>>
      parentAdjSkeletonFunctionDim1;
};
 
} // namespace MoFEM
 
#endif // __SIMPLE_HPP__
 
/**
 * \defgroup mofem_simple_interface Simple interface
 * \brief Implementation of simple interface for fast problem set-up.
 *
 * \ingroup mofem
 **/