DMMoFEM.hpp

Go to the documentation of this file.

/** \file DMMoFEM.hpp
  \brief Discrete manager interface for MoFEM
  
  This file implements PETSc's DM (Distributed Mesh) interface for MoFEM,
  providing integration between MoFEM's Problem data structures and PETSc
  solvers (KSP, SNES, TS, TAO).
  
  \section dm_problem_connection DM-Problem Integration
  
  The DM serves as a bridge between MoFEM's core Problem structure and PETSc's
  algebraic objects (vectors and matrices). Each DM instance contains:
  - Reference to MoFEM::Interface for mesh and field access
  - Reference to MoFEM::Problem for DOF management and finite element definitions
  - Context structures for solver-specific data (KspCtx, SnesCtx, TsCtx)
  
  \section dm_solver_usage DM Usage with PETSc Solvers
  
  **KSP (Linear Solvers):**
  - DM provides matrix and vector creation routines
  - Integrates with MoFEM finite element assembly for stiffness matrices
  - Handles DOF mapping between local and global indexing
  
  **SNES (Nonlinear Solvers):**
  - DM manages residual and Jacobian evaluation through FEMethod callbacks
  - Coordinates finite element loops for nonlinear assembly operations
  - Provides context for Newton-based iteration data
  
  **TS (Time Steppers):**
  - DM handles time-dependent mass matrices and residual evaluation
  - Manages implicit and explicit time integration schemes
  - Coordinates with MoFEM field history and time-dependent boundary conditions
  
  \section dm_workflow Typical Workflow
  
  1. Create DM with DMMoFEMCreateMoFEM() linking to Problem
  2. Configure solver contexts (KSP/SNES/TS) via DMMoFEMSetXxxCtx()
  3. Set finite element methods for matrix/vector assembly
  4. Use standard PETSc solver interfaces with MoFEM-backed DM
  5. DM automatically handles MoFEM-specific operations during solve
  
  The DM abstraction enables seamless use of PETSc's solver ecosystem while
  maintaining access to MoFEM's mesh management and finite element capabilities.
  */
 
 
 
#ifndef __DMMOFEM_H
#define __DMMOFEM_H
 
#define DM_NO_ELEMENT "DMNONEFE"
 
namespace MoFEM {
 
/**
 * \brief Register MoFEM problem
 * \ingroup dm
 * 
 * @param sname Solver name string identifier
 * @return Error code
 */
PetscErrorCode DMRegister_MoFEM(const char sname[]);
 
/**
 * @brief Clear solver and Schur data owned by a MoFEM DM context
 *
 * This operation is idempotent. Call it before releasing an application-owned
 * DM when objects retained by the context can hold PETSc references back to
 * that DM. DMDestroy_MoFEM also calls it once PETSc reaches the type-specific
 * destroy callback.
 *
 * @param dm Distributed mesh manager
 * @return Error code
 */
PetscErrorCode DMMoFEMClearDMCtx(DM dm);
 
/**
 * \brief Must be called by user to set MoFEM data structures
 *
 * \note If problem exist function create DM for it. If you set bit levels,
 * those bits are to existing bits. Thus if you do not like to change bit ref
 * level for existing problem, set bits to zero.
 *
 * \ingroup dm
 */
PetscErrorCode DMMoFEMCreateMoFEM(
    DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[],
    const MoFEM::BitRefLevel bit_level,
    const MoFEM::BitRefLevel bit_mask = MoFEM::BitRefLevel().set());
 
/**
 * @brief Duplicate internal data structure
 * 
 * @param dm 
 * @param dm_duplicate 
 * @return PetscErrorCode 
 */
PetscErrorCode DMMoFEMDuplicateDMCtx(DM dm, DM dm_duplicate);
 
/**
 * @brief Swap internal data structure
 *
 * @param dm
 * @param dm_swap
 * @return PetscErrorCode
 */
PetscErrorCode DMMoFEMSwapDMCtx(DM dm, DM dm_swap);
 
/**
 * \brief Must be called by user to set Sub DM MoFEM data structures
 * \ingroup dm
 * 
 * @param subdm Sub distributed mesh manager to create
 * @param dm Parent distributed mesh manager
 * @param problem_name Name of the problem to create sub DM for
 * @return Error code
 */
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[]);
 
/**
 * \brief Get pointer to MoFEM::Interface
 * @param  dm          Distributed mesh manager
 * @param  m_field_ptr Pointer to pointer of field interface
 * @return             Error code
 * \ingroup dm
 */
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr);
 
/**
 * \brief Get pointer to problem data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param problem_ptr Pointer to pointer of problem data structure
 * @return Error code
 */
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr);
 
/**
 * If this is set to PETSC_TRUE problem is deleted with DM
 * @param  dm             the DM object
 * @param  destroy        if PETSC_TRUE problem is destroyed
 * @return                error code
 */
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem);
 
