v0.15.0
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
14#include <MoFEM.hpp>
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 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
85 "none");
86 // Set approximation order
87 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
88 PETSC_NULLPTR);
89 // Set testing (used by CTest)
90 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
91 &flg_test, PETSC_NULLPTR);
92 PetscOptionsEnd();
93
94 // Create MoFEM database and link it to MoAB
95 MoFEM::Core mofem_core(moab); // create database
96 MoFEM::Interface &m_field = mofem_core; // create interface to database
97 // Register DM Manager
98 CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
99
100 // Create vector to store approximation global error
101 Vec global_error;
102 CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
103
104 // First we crate elements, implementation of elements is problem independent,
105 // we do not know yet what fields are present in the problem, or
106 // even we do not decided yet what approximation base or spaces we
107 // are going to use. Implementation of element is free from
108 // those constrains and can be used in different context.
109
110 // Elements used by KSP & DM to assemble system of equations
111 boost::shared_ptr<ForcesAndSourcesCore> domain_lhs_fe; ///< Volume element for the matrix
112 boost::shared_ptr<ForcesAndSourcesCore> boundary_lhs_fe; ///< Boundary element for the matrix
113 boost::shared_ptr<ForcesAndSourcesCore> domain_rhs_fe; ///< Volume element to assemble vector
114 boost::shared_ptr<ForcesAndSourcesCore> boundary_rhs_fe; ///< Volume element to assemble vector
115 boost::shared_ptr<ForcesAndSourcesCore> domain_error; ///< Volume element evaluate error
116 boost::shared_ptr<PoissonExample::PostProcFE>
117 post_proc_volume; ///< Volume element to Post-process results
118 boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
119 {
120 // Add problem specific operators the generic finite elements to calculate matrices and vectors.
123 ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
124 boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe);
125 // Add problem specific operators the generic finite elements to calculate error on elements and global error
126 // in H1 norm
129 global_error, domain_error);
130 // Post-process results
132 .creatFEToPostProcessResults(post_proc_volume);
133 }
134
135 // Get simple interface is simplified version enabling quick and
136 // easy construction of problem.
137 Simple *simple_interface;
138 // Query interface and get pointer to Simple interface
139 CHKERR m_field.getInterface(simple_interface);
140
141 // Build problem with simple interface
142 {
143
144 // Get options for simple interface from command line
145 CHKERR simple_interface->getOptions();
146 // Load mesh file to database
147 CHKERR simple_interface->loadFile();
148
149 // Add field on domain and boundary. Field is declared by space and base and rank. space
150 // can be H1. Hcurl, Hdiv and L2, base can be AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more,
151 // where rank is number of coefficients for dof.
152 //
153 // Simple interface assumes that there are four types of field; 1) defined
154 // on body domain, 2) fields defined on body boundary, 3) skeleton field defined
155 // on finite element skeleton. Finally data field is defined on body domain. Data field
156 // is not solved but used for post-process or to keep material parameters, etc. Here
157 // we using it to calculate approximation error on elements.
158
159 // Add domain field "U" in space H1 and Legendre base, Ainsworth recipe is used
160 // to construct base functions.
161 CHKERR simple_interface->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,1);
162 // Add Lagrange multiplier field on body boundary
163 CHKERR simple_interface->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,1);
164 // Add error (data) field, we need only L2 norm. Later order is set to 0, so this
165 // is piecewise discontinuous constant approx., i.e. 1 DOF for element. You can use
166 // more DOFs and collate moments of error to drive anisotropic h/p-adaptivity, however
167 // this is beyond this example.
168 CHKERR simple_interface->addDataField("ERROR",L2,AINSWORTH_LEGENDRE_BASE,1);
169
170 // Set fields order domain and boundary fields.
171 CHKERR simple_interface->setFieldOrder("U",order); // to approximate function
172 CHKERR simple_interface->setFieldOrder("L",order); // to Lagrange multipliers
173 CHKERR simple_interface->setFieldOrder("ERROR",0); // approximation order for error
174
175 // Setup problem. At that point database is constructed, i.e. fields, finite elements and
176 // problem data structures with local and global indexing.
177 CHKERR simple_interface->setUp();
178
179 }
180
181 // Get access to PETSC-MoFEM DM manager.
182 // or more derails see <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
183 // Form that point internal MoFEM data structures are linked with PETSc interface. If
184 // DM functions contains string MoFEM is is MoFEM specific DM interface function,
185 // otherwise DM function part of standard PETSc interface.
186
187 DM dm;
188 // Get dm
189 CHKERR simple_interface->getDM(&dm);
190
191 // Set KSP context for DM. At that point only elements are added to DM operators.
192 // Calculations of matrices and vectors is executed by KSP solver. This part
193 // of the code makes connection between implementation of finite elements and
194 // data operators with finite element declarations in DM manager. From that
195 // point DM takes responsibility for executing elements, calculations of
196 // matrices and vectors, and solution of the problem.
197 {
198 // Set operators for KSP solver
200 dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
202 dm, simple_interface->getBoundaryFEName(), boundary_lhs_fe, null,
203 null);
204 // Set calculation of the right hand side vector for KSP solver
205 CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getDomainFEName(),
206 domain_rhs_fe, null, null);
207 CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getBoundaryFEName(),
208 boundary_rhs_fe, null, null);
209 }
210
211 // Solve problem, only PETEc interface here.
212 {
213
214 // Create the right hand side vector and vector of unknowns
215 Vec F,D;
216 CHKERR DMCreateGlobalVector(dm,&F);
217 // Create unknown vector by creating duplicate copy of F vector. only
218 // structure is duplicated no values.
219 CHKERR VecDuplicate(F,&D);
220
221 // Create solver and link it to DM
222 KSP solver;
223 CHKERR KSPCreate(PETSC_COMM_WORLD,&solver);
224 CHKERR KSPSetFromOptions(solver);
225 CHKERR KSPSetDM(solver,dm);
226 // Set-up solver, is type of solver and pre-conditioners
227 CHKERR KSPSetUp(solver);
228 // At solution process, KSP solver using DM creates matrices, Calculate
229 // values of the left hand side and the right hand side vector. then
230 // solves system of equations. Results are stored in vector D.
231 CHKERR KSPSolve(solver,F,D);
232
233 // Scatter solution on the mesh. Stores unknown vector on field on the mesh.
234 CHKERR DMoFEMMeshToGlobalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE);
235
236 // Clean data. Solver and vector are not needed any more.
237 CHKERR KSPDestroy(&solver);
238 CHKERR VecDestroy(&D);
239 CHKERR VecDestroy(&F);
240 }
241
242 // Calculate error
243 {
244 // Loop over all elements in mesh, and run users operators on each element.
245 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
246 domain_error);
248 global_error);
249 CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
250 if (flg_test == PETSC_TRUE) {
251 CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
252 }
253 }
254
255 {
256 // Loop over all elements in the mesh and for each execute post_proc_volume
257 // element and operators on it.
258 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
259 post_proc_volume);
260 // Write results
261 post_proc_volume->writeFile("out_vol.h5m");
262 }
263
264 // Destroy DM, no longer needed.
265 CHKERR DMDestroy(&dm);
266
267 // Destroy ghost vector
268 CHKERR VecDestroy(&global_error);
269 }
271
272 // finish work cleaning memory, getting statistics, etc.
274
275 return 0;
276
277}
int main()
static char help[]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
#define CHKERR
Inline error check.
constexpr int order
@ F
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set compute operator for KSP solver via sub-matrix and IS.
Definition DMMoFEM.cpp:627
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:576
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:525
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:668
double D
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:118
Deprecated interface functions.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
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:261
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:396
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
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:355
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
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:393
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:736
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:389
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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.