v0.13.1
Loading...
Searching...
No Matches
analytical_poisson.cpp

For more information and detailed explain of this example see COR-2: Solving the Poisson equation

/**
* \file analytical_poisson.cpp
* \ingroup mofem_simple_interface
* \example analytical_poisson.cpp
*
* For more information and detailed explain of this
* example see \ref poisson_tut1
*
*
*/
#include <PoissonOperators.hpp>
static char help[] = "...\n\n";
/**
* \brief Function
*
* This is prescribed exact function. If this function is given by polynomial
* order of "p" and order of approximation is "p" or higher, solution of
* finite element method is exact (with machine precision).
*
* \f[
* u = 1+x^2+y^2+z^3
* \f]
*
*/
struct ExactFunction {
double operator()(const double x, const double y, const double z) const {
return 1 + x * x + y * y + z * z * z;
}
};
/**
* \brief Exact gradient
*/
FTensor::Tensor1<double, 3> operator()(const double x, const double y,
const double z) const {
grad(0) = 2 * x;
grad(1) = 2 * y;
grad(2) = 3 * z * z;
return grad;
}
};
/**
* \brief Laplacian of function.
*
* This is Laplacian of \f$u\f$, it is calculated using formula
* \f[
* \nabla^2 u(x,y,z) = \nabla \cdot \nabla u
* \frac{\partial^2 u}{\partial x^2}+
* \frac{\partial^2 u}{\partial y^2}+
* \frac{\partial^2 u}{\partial z^2}
* \f]
*
*/
double operator()(const double x, const double y, const double z) const {
return 4 + 6 * z;
}
};
int main(int argc, char *argv[]) {
// 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
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();
// 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;
// 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.
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
global_error, domain_error);
// Post-process results
.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 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);
// 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
dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
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);
global_error);
if (flg_test == PETSC_TRUE) {
}
}
{
// 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);
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMMoFEM.cpp:584
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMMoFEM.cpp:625
double D
const FTensor::Tensor2< T, Dim, Dim > Vec
Exact gradient.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double z) const
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Create finite elements instances.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.