/**
 * Get if problem will be destroyed with DM
 * @param  dm             the DM object
 * @param  destroy        return if PETSC_TRUE problem is destroyed
 * @return                error code
 */
PetscErrorCode DMMoFEMGetDestroyProblem(DM dm, PetscBool *destroy_problem);
 
/**
  * \brief set squared problem
  * \ingroup dm
 
  If true, it is assumed that matrix has the same indexing on rows and
  columns. This reduces interprocessor communication.
 
  */
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem);
 
/**
  * \brief get squared problem
  * \ingroup dm
 
  If true, it is assumed that matrix has the same indexing on rows and
  columns. This reduces interprocessor communication.
 
  */
PetscErrorCode DMMoFEMGetSquareProblem(DM dm, PetscBool *square_problem);
 
/**
 * \brief Resolve shared entities
 *
 * @param  dm      dm
 * @param  fe_name finite element for which shared entities are resolved
 * @return         error code
 *
 * \note This function is valid for parallel algebra and serial mesh. It
 * should be run collectively, i.e. on all processors.
 *
 * This allows for tag reduction or tag exchange, f.e.
 *
 * \code
 * CHKERR DMMoFEMResolveSharedFiniteElements(dm,"SHELL_ELEMENT");
 * Tag th;
 * CHKERR mField.get_moab().tag_get_handle("ADAPT_ORDER",th);
 * ParallelComm* pcomm =
 * ParallelComm::get_pcomm(&mField.get_moab(),MYPCOMM_INDEX);
 * // CHKERR pcomm->reduce_tags(th,MPI_SUM,prisms);
 * CHKERR pcomm->exchange_tags(th,prisms);
 * \endcode
 *
 * \ingroup dm
 */
PetscErrorCode DMMoFEMResolveSharedFiniteElements(DM dm, std::string fe_name);
 
/**
 * \brief Get finite elements layout in the problem
 *
 * In layout is stored information how many elements are on each processor, for
 * more information look at petsc documentation
 * <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscLayoutCreate.html#PetscLayoutCreate>
 *
 * @param  dm     discrete manager for this problem
 * @param  fe_name finite element name
 * @param  layout pointer to layout, for created layout user takes
 * responsibility for destroying it.
 * @return        error code
 *
 * \ingroup dm
 */
PetscErrorCode DMMoFEMGetProblemFiniteElementLayout(DM dm, std::string fe_name,
                                                    PetscLayout *layout);
 
/**
 * \brief add element to dm
 * \ingroup dm
 *
 * \note add_file is a collective, should be executed on all processors.
 * Otherwise could lead to deadlock.
 *
 * @param dm Distributed mesh manager
 * @param fe_name Finite element name
 * @return Error code
 */
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name);
 
/**
 * \brief add element to dm
 * \ingroup dm
 *
 * \note add_file is a collective, should be executed on all processors.
 * Otherwise could lead to deadlock.
 *
 * @param dm Distributed mesh manager
 * @param fe_name Vector of finite element names
 * @return Error code
 */
PetscErrorCode DMMoFEMAddElement(DM dm, std::vector<std::string> fe_name);
 
/**
 * \brief unset element from dm
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param fe_name Finite element name to remove
 * @return Error code
 */
PetscErrorCode DMMoFEMUnSetElement(DM dm, std::string fe_name);
 
/**
  * \brief set local (or ghosted) vector values on mesh for partition only
  * \ingroup dm
 
  * \param l vector
  * \param mode see petsc manual for VecSetValue (ADD_VALUES or INSERT_VALUES)
  * \param scatter_mode see petsc manual for ScatterMode (The available modes
  are: SCATTER_FORWARD or SCATTER_REVERSE)
  *
  * SCATTER_REVERSE set data to field entities from V vector.
  * SCATTER_FORWARD set vector V from data field entities
  * 
  * This code will update ghosted entities on all processors
   \code 
  * auto d_vector =  createDMVector(dm);
  * // fill d_vector on owned entities
    CHKERR DMoFEMMeshToLocalVector(dm, d_vector, INSERT_VALUES,
                                       SCATTER_FORWARD);
    // update ghost values in d_vector
    CHKERR VecGhostUpdateBegin(d_vector, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d_vector, INSERT_VALUES, SCATTER_FORWARD);
    // now d_vector has updated ghost values, and following save values of 
    // vector to fields entities
    CHKERR DMoFEMMeshToLocalVector(dm, d_vector, INSERT_VALUES,
                                       SCATTER_REVERSE);
  * \endcode
  *
  * \note you use DMoFEMMeshToGlobalVector as scatter reverse
  *
  * For implementations see VecManager::setLocalGhostVector.
  */
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode,
                                       ScatterMode scatter_mode,
                                       RowColData rc = RowColData::COL);
 
