#include <MoFEM.hpp>
#include <PoissonOperators.hpp>
#include <AuxPoissonFunctions.hpp>
Classes
struct	ExactFunction
	Function. More...

struct	ExactFunctionGrad
	Exact gradient. More...

struct	ExactLaplacianFunction
	Laplacian of function. More...

struct	OpS
Functions
int	main (int argc, char *argv[])
Variables
static char	help [] = "...\n\n"

static const bool	debug = false
Function Documentation

◆ main()

int main	(	int	argc,
		char *	argv[]
	)
< Volume element for the matrix
< Boundary element for the matrix
< Volume element to assemble vector
< Volume element evaluate error
< Volume element to Post-process results
< Null element do nothing
Examples: analytical_poisson_field_split.cpp.
Definition at line 207 of file analytical_poisson_field_split.cpp.
                                  {
  
   // Initialize PETSc
   MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
   // Create MoAB database
   moab::Core moab_core;              // create database
   moab::Interface &moab = moab_core; // create interface to database
  
   try {
  
     // Get command line options
     int order = 3;                    // default approximation order
     PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
     CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
                              "none");
     // Set approximation order
     CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
                            PETSC_NULL);
     // Set testing (used by CTest)
     CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
                             &flg_test, PETSC_NULL);
     ierr = PetscOptionsEnd();
     CHKERRG(ierr);
  
     // Create MoFEM database and link it to MoAB
     MoFEM::Core mofem_core(moab);           // create database
     MoFEM::Interface &m_field = mofem_core; // create interface to database
     // Register DM Manager
     CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
  
     // Create vector to store approximation global error
     Vec global_error;
     CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
  
     // First we crate elements, implementation of elements is problem
     // independent, we do not know yet what fields are present in the problem,
     // or even we do not decided yet what approximation base or spaces we are
     // going to use. Implementation of element is free from those constrains and
     // can be used in different context.
  
     // Elements used by KSP & DM to assemble system of equations
     boost::shared_ptr<ForcesAndSourcesCore>
         domain_lhs_fe; ///< Volume element for the matrix
     boost::shared_ptr<ForcesAndSourcesCore>
         boundary_lhs_fe; ///< Boundary element for the matrix
     boost::shared_ptr<ForcesAndSourcesCore>
         domain_rhs_fe; ///< Volume element to assemble vector
     boost::shared_ptr<ForcesAndSourcesCore>
         boundary_rhs_fe; ///< Volume element to assemble vector
     boost::shared_ptr<ForcesAndSourcesCore>
         domain_error; ///< Volume element evaluate error
     boost::shared_ptr<PoissonExample::PostProcFE>
         post_proc_volume; ///< Volume element to Post-process results
     boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
     boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
     {
       // Add problem specific operators the generic finite elements to calculate
       // matrices and vectors.
       CHKERR PoissonExample::CreateFiniteElements(m_field)
           .createFEToAssembleMatrixAndVector(
               ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
               boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
       // Add problem specific operators the generic finite elements to calculate
       // error on elements and global error in H1 norm
       CHKERR PoissonExample::CreateFiniteElements(m_field)
           .createFEToEvaluateError(ExactFunction(), ExactFunctionGrad(),
                                    global_error, domain_error);
       // Post-process results
       CHKERR PoissonExample::CreateFiniteElements(m_field)
           .creatFEToPostProcessResults(post_proc_volume);
  
       const double beta = 1;
       boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
           new FaceElementForcesAndSourcesCore(m_field));
       boundary_penalty_lhs_fe->getRuleHook = PoissonExample::FaceRule();
       boundary_penalty_lhs_fe->getOpPtrVector().push_back(new OpS(beta));
       boundary_rhs_fe->getOpPtrVector().push_back(
           new PoissonExample::Op_g(ExactFunction(), "U", beta));
     }
  
     // Get simple interface is simplified version enabling quick and
     // easy construction of problem.
     Simple *simple_interface;
     // Query interface and get pointer to Simple interface
     CHKERR m_field.getInterface(simple_interface);
  
     // Build problem with simple interface
     {
  
       // Get options for simple interface from command line
       CHKERR simple_interface->getOptions();
       // Load mesh file to database
       CHKERR simple_interface->loadFile();
  
       // Add field on domain and boundary. Field is declared by space and base
       // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
       // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
       // number of coefficients for dof.
       //
       // Simple interface assumes that there are four types of field; 1) defined
       // on body domain, 2) fields defined on body boundary, 3) skeleton field
       // defined on finite element skeleton. Finally data field is defined on
       // body domain. Data field is not solved but used for post-process or to
       // keep material parameters, etc. Here we using it to calculate
       // approximation error on elements.
  
       // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
       // used to construct base functions.
       CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
                                               1);
       // Add Lagrange multiplier field on body boundary
       CHKERR simple_interface->addBoundaryField("L", H1,
                                                 AINSWORTH_LEGENDRE_BASE, 1);
       // Add error (data) field, we need only L2 norm. Later order is set to 0,
       // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
       // element. You can use more DOFs and collate moments of error to drive
       // anisotropic h/p-adaptivity, however this is beyond this example.
       CHKERR simple_interface->addDataField("ERROR", L2,
                                             AINSWORTH_LEGENDRE_BASE, 1);
  
       // Set fields order domain and boundary fields.
       CHKERR simple_interface->setFieldOrder("U",
                                              order); // to approximate function
       CHKERR simple_interface->setFieldOrder("L",
                                              order); // to Lagrange multipliers
       CHKERR simple_interface->setFieldOrder(
           "ERROR", 0); // approximation order for error
  
       // Setup problem. At that point database is constructed, i.e. fields,
       // finite elements and problem data structures with local and global
       // indexing.
       CHKERR simple_interface->setUp();
     }
  
     // Get access to PETSC-MoFEM DM manager.
     // or more derails see
     // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
     // Form that point internal MoFEM data structures are linked with PETSc
     // interface. If DM functions contains string MoFEM is is MoFEM specific DM
     // interface function, otherwise DM function part of standard PETSc
     // interface.
  
     DM dm;
     // Get dm
     CHKERR simple_interface->getDM(&dm);
  
     // Solve problem, only PETEc interface here.
     {
  
       // Create the right hand side vector and vector of unknowns
       Vec F, D;
       CHKERR DMCreateGlobalVector(dm, &F);
       // Create unknown vector by creating duplicate copy of F vector. only
       // structure is duplicated no values.
       CHKERR VecDuplicate(F, &D);
  
       DM dm_sub_KK, dm_sub_LU;
       ublas::matrix<Mat> nested_matrices(2, 2);
       ublas::vector<IS> nested_is(2);
  
       CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
       CHKERR DMSetType(dm_sub_KK, "DMMOFEM");
       CHKERR DMMoFEMCreateSubDM(dm_sub_KK, dm, "SUB_KK");
       CHKERR DMSetFromOptions(dm_sub_KK);
       CHKERR DMMoFEMSetSquareProblem(dm_sub_KK, PETSC_TRUE);
       CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
       CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
       CHKERR DMMoFEMAddElement(dm_sub_KK,
                                simple_interface->getDomainFEName().c_str());
       CHKERR DMMoFEMAddElement(dm_sub_KK,
                                simple_interface->getBoundaryFEName().c_str());
       CHKERR DMSetUp(dm_sub_KK);
       CHKERR DMMoFEMGetSubRowIS(dm_sub_KK, &nested_is[0]);
       CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
       domain_lhs_fe->ksp_B = nested_matrices(0, 0);
       CHKERR DMoFEMLoopFiniteElements(
           dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
       boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
       CHKERR DMoFEMLoopFiniteElements(dm_sub_KK,
                                       simple_interface->getBoundaryFEName(),
                                       boundary_penalty_lhs_fe);
       CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
       CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
       CHKERR DMDestroy(&dm_sub_KK);
       //
       CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
       CHKERR DMSetType(dm_sub_LU, "DMMOFEM");
       CHKERR DMMoFEMCreateSubDM(dm_sub_LU, dm, "SUB_LU");
       CHKERR DMSetFromOptions(dm_sub_LU);
       CHKERR DMMoFEMSetSquareProblem(dm_sub_LU, PETSC_FALSE);
       CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
       CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
       CHKERR DMMoFEMAddElement(dm_sub_LU,
                                simple_interface->getBoundaryFEName().c_str());
       CHKERR DMSetUp(dm_sub_LU);
       CHKERR DMMoFEMGetSubRowIS(dm_sub_LU, &nested_is[1]);
       CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
       boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
       CHKERR DMoFEMLoopFiniteElements(
           dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
       CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
       CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
       // CHKERR MatCreateTranspose(nested_matrices(1,0),&nested_matrices(0,1));
       CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
                           &nested_matrices(0, 1));
       CHKERR DMDestroy(&dm_sub_LU);
  
       domain_rhs_fe->ksp_f = F;
       CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                       domain_rhs_fe);
       boundary_rhs_fe->ksp_f = F;
       CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
                                       boundary_rhs_fe);
       CHKERR VecAssemblyBegin(F);
       CHKERR VecAssemblyEnd(F);
  
       Mat A;
       nested_matrices(1, 1) = PETSC_NULL;
  
       if (debug) {
         MatType type;
         MatGetType(nested_matrices(0, 0), &type);
         cerr << "K " << type << endl;
         MatGetType(nested_matrices(1, 0), &type);
         cerr << "C " << type << endl;
         MatGetType(nested_matrices(0, 1), &type);
         cerr << "CT " << type << endl;
         std::string wait;
         cerr << "UU" << endl;
         MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
         std::cin >> wait;
         cerr << "LU" << endl;
         MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
         std::cin >> wait;
         cerr << "UL" << endl;
         MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
         std::cin >> wait;
       }
  
       CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
                            &nested_matrices(0, 0), &A);
  
       // Create solver and link it to DM
       KSP solver;
       CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
       CHKERR KSPSetFromOptions(solver);
       // Set operators
       CHKERR KSPSetOperators(solver, A, A);
       PC pc;
       CHKERR KSPGetPC(solver, &pc);
       PetscBool is_pcfs = PETSC_FALSE;
       PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
       if (is_pcfs) {
         CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
         CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
       } else {
         SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                 "Only works with pre-conditioner PCFIELDSPLIT");
       }
       // Set-up solver, is type of solver and pre-conditioners
       CHKERR KSPSetUp(solver);
       // At solution process, KSP solver using DM creates matrices, Calculate
       // values of the left hand side and the right hand side vector. then
       // solves system of equations. Results are stored in vector D.
       CHKERR KSPSolve(solver, F, D);
  
       // Scatter solution on the mesh. Stores unknown vector on field on the
       // mesh.
       CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
  
       // Clean data. Solver and vector are not needed any more.
       CHKERR KSPDestroy(&solver);
       for (int i = 0; i != 2; i++) {
         CHKERR ISDestroy(&nested_is[i]);
         for (int j = 0; j != 2; j++) {
           CHKERR MatDestroy(&nested_matrices(i, j));
         }
       }
       CHKERR MatDestroy(&A);
       CHKERR VecDestroy(&D);
       CHKERR VecDestroy(&F);
     }
  
     // Calculate error
     {
       // Loop over all elements in mesh, and run users operators on each
       // element.
       CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                       domain_error);
       CHKERR PoissonExample::AuxFunctions(m_field).assembleGhostVector(
           global_error);
       CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
       if (flg_test == PETSC_TRUE) {
         CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
       }
     }
  
     {
       // Loop over all elements in the mesh and for each execute
       // post_proc_volume element and operators on it.
       CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                       post_proc_volume);
       // Write results
       post_proc_volume->writeFile("out_vol.h5m");
     }
  
     // Destroy DM, no longer needed.
     CHKERR DMDestroy(&dm);
  
     // Destroy ghost vector
     CHKERR VecDestroy(&global_error);
   }
   CATCH_ERRORS;
  
   // finish work cleaning memory, getting statistics, etc.
   MoFEM::Core::Finalize();
  
   return 0;
 }