v0.14.0
analytical_poisson.cpp
Go to the documentation of this file.
1 /**
2  * \file analytical_poisson.cpp
3  * \ingroup mofem_simple_interface
4  * \example analytical_poisson.cpp
5  *
6  * For more information and detailed explain of this
7  * example see \ref poisson_tut1
8  *
9  *
10  */
11 
12 
13 
14 #include <BasicFiniteElements.hpp>
15 #include <PoissonOperators.hpp>
16 #include <AuxPoissonFunctions.hpp>
17 
18 static char help[] = "...\n\n";
19 
20 /**
21  * \brief Function
22  *
23  * This is prescribed exact function. If this function is given by polynomial
24  * order of "p" and order of approximation is "p" or higher, solution of
25  * finite element method is exact (with machine precision).
26  *
27  * \f[
28  * u = 1+x^2+y^2+z^3
29  * \f]
30  *
31  */
32 struct ExactFunction {
33  double operator()(const double x, const double y, const double z) const {
34  return 1 + x * x + y * y + z * z * z;
35  }
36 };
37 
38 /**
39  * \brief Exact gradient
40  */
41 struct ExactFunctionGrad {
42  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
43  const double z) const {
45  grad(0) = 2 * x;
46  grad(1) = 2 * y;
47  grad(2) = 3 * z * z;
48  return grad;
49  }
50 };
51 
52 /**
53  * \brief Laplacian of function.
54  *
55  * This is Laplacian of \f$u\f$, it is calculated using formula
56  * \f[
57  * \nabla^2 u(x,y,z) = \nabla \cdot \nabla u
58  * \frac{\partial^2 u}{\partial x^2}+
59  * \frac{\partial^2 u}{\partial y^2}+
60  * \frac{\partial^2 u}{\partial z^2}
61  * \f]
62  *
63  */
65  double operator()(const double x, const double y, const double z) const {
66  return 4 + 6 * z;
67  }
68 };
69 
70 int main(int argc, char *argv[]) {
71 
72  // Initialize PETSc
73  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
74 
75  try {
76 
77  // Create MoAB database
78  moab::Core moab_core; // create database
79  moab::Interface &moab = moab_core; // create interface to database
80 
81  // Get command line options
82  int order = 3; // default approximation order
83  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
84  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
85  "none");
86  // Set approximation order
87  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
88  PETSC_NULL);
89  // Set testing (used by CTest)
90  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
91  &flg_test, PETSC_NULL);
92  ierr = PetscOptionsEnd();
93  CHKERRG(ierr);
94 
95  // Create MoFEM database and link it to MoAB
96  MoFEM::Core mofem_core(moab); // create database
97  MoFEM::Interface &m_field = mofem_core; // create interface to database
98  // Register DM Manager
99  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
100 
101  // Create vector to store approximation global error
102  Vec global_error;
103  CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
104 
105  // First we crate elements, implementation of elements is problem independent,
106  // we do not know yet what fields are present in the problem, or
107  // even we do not decided yet what approximation base or spaces we
108  // are going to use. Implementation of element is free from
109  // those constrains and can be used in different context.
110 
111  // Elements used by KSP & DM to assemble system of equations
112  boost::shared_ptr<ForcesAndSourcesCore> domain_lhs_fe; ///< Volume element for the matrix
113  boost::shared_ptr<ForcesAndSourcesCore> boundary_lhs_fe; ///< Boundary element for the matrix
114  boost::shared_ptr<ForcesAndSourcesCore> domain_rhs_fe; ///< Volume element to assemble vector
115  boost::shared_ptr<ForcesAndSourcesCore> boundary_rhs_fe; ///< Volume element to assemble vector
116  boost::shared_ptr<ForcesAndSourcesCore> domain_error; ///< Volume element evaluate error
117  boost::shared_ptr<PoissonExample::PostProcFE>
118  post_proc_volume; ///< Volume element to Post-process results
119  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
120  {
121  // Add problem specific operators the generic finite elements to calculate matrices and vectors.
124  ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
125  boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe);
126  // Add problem specific operators the generic finite elements to calculate error on elements and global error
127  // in H1 norm
130  global_error, domain_error);
131  // Post-process results
133  .creatFEToPostProcessResults(post_proc_volume);
134  }
135 
136  // Get simple interface is simplified version enabling quick and
137  // easy construction of problem.
138  Simple *simple_interface;
139  // Query interface and get pointer to Simple interface
140  CHKERR m_field.getInterface(simple_interface);
141 
142  // Build problem with simple interface
143  {
144 
145  // Get options for simple interface from command line
146  CHKERR simple_interface->getOptions();
147  // Load mesh file to database
148  CHKERR simple_interface->loadFile();
149 
150  // Add field on domain and boundary. Field is declared by space and base and rank. space
151  // can be H1. Hcurl, Hdiv and L2, base can be AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more,
152  // where rank is number of coefficients for dof.
153  //
154  // Simple interface assumes that there are four types of field; 1) defined
155  // on body domain, 2) fields defined on body boundary, 3) skeleton field defined
156  // on finite element skeleton. Finally data field is defined on body domain. Data field
157  // is not solved but used for post-process or to keep material parameters, etc. Here
158  // we using it to calculate approximation error on elements.
159 
160  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is used
161  // to construct base functions.
162  CHKERR simple_interface->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,1);
163  // Add Lagrange multiplier field on body boundary
164  CHKERR simple_interface->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,1);
165  // Add error (data) field, we need only L2 norm. Later order is set to 0, so this
166  // is piecewise discontinuous constant approx., i.e. 1 DOF for element. You can use
167  // more DOFs and collate moments of error to drive anisotropic h/p-adaptivity, however
168  // this is beyond this example.
169  CHKERR simple_interface->addDataField("ERROR",L2,AINSWORTH_LEGENDRE_BASE,1);
170 
171  // Set fields order domain and boundary fields.
172  CHKERR simple_interface->setFieldOrder("U",order); // to approximate function
173  CHKERR simple_interface->setFieldOrder("L",order); // to Lagrange multipliers
174  CHKERR simple_interface->setFieldOrder("ERROR",0); // approximation order for error
175 
176  // Setup problem. At that point database is constructed, i.e. fields, finite elements and
177  // problem data structures with local and global indexing.
178  CHKERR simple_interface->setUp();
179 
180  }
181 
182  // Get access to PETSC-MoFEM DM manager.
183  // or more derails see <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
184  // Form that point internal MoFEM data structures are linked with PETSc interface. If
185  // DM functions contains string MoFEM is is MoFEM specific DM interface function,
186  // otherwise DM function part of standard PETSc interface.
187 
188  DM dm;
189  // Get dm
190  CHKERR simple_interface->getDM(&dm);
191 
192  // Set KSP context for DM. At that point only elements are added to DM operators.
193  // Calculations of matrices and vectors is executed by KSP solver. This part
194  // of the code makes connection between implementation of finite elements and
195  // data operators with finite element declarations in DM manager. From that
196  // point DM takes responsibility for executing elements, calculations of
197  // matrices and vectors, and solution of the problem.
198  {
199  // Set operators for KSP solver
201  dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
203  dm, simple_interface->getBoundaryFEName(), boundary_lhs_fe, null,
204  null);
205  // Set calculation of the right hand side vector for KSP solver
206  CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getDomainFEName(),
207  domain_rhs_fe, null, null);
208  CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getBoundaryFEName(),
209  boundary_rhs_fe, null, null);
210  }
211 
212  // Solve problem, only PETEc interface here.
213  {
214 
215  // Create the right hand side vector and vector of unknowns
216  Vec F,D;
217  CHKERR DMCreateGlobalVector(dm,&F);
218  // Create unknown vector by creating duplicate copy of F vector. only
219  // structure is duplicated no values.
220  CHKERR VecDuplicate(F,&D);
221 
222  // Create solver and link it to DM
223  KSP solver;
224  CHKERR KSPCreate(PETSC_COMM_WORLD,&solver);
225  CHKERR KSPSetFromOptions(solver);
226  CHKERR KSPSetDM(solver,dm);
227  // Set-up solver, is type of solver and pre-conditioners
228  CHKERR KSPSetUp(solver);
229  // At solution process, KSP solver using DM creates matrices, Calculate
230  // values of the left hand side and the right hand side vector. then
231  // solves system of equations. Results are stored in vector D.
232  CHKERR KSPSolve(solver,F,D);
233 
234  // Scatter solution on the mesh. Stores unknown vector on field on the mesh.
235  CHKERR DMoFEMMeshToGlobalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE);
236 
237  // Clean data. Solver and vector are not needed any more.
238  CHKERR KSPDestroy(&solver);
239  CHKERR VecDestroy(&D);
240  CHKERR VecDestroy(&F);
241  }
242 
243  // Calculate error
244  {
245  // Loop over all elements in mesh, and run users operators on each element.
246  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
247  domain_error);
249  global_error);
250  CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
251  if (flg_test == PETSC_TRUE) {
252  CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
253  }
254  }
255 
256  {
257  // Loop over all elements in the mesh and for each execute post_proc_volume
258  // element and operators on it.
259  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
260  post_proc_volume);
261  // Write results
262  post_proc_volume->writeFile("out_vol.h5m");
263  }
264 
265  // Destroy DM, no longer needed.
266  CHKERR DMDestroy(&dm);
267 
268  // Destroy ghost vector
269  CHKERR VecDestroy(&global_error);
270  }
271  CATCH_ERRORS;
272 
273  // finish work cleaning memory, getting statistics, etc.
275 
276  return 0;
277 
278 }
PoissonExample::AuxFunctions
Definition: AuxPoissonFunctions.hpp:14
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMKSPSetComputeRHS
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: DMMoFEM.cpp:637
FTensor::Tensor1< double, 3 >
AuxPoissonFunctions.hpp
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
PoissonExample::CreateFiniteElements
Create finite elements instances.
Definition: PoissonOperators.hpp:858
MoFEM::DMMoFEMKSPSetComputeOperators
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: DMMoFEM.cpp:678
ExactFunctionGrad
Exact gradient.
Definition: analytical_nonlinear_poisson.cpp:43
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
PoissonExample::AuxFunctions::createGhostVec
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
Definition: AuxPoissonFunctions.hpp:30
ExactLaplacianFunction::operator()
double operator()(const double x, const double y, const double z) const
Definition: analytical_poisson.cpp:65
ExactFunction::operator()
double operator()(const double x, const double y, const double z) const
Definition: analytical_poisson.cpp:33
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:693
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
PoissonExample::CreateFiniteElements::creatFEToPostProcessResults
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
Definition: PoissonOperators.hpp:944
ExactLaplacianFunction
Laplacian of function.
Definition: analytical_nonlinear_poisson.cpp:58
PoissonExample::AuxFunctions::testError
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
Definition: AuxPoissonFunctions.hpp:73
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::Simple::addDomainField
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)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
PoissonExample::AuxFunctions::printError
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Definition: AuxPoissonFunctions.hpp:60
main
int main(int argc, char *argv[])
Definition: analytical_poisson.cpp:70
MoFEM::Simple::addDataField
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)
Add field on domain.
Definition: Simple.cpp:338
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:491
ExactFunction
Function.
Definition: analytical_nonlinear_poisson.cpp:34
MoFEM::CoreTmp< 0 >::Initialize
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
help
static char help[]
Definition: analytical_poisson.cpp:18
PoissonOperators.hpp
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PoissonExample::CreateFiniteElements::createFEToEvaluateError
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.
Definition: PoissonOperators.hpp:908
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::Simple::getBoundaryFEName
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:375
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::Simple::addBoundaryField
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)
Add field on boundary.
Definition: Simple.cpp:300
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
PoissonExample::CreateFiniteElements::createFEToAssembleMatrixAndVector
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.
Definition: PoissonOperators.hpp:865
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:629
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
PoissonExample::AuxFunctions::assembleGhostVector
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
Definition: AuxPoissonFunctions.hpp:43
F
@ F
Definition: free_surface.cpp:394
ExactFunctionGrad::operator()
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double z) const
Definition: analytical_poisson.cpp:42