/**
  * \brief set ghosted vector values on all existing mesh entities
  * \ingroup dm
 
  * \param g vector
  * \param mode see petsc manual for VecSetValue (ADD_VALUES or INSERT_VALUES)
  * \param scatter_mode see petsc manual for ScatterMode (The available modes
  are: SCATTER_FORWARD or SCATTER_REVERSE)
  *
  * SCATTER_REVERSE set data to field entities from V vector.
  *
  * SCATTER_FORWARD set vector V from data field entities
 
  */
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode,
                                        ScatterMode scatter_mode,
                                        RowColData rc = RowColData::COL);
 
/**
 * \brief execute finite element method for each element in dm (problem)
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param method Finite element method to execute for preprocessing
 * @return Error code
 */
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method);
 
/**
 * \brief execute finite element method for each element in dm (problem)
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param method Finite element method to execute for postprocessing
 * @return Error code
 */
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method);
 
/**
 * \brief Executes FEMethod for finite elements in DM
 * @param  dm       MoFEM discrete manager
 * @param  fe_name  name of finite element
 * @param  method   pointer to \ref MoFEM::FEMethod
 * @param  low_rank lowest rank of processor
 * @param  up_rank  upper run of processor
 * @return          Error code
 * \ingroup dm
 */
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(
    DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank,
    int up_rank, CacheTupleWeakPtr cache_ptr = CacheTupleSharedPtr());
 
/**
 * \brief Executes FEMethod for finite elements in DM
 * @param  dm       MoFEM discrete manager
 * @param  fe_name  name of finite element
 * @param  method   pointer to \ref MoFEM::FEMethod
 * @param  low_rank lowest rank of processor
 * @param  up_rank  upper run of processor
 * @return          Error code
 * \ingroup dm
 */
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(
    DM dm, const std::string fe_name, boost::shared_ptr<MoFEM::FEMethod> method,
    int low_rank, int up_rank,
    CacheTupleWeakPtr cache_ptr = CacheTupleSharedPtr());
 
/**
 * \brief Executes FEMethod for finite elements in DM
 * @param  dm      MoFEM discrete manager
 * @param  fe_name name of element
 * @param  method  pointer to \ref MOFEM::FEMethod
 * @return         Error code
 * \ingroup dm
 */
PetscErrorCode
DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method,
                         CacheTupleWeakPtr cache_ptr = CacheTupleSharedPtr());
 
/**
 * \brief Executes FEMethod for finite elements in DM
 * @param  dm      MoFEM discrete manager
 * @param  fe_name name of element
 * @param  method  pointer to \ref MOFEM::FEMethod
 * @return         Error code
 * \ingroup dm
 */
PetscErrorCode
DMoFEMLoopFiniteElements(DM dm, const std::string fe_name,
                         boost::shared_ptr<MoFEM::FEMethod> method,
                         CacheTupleWeakPtr cache_ptr = CacheTupleSharedPtr());
 
/**
 * \brief execute method for dofs on field in problem
 * \ingroup dm
 */
PetscErrorCode DMoFEMLoopDofs(DM dm, const char field_name[],
                              MoFEM::DofMethod *method);
 
// /**
//  * \brief Set compute operator for KSP solver via sub-matrix and IS
//  *
//  * @param  dm   DM
//  * @return      error code
//  *
//  * \ingroup dm
//  */
// PetscErrorCode DMMoFEMKSPSetComputeOperatorsViaSubMatrixbByIs(DM dm);
 
/**
 * \brief set KSP right hand side evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[],
                                       MoFEM::FEMethod *method,
                                       MoFEM::BasicMethod *pre_only,
                                       MoFEM::BasicMethod *post_only);
 
/**
 * \brief set KSP right hand side evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMKSPSetComputeRHS(DM dm, const std::string fe_name,
                        boost::shared_ptr<MoFEM::FEMethod> method,
                        boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                        boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief Set KSP operators and push mofem finite element methods
 *
 * @param  dm        DM
 * @param  fe_name   finite element name
 * @param  method    method on the element (executed for each element in the
 * problem which given name)
 * @param  pre_only  method for pre-process before element method
 * @param  post_only method for post-process after element method
 * @return           error code
 *
 * \ingroup dm
 */
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[],
                                             MoFEM::FEMethod *method,
                                             MoFEM::BasicMethod *pre_only,
                                             MoFEM::BasicMethod *post_only);
 
