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