v0.13.2
Loading...
Searching...
No Matches
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
15#include <PoissonOperators.hpp>
17
18static 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 */
32struct 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 */
41struct 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
70int 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();
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 }
272
273 // finish work cleaning memory, getting statistics, etc.
275
276 return 0;
277
278}
int main()
Definition: adol-c_atom.cpp:46
static char help[]
#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: DMMoFEM.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: DMMoFEM.cpp:605
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:554
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:503
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:646
double D
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
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: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.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:327
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:320
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:668
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
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:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:320
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.