/**
 * \brief Set KSP operators and push mofem finite element methods
 *
 * @param  dm        DM
 * @param  fe_name   finite element name
 * @param  method    method on the element (executed for each element in the
 * problem which given name)
 * @param  pre_only  method for pre-process before element method
 * @param  post_only method for post-process after element method
 * @return           error code
 *
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMKSPSetComputeOperators(DM dm, const std::string fe_name,
                              boost::shared_ptr<MoFEM::FEMethod> method,
                              boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                              boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief set SNES residual evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[],
                                      MoFEM::FEMethod *method,
                                      MoFEM::BasicMethod *pre_only,
                                      MoFEM::BasicMethod *post_only);
 
/**
 * \brief set SNES residual evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMSNESSetFunction(DM dm, const std::string fe_name,
                       boost::shared_ptr<MoFEM::FEMethod> method,
                       boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                       boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief set SNES Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[],
                                      MoFEM::FEMethod *method,
                                      MoFEM::BasicMethod *pre_only,
                                      MoFEM::BasicMethod *post_only);
 
/**
 * \brief set SNES Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMSNESSetJacobian(DM dm, const std::string fe_name,
                       boost::shared_ptr<MoFEM::FEMethod> method,
                       boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                       boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief set TS implicit function evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[],
                                     MoFEM::FEMethod *method,
                                     MoFEM::BasicMethod *pre_only,
                                     MoFEM::BasicMethod *post_only);
 
/**
 * \brief set TS implicit function evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMTSSetIFunction(DM dm, const std::string fe_name,
                      boost::shared_ptr<MoFEM::FEMethod> method,
                      boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                      boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief set TS Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name,
                      boost::shared_ptr<MoFEM::FEMethod> method,
                      boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                      boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * \brief set TS Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const char fe_name[],
                                     MoFEM::FEMethod *method,
                                     MoFEM::BasicMethod *pre_only,
                                     MoFEM::BasicMethod *post_only);
 
/**
 * @brief set TS the right hand side function
 *
 * <a
 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSSetRHSFunction.html#TSSetRHSFunction>See
 * petsc documentation</a>
 *
 * @param dm
 * @param fe_name
 * @param method
 * @param pre_only
 * @param post_only
 * @return PetscErrorCode
 */
PetscErrorCode
DMMoFEMTSSetRHSFunction(DM dm, const std::string fe_name,
                        boost::shared_ptr<MoFEM::FEMethod> method,
                        boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                        boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, const char fe_name[],
                                       MoFEM::FEMethod *method,
                                       MoFEM::BasicMethod *pre_only,
                                       MoFEM::BasicMethod *post_only);
 
/**
 * @brief set TS the right hand side jacobian
 *
 * <a
 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSSetRHSJacobian.html>See
 * petsc documentation</a>
 *
 * @param dm
 * @param fe_name
 * @param method
 * @param pre_only
 * @param post_only
 * @return PetscErrorCode
 */
PetscErrorCode
DMMoFEMTSSetRHSJacobian(DM dm, const std::string fe_name,
                        boost::shared_ptr<MoFEM::FEMethod> method,
                        boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                        boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
PetscErrorCode DMMoFEMTSSetRHSJacobian(DM dm, const char fe_name[],
                                       MoFEM::FEMethod *method,
                                       MoFEM::BasicMethod *pre_only,
                                       MoFEM::BasicMethod *post_only);
 
/**
 * \brief set TS implicit function evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMTSSetI2Function(DM dm, const std::string fe_name,
                       boost::shared_ptr<MoFEM::FEMethod> method,
                       boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                       boost::shared_ptr<MoFEM::BasicMethod> post_only);
/**
 * \brief set TS implicit function evaluation function
 * \ingroup dm
 */
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const char fe_name[],
                                      MoFEM::FEMethod *method,
                                      MoFEM::BasicMethod *pre_only,
                                      MoFEM::BasicMethod *post_only);
 
/**
 * \brief set TS Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name,
                       boost::shared_ptr<MoFEM::FEMethod> method,
                       boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                       boost::shared_ptr<MoFEM::BasicMethod> post_only);
/**
 * \brief set TS Jacobian evaluation function
 * \ingroup dm
 */
PetscErrorCode
DMMoFEMTSSetI2Jacobian(DM dm, const char fe_name[],
                                     MoFEM::FEMethod *method,
                                     MoFEM::BasicMethod *pre_only,
                                     MoFEM::BasicMethod *post_only);
 
/**
 * @brief Set Monitor To TS solver
 *
 * <a
 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSMonitorSet.html>See
 * PETSc documentation here</a>
 *
 * @param dm
 * @param ts time solver
 * @param fe_name
 * @param method
 * @param pre_only
 * @param post_only
 * @return PetscErrorCode
 */
PetscErrorCode
DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name,
                    boost::shared_ptr<MoFEM::FEMethod> method,
                    boost::shared_ptr<MoFEM::BasicMethod> pre_only,
                    boost::shared_ptr<MoFEM::BasicMethod> post_only);
 
/**
 * @brief Set Monitor To TS solver
 *
 * <a
 * href=https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSMonitorSet.html>See
 * PETSc documentation here</a>
 *
 * @param dm
 * @param ts time solver
 * @param fe_name
 * @param method
 * @param pre_only
 * @param post_only
 * @return PetscErrorCode
 */
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const char fe_name[],
                                   MoFEM::FEMethod *method,
                                   MoFEM::BasicMethod *pre_only,
                                   MoFEM::BasicMethod *post_only);
 
