analytical_poisson.cpp File Reference

#include <MoFEM.hpp>
#include <PoissonOperators.hpp>
#include <AuxPoissonFunctions.hpp>

Go to the source code of this file.

Classes
struct	ExactFunction
	Function. More...
struct	ExactFunctionGrad
	Exact gradient. More...
struct	ExactLaplacianFunction
	Laplacian of function. More...

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

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

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 to assemble vector

< Volume element evaluate error

< Volume element to Post-process results

< Null element do nothing

Definition at line 70 of file analytical_poisson.cpp.

                                 {
 
  // Initialize PETSc
  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
  try {
 
    // Create MoAB database
    moab::Core moab_core;              // create database
    moab::Interface &moab = moab_core; // create interface to database
 
    // Get command line options
    int order          = 3;           // default approximation order
    PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
                             "none");
    // Set approximation order
    CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
                           PETSC_NULLPTR);
    // Set testing (used by CTest)
    CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
                            &flg_test, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    // 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
    {
      // 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);
      // 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);
    }
 
    // 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 field "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); 
 
    // Set KSP context for DM. At that point only elements are added to DM operators.
    // Calculations of matrices and vectors is executed by KSP solver. This part
    // of the code makes connection between implementation of finite elements and
    // data operators with finite element declarations in DM manager. From that
    // point DM takes responsibility for executing elements, calculations of
    // matrices and vectors, and solution of the problem.
    {
      // Set operators for KSP solver
      CHKERR DMMoFEMKSPSetComputeOperators(
          dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
      CHKERR DMMoFEMKSPSetComputeOperators(
          dm, simple_interface->getBoundaryFEName(), boundary_lhs_fe, null,
          null);
      // Set calculation of the right hand side vector for KSP solver
      CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getDomainFEName(),
                                     domain_rhs_fe, null, null);
      CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getBoundaryFEName(),
                                     boundary_rhs_fe, null, null);
    }
 
    // 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); 
 
      // Create solver and link it to DM
      KSP solver;
      CHKERR KSPCreate(PETSC_COMM_WORLD,&solver); 
      CHKERR KSPSetFromOptions(solver); 
      CHKERR KSPSetDM(solver,dm); 
      // 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); 
      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;
 
}

Variable Documentation

help

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

static

Definition at line 18 of file analytical_poisson.cpp.