/**
 * \brief get MoFEM::KspCtx data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param ksp_ctx Pointer to KSP context structure
 * @return Error code
 */
PetscErrorCode DMMoFEMGetKspCtx(DM dm, MoFEM::KspCtx **ksp_ctx);
 
/**
 * \brief get MoFEM::KspCtx data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param ksp_ctx Shared pointer to KSP context structure
 * @return Error code
 */
PetscErrorCode
DMMoFEMGetKspCtx(DM dm, const boost::shared_ptr<MoFEM::KspCtx> &ksp_ctx);
 
/**
 * \brief set MoFEM::KspCtx data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param ksp_ctx Shared pointer to KSP context structure to set
 * @return Error code
 */
PetscErrorCode DMMoFEMSetKspCtx(DM dm,
                                boost::shared_ptr<MoFEM::KspCtx> ksp_ctx);
 
/**
 * \brief get MoFEM::SnesCtx data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param snes_ctx Pointer to SNES context structure
 * @return Error code
 */
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx);
 
/**
 * \brief get MoFEM::SnesCtx data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param snes_ctx Shared pointer to SNES context structure
 * @return Error code
 */
PetscErrorCode
DMMoFEMGetSnesCtx(DM dm, const boost::shared_ptr<MoFEM::SnesCtx> &snes_ctx);
 
/**
  * \brief Set MoFEM::SnesCtx data structure
  * \ingroup dm
 
  */
PetscErrorCode DMMoFEMSetSnesCtx(DM dm,
                                 boost::shared_ptr<MoFEM::SnesCtx> snes_ctx);
 
/**
 * \brief get MoFEM::TsCtx data structure
 * \ingroup dm
 */
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx);
 
/**
 * \brief get MoFEM::TsCtx data structure
 * \ingroup dm
 */
PetscErrorCode DMMoFEMGetTsCtx(DM dm,
                               const boost::shared_ptr<MoFEM::TsCtx> &ts_ctx);
 
/**
  * \brief Set MoFEM::TsCtx data structure
  * \ingroup dm
 
  It take over pointer, do not delete it, DM will destroy pointer
  when is destroyed.
 
  */
PetscErrorCode DMMoFEMSetTsCtx(DM dm, boost::shared_ptr<MoFEM::TsCtx> ts_ctx);
 
/** sets if read mesh is partitioned
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param is_partitioned Boolean flag indicating if mesh is partitioned
 * @return Error code
 */
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned);
 
/** get if read mesh is partitioned
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param is_partitioned Pointer to boolean flag for partition status
 * @return Error code
 */
PetscErrorCode DMMoFEMGetIsPartitioned(DM dm, PetscBool *is_partitioned);
 
/**
 * \brief Set operators for MoFEM dm
 * @param dm Distributed mesh manager
 * @return Error code
 * \ingroup dm
 */
PetscErrorCode DMSetOperators_MoFEM(DM dm);
 
/**
 * \brief Create dm data structure with MoFEM data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager to create
 * @return Error code
 */
PetscErrorCode DMCreate_MoFEM(DM dm);
 
/**
 * \brief Destroys dm with MoFEM data structure
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager to destroy
 * @return Error code
 */
PetscErrorCode DMDestroy_MoFEM(DM dm);
 
/**
 * \brief DMShellSetCreateGlobalVector
 * \ingroup dm
 *
 * sets the routine to create a global vector
 * associated with the shell DM
 * 
 * @param dm Distributed mesh manager
 * @param g Pointer to global vector to create
 * @return Error code
 */
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g);
 
/**
 * \brief DMShellSetCreateGlobalVector
 * \ingroup dm
 *
 * sets the routine to create a global vector
 * associated with the shell DM
 */
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, SmartPetscObj<Vec> &g_ptr,
                                          RowColData rc = RowColData::COL);
 
/**
 * \brief DMShellSetCreateLocalVector
 * \ingroup dm
 *
 * sets the routine to create a local vector
 * associated with the shell DM
 */
PetscErrorCode DMCreateLocalVector_MoFEM(DM dm, Vec *l);
 
/**
 * DMShellSetCreateMatrix
 * \ingroup dm
 *
 * sets the routine to create a matrix associated with the shell DM
 */
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M);
 
/**
 * DMShellSetCreateMatrix
 * \ingroup dm
 *
 * sets the routine to create a matrix associated with the shell DM
 */
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, SmartPetscObj<Mat> &M);
 
/**
 * Set options for MoFEM DM
 * \ingroup dm
 */
#if PETSC_VERSION_GE(3, 17, 0)
PetscErrorCode DMSetFromOptions_MoFEM(DM dm,
                                      PetscOptionItems *PetscOptionsObject);
#elif PETSC_VERSION_GE(3, 7, 0)
PetscErrorCode DMSetFromOptions_MoFEM(PetscOptionItems *PetscOptionsObject,
                                      DM dm);
#elif PETSC_VERSION_GE(3, 5, 3)
PetscErrorCode DMSetFromOptions_MoFEM(PetscOptions *PetscOptionsObject, DM dm);
#else
PetscErrorCode DMSetFromOptions_MoFEM(DM dm);
#endif
 
/**
 * sets up the MoFEM structures inside a DM object
 * \ingroup dm
 */
PetscErrorCode DMSetUp_MoFEM(DM dm);
 
/**
 * Sets up the MoFEM structures inside a DM object for sub dm
 * \ingroup dm
 * 
 * @param subdm Sub distributed mesh manager to set up
 * @return Error code
 */
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm);
 
/**
 * Add field to sub dm problem on rows
 * \ingroup dm
 * 
 * @param dm Distributed mesh manager
 * @param field_name Name of field to add to rows
 * @return Error code
 */
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[]);
 
/** @copydoc DMMoFEMAddSubFieldRow(DM dm, const char field_name[]) */
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, std::string field_name);
 
/**
 * @brief Add field to sub dm problem on rows
 * \ingroup dm
 * 
 * @param dm 
 * @param field_name 
 * @param m 
 * @return PetscErrorCode 
 */
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[],
                                     boost::shared_ptr<Range> r_ptr);
 
/** @copydoc DMMoFEMAddSubFieldRow(DM dm, const char field_name[], boost::shared_ptr<Range> r_ptr) */
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, std::string field_name,
                                     boost::shared_ptr<Range> r_ptr);
 
/**
 * Add field to sub dm problem on columns
 * \ingroup dm
 */
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[]);
 
/** @copydoc DMMoFEMAddSubFieldCol(DM dm, const char field_name[]) */
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, std::string field_name);
 
/**
 * @brief Add field to sub dm problem on columns
 *
 * @param dm
 * @param field_name
 * @param range of entities
 * @return PetscErrorCode
 */
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[],
                                     boost::shared_ptr<Range> r_ptr);
 
 
/** @copydoc DMMoFEMAddSubFieldCol(DM dm, const char field_name[], boost::shared_ptr<Range> r_ptr); */
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, std::string field_name,
                                     boost::shared_ptr<Range> r_ptr);
 
/**
 * Return true if this DM is sub problem
  \ingroup dm
 * @param  dm            the DM object
 * @param  is_subproblem true if subproblem
 * @return               error code
 */
PetscErrorCode DMMoFEMGetIsSubDM(DM dm, PetscBool *is_sub_dm);
 
/**
 * \brief get sub problem is
 * @param  dm has to be created with DMMoFEMSetSquareProblem
 * @param  is return is on the row
 * @return    error code
 *
 * Returns IS with global indices of the DM used to create SubDM
 *
 */
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is);
 
/**
 * \brief get sub problem is
 * @param  dm has to be created with DMMoFEMSetSquareProblem
 * @param  is return is on the row
 * @return    error code
 *
 * Returns IS with global indices of the DM used to create SubDM
 *
 */
PetscErrorCode DMMoFEMGetSubColIS(DM dm, IS *is);
 
/**
 * \brief Add problem to composite DM on row
  \ingroup dm
 *
 * This create block on row with DOFs from problem of given name
 *
 * @param  dm       the DM object
 * @param  prb_name add problem name
 * @return          error code
 */
PetscErrorCode DMMoFEMAddRowCompositeProblem(DM dm, const char prb_name[]);
 
/**
 * \brief Add problem to composite DM on col
 * \ingroup dm
 *
 * This create block on col with DOFs from problem of given name
 *
 * @param  dm       the DM object
 * @param  prb_name add problem name
 * @return          error code
 */
PetscErrorCode DMMoFEMAddColCompositeProblem(DM dm, const char prb_name[]);
 
/**
 * \brief Get if this DM is composite DM
 * \ingroup dm
 *
 * @param  dm         the DM object
 * @param  is_comp_dm return true if composite problem here
 * @return            error code
 */
PetscErrorCode DMMoFEMGetIsCompDM(DM dm, PetscBool *is_comp_dm);
 
/**
 * destroy the MoFEM structure
 * \ingroup dm
 */
PetscErrorCode DMDestroy_MoFEM(DM dm);
 
/**
 * DMShellSetGlobalToLocal
 * \ingroup dm
 *
 * the routine that begins the global to local scatter
 */
PetscErrorCode DMGlobalToLocalBegin_MoFEM(DM dm, Vec, InsertMode, Vec);
 
/**
 * DMShellSetGlobalToLocal
 * \ingroup dm
 *
 * the routine that begins the global to local scatter
 */
PetscErrorCode DMGlobalToLocalEnd_MoFEM(DM dm, Vec, InsertMode, Vec);
 
/**
 * DMShellSetLocalToGlobal
 * \ingroup dm
 *
 * the routine that begins the local to global scatter
 */
PetscErrorCode DMLocalToGlobalBegin_MoFEM(DM, Vec, InsertMode, Vec);
 
/**
 * DMShellSetLocalToGlobal
 * \ingroup dm
 *
 * the routine that ends the local to global scatter
 */
PetscErrorCode DMLocalToGlobalBegin_MoFEM(DM, Vec, InsertMode, Vec);
 
/**
 * DMShellSetLocalToGlobal
 * \ingroup dm
 *
 * the routine that ends the local to global scatter
 */
PetscErrorCode DMLocalToGlobalEnd_MoFEM(DM, Vec, InsertMode, Vec);
 
/**
 * Creates a set of IS objects with the global indices of dofs for each field
 * @param  dm         The number of fields (or NULL if not requested)
 *
 * Output:
 * @param  numFields  The number of fields (or NULL if not requested)
 * @param  fieldNames The name for each field (or NULL if not requested)
 * @param  fields     The global indices for each field (or NULL if not
 requested)
 *
 * @return            error code
 
 * \note The user is responsible for freeing all requested arrays. In
 particular,
 * every entry of names should be freed with PetscFree(), every entry of fields
 * should be destroyed with ISDestroy(), and both arrays should be freed with
 * PetscFree().
 
  \ingroup dm
 
 */
PetscErrorCode DMCreateFieldIS_MoFEM(DM dm, PetscInt *numFields,
                                     char ***fieldNames, IS **fields);
 
/**
 * \brief get field is in the problem
 * @param  dm         the DM objects
 * @param  rc         ROW or COL (no difference is problem is squared)
 * @param  field_name name of the field
 * @param  is         returned the IS object
 * @return            error code
 *
 * \code
 * IS is;
 * ierr = DMMoFEMGetFieldIS(dm,ROW,"DISP",&is_disp); CHKERRG(ierr);
 * \endcode
 *
 *
  \ingroup dm
 */
PetscErrorCode DMMoFEMGetFieldIS(DM dm, RowColData rc, const char field_name[],
                                 IS *is);
 
/**
 * @brief Set verbosity level
 *
 * @param dm
 * @param verb see VERBOSITY_LEVELS for list of the levels
 * @return PetscErrorCode
 */
PetscErrorCode DMMoFEMSetVerbosity(DM dm, const int verb);
 
struct BlockStructure;
 
/** @brief Set data for block mat
 *
 * \note You can reset data by setting nullptr
 *
 */
MoFEMErrorCode DMMoFEMSetBlocMatData(DM dm,
                                     boost::shared_ptr<BlockStructure>);
 
/**
 * @brief Get data for block mat
 * 
 * @param dm 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode DMMoFEMGetBlocMatData(DM dm, boost::shared_ptr<BlockStructure> &);
 
/**
 * @brief Create block matrix
  \ingroup dm
 * 
 * @param dm 
 * @param mat 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode DMMoFEMCreateBlockMat(DM dm, Mat *mat);
 
/**
 * @brief Create block matrix
  \ingroup dm
 *
 * @param dm
 * @param mat smart pointer
 * @return MoFEMErrorCode
 */
MoFEMErrorCode DMMoFEMCreateBlockMat(DM dm, SmartPetscObj<Mat> &mat);
 
/**
 * @brief Set data for nest schur (see specialisation in Schur.hpp)
 * 
 * \note You can reset data by setting nullptr
 * 
 * @tparam T = NestSchurData
 * @param dm 
 * @return MoFEMErrorCode 
 */
template <typename T>
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr<T>);
 
/**
 * @brief Create nest schur matrix
 * 
 * @param dm 
 * @param mat 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode DMMoFEMCreateNestSchurMat(DM dm, Mat *mat);
 
/**
 * @brief Create matrix for hybridised system
 * 
 * @param dm 
 * @param mat 
 * @return MoFEMErrorCode 
 */
MoFEMErrorCode DMMoFEMCreateHybridL2Mat(DM dm, SmartPetscObj<Mat> &mat);
 
/**
 * \brief PETSc  Discrete Manager data structure
 *
 * This structure should not be accessed or modified by user. Is not
 * available from outside MoFEM DM manager. However user can inherit dat
 * class and add data for additional functionality.
 *
 * This is part of implementation for PETSc interface, see more details in
 * <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
 *
 * \ingroup dm
 *
 */
struct DMCtx : public UnknownInterface {
 
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  virtual ~DMCtx() = default;
  
  boost::shared_ptr<KspCtx> kspCtx;   ///< data structure KSP
  boost::shared_ptr<SnesCtx> snesCtx; ///< data structure SNES
  boost::shared_ptr<TsCtx> tsCtx;     ///< data structure for TS solver
 
protected:
  DMCtx() = default;
};
 
/**
 * @brief Get the Interface Ptr object
 * 
 * @param dm 
 * @return auto 
 */
inline auto getInterfacePtr(DM dm) {
  MoFEM::Interface *m_field_ptr;
  CHK_THROW_MESSAGE(DMoFEMGetInterfacePtr(dm, &m_field_ptr),
                    "Can not get interface ptr from DM");
  return m_field_ptr;
};
 
/**
 * @brief get problem pointer from DM
 * 
 */
inline auto getProblemPtr(DM dm) {
  const MoFEM::Problem *problem_ptr;
  CHK_THROW_MESSAGE(DMMoFEMGetProblemPtr(dm, &problem_ptr),
                    "Get cot get problem ptr from DM");
  return problem_ptr;
};
 
/**
 * @brief Get smart matrix from DM
 * \ingroup dm
 * 
 */
inline auto createDMMatrix(DM dm) {
  SmartPetscObj<Mat> a;
  ierr = DMCreateMatrix_MoFEM(dm, a);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return a;
};
 
/**
 * @brief Get smart hybridised L2 matrix from DM
 *
 * @param dm
 * @return auto
 */
inline auto createDMHybridisedL2Matrix(DM dm) {
  SmartPetscObj<Mat> a;
  ierr = DMMoFEMCreateHybridL2Mat(dm, a);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return a;
};
 
inline auto createDMBlockMat(DM dm) {
  Mat raw_a;
  ierr = DMMoFEMCreateBlockMat(dm, &raw_a);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return SmartPetscObj<Mat>(raw_a);
};
 
inline auto createDMNestSchurMat(DM dm) {
  Mat raw_a;
  ierr = DMMoFEMCreateNestSchurMat(dm, &raw_a);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return SmartPetscObj<Mat>(raw_a);
};
 
 
/** @deprecated use createDMMatrix */
DEPRECATED inline auto smartCreateDMMatrix(DM dm) { return createDMMatrix(dm); }
 
/**
 * @brief Get smart vector from DM
 * \ingroup dm
 * 
 */
inline auto createDMVector(DM dm, RowColData rc = RowColData::COL) {
  SmartPetscObj<Vec> v;
  ierr = DMCreateGlobalVector_MoFEM(dm, v, rc);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return v;
};
 
/** @deprecated use createDMVector*/
DEPRECATED inline auto smartCreateDMVector(DM dm) { return createDMVector(dm); }
 
/**
 * @brief Get KSP context data structure used by DM
 * 
 */
inline auto getDMKspCtx(DM dm) {
  boost::shared_ptr<MoFEM::KspCtx> ksp_ctx;
  ierr = DMMoFEMGetKspCtx(dm, ksp_ctx);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return ksp_ctx;
};
 
/** @deprecated getDMKspCtx  */
DEPRECATED inline auto smartGetDMKspCtx(DM dm) { return getDMKspCtx(dm); }
 
/**
 * @brief Get SNES context data structure used by DM
 * 
 */
inline auto getDMSnesCtx(DM dm) {
  boost::shared_ptr<MoFEM::SnesCtx> snes_ctx;
  ierr = DMMoFEMGetSnesCtx(dm, snes_ctx);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return snes_ctx;
};
 
/** @deprecated use smartGetDMSnesCtx*/
DEPRECATED inline auto smartGetDMSnesCtx(DM dm) { return getDMSnesCtx(dm); }
 
/**
 * @brief Get TS context data structure used by DM
 * 
 */
inline auto getDMTsCtx(DM dm) {
  boost::shared_ptr<MoFEM::TsCtx> ts_ctx;
  ierr = DMMoFEMGetTsCtx(dm, ts_ctx);
  CHKERRABORT(getCommFromPetscObject(reinterpret_cast<PetscObject>(dm)), ierr);
  return ts_ctx;
};
 
/** @deprecated use getDMTsCtx */
DEPRECATED inline auto smartGetDMTsCtx(DM dm) { return getDMTsCtx(dm); }
 
/**
 * @brief  Get sub problem data structure
 *
 * @param dm
 * @return auto
 */
inline auto getDMSubData(DM dm) {
  auto prb_ptr = getProblemPtr(dm);
  return prb_ptr->getSubData();
};
 
 
} // namespace MoFEM
 
#endif //__DMMOFEM_H
 
/**
 * \defgroup dm Distributed mesh manager
 * \brief Implementation of PETSc DM, managing interactions between mesh data
 *structures and vectors and matrices
 *
 * DM objects are used to manage communication between the algebraic structures
 *in PETSc (Vec and Mat) and mesh data structures in PDE-based (or other)
 * simulations.
 *
 * DM is abstract interface, here is it particular implementation for MoFEM
 *code.
 *
 * \ingroup mofem
 **/