v0.14.0
Loading...
Searching...
No Matches
analytical_nonlinear_poisson.cpp
Go to the documentation of this file.
1/**
2 * \file analytical_nonlinear_poisson.cpp
3 * \ingroup mofem_simple_interface
4 * \example analytical_nonlinear_poisson.cpp
5 *
6 * For more information and detailed explain of this
7 * example see \ref poisson_tut4
8 *
9 *
10 */
11
12
13
15
16#include <PoissonOperators.hpp>
17
19
20static char help[] = "...\n\n";
21
22/**
23 * \brief Function
24 *
25 * This is prescribed exact function. If this function is given by polynomial
26 * order of "p" and order of approximation is "p" or higher, solution of
27 * finite element method is exact (with machine precision).
28 *
29 * \f[
30 * u = 1+x+2y+3z
31 * \f]
32 *
33 */
34struct ExactFunction {
35 double operator()(const double x, const double y, const double z) const {
36 return 1 + x + y + pow(z, 3);
37 }
38};
39
40/**
41 * \brief Exact gradient
42 */
43struct ExactFunctionGrad {
44 FTensor::Tensor1<double, 3> operator()(const double x, const double y,
45 const double z) const {
47 grad(0) = 1;
48 grad(1) = 1;
49 grad(2) = 3 * z * z;
50 return grad;
51 }
52};
53
54/**
55 * \brief Laplacian of function.
56 *
57 */
59 double operator()(const double x, const double y, const double z) const {
60 return 0.4e1 + (double)(4 * x) + (double)(4 * y) + 0.4e1 * pow(z, 0.3e1) +
61 0.3e1 *
62 (0.6e1 * z * z + 0.6e1 * (double)x * z * z +
63 0.6e1 * (double)y * z * z + 0.6e1 * pow(z, 0.5e1)) *
64 z * z +
65 0.6e1 *
66 (0.2e1 + (double)(2 * x) + (double)(2 * y) +
67 0.2e1 * pow(z, 0.3e1) + (double)(x * x) + (double)(2 * x * y) +
68 0.2e1 * (double)x * pow(z, 0.3e1) + (double)(y * y) +
69 0.2e1 * (double)y * pow(z, 0.3e1) + pow(z, 0.6e1)) *
70 z;
71 }
72};
73
74struct FunA {
75 double operator()(const double u) { return 1 + u * u; }
76};
77
78struct DiffFunA {
79 double operator()(const double u) { return 2 * u; }
80};
81
83 int operator()(int, int, int p) const { return 2 * (p + 1); }
84};
85
86int main(int argc, char *argv[]) {
87
88 // Initialize PETSc
89 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
90 // Create MoAB database
91 moab::Core moab_core; // create database
92 moab::Interface &moab = moab_core; // create interface to database
93
94 try {
95
96 // Get command line options
97 int order = 3; // default approximation order
98 PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
99 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
100 "none");
101 // Set approximation order
102 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
103 PETSC_NULL);
104 // Set testing (used by CTest)
105 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
106 &flg_test, PETSC_NULL);
107 ierr = PetscOptionsEnd();
109
110 // Create MoFEM database and link it to MoAB
111 MoFEM::Core mofem_core(moab); // create database
112 MoFEM::Interface &m_field = mofem_core; // create interface to database
113 // Register DM Manager
114 CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
115
116 // Create vector to store approximation global error
117 Vec global_error;
118 CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
119
120 // First we crate elements, implementation of elements is problem
121 // independent, we do not know yet what fields are present in the problem,
122 // or even we do not decided yet what approximation base or spaces we are
123 // going to use. Implementation of element is free from those constrains and
124 // can be used in different context.
125
126 // Elements used by KSP & DM to assemble system of equations
127 boost::shared_ptr<ForcesAndSourcesCore>
128 domain_lhs_fe; ///< Volume element for the matrix
129 boost::shared_ptr<ForcesAndSourcesCore>
130 boundary_lhs_fe; ///< Boundary element for the matrix
131 boost::shared_ptr<ForcesAndSourcesCore>
132 domain_rhs_fe; ///< Volume element to assemble vector
133 boost::shared_ptr<ForcesAndSourcesCore>
134 boundary_rhs_fe; ///< Volume element to assemble vector
135 boost::shared_ptr<ForcesAndSourcesCore>
136 domain_error; ///< Volume element evaluate error
137 boost::shared_ptr<PoissonExample::PostProcFE>
138 post_proc_volume; ///< Volume element to Post-process results
139 boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
140 {
141 // Add problem specific operators the generic finite elements to calculate
142 // matrices and vectors.
146 domain_lhs_fe, boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe,
148 // Add problem specific operators the generic finite elements to calculate
149 // error on elements and global error in H1 norm
152 global_error, domain_error);
153 // Post-process results
155 .creatFEToPostProcessResults(post_proc_volume);
156 }
157
158 // Get simple interface is simplified version enabling quick and
159 // easy construction of problem.
160 Simple *simple_interface;
161 // Query interface and get pointer to Simple interface
162 CHKERR m_field.getInterface(simple_interface);
163
164 // Build problem with simple interface
165 {
166
167 // Get options for simple interface from command line
168 CHKERR simple_interface->getOptions();
169 // Load mesh file to database
170 CHKERR simple_interface->loadFile();
171
172 // Add field on domain and boundary. Field is declared by space and base
173 // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
174 // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
175 // number of coefficients for dof.
176 //
177 // Simple interface assumes that there are four types of field; 1) defined
178 // on body domain, 2) fields defined on body boundary, 3) skeleton field
179 // defined on finite element skeleton. Finally data field is defined on
180 // body domain. Data field is not solved but used for post-process or to
181 // keep material parameters, etc. Here we using it to calculate
182 // approximation error on elements.
183
184 // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
185 // used to construct base functions.
186 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
187 1);
188 // Add Lagrange multiplier field on body boundary
189 CHKERR simple_interface->addBoundaryField("L", H1,
191 // Add error (data) field, we need only L2 norm. Later order is set to 0,
192 // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
193 // element. You can use more DOFs and collate moments of error to drive
194 // anisotropic h/p-adaptivity, however this is beyond this example.
195 CHKERR simple_interface->addDataField("ERROR", L2,
197
198 // Set fields order domain and boundary fields.
199 CHKERR simple_interface->setFieldOrder("U",
200 order); // to approximate function
201 CHKERR simple_interface->setFieldOrder("L",
202 order); // to Lagrange multipliers
203 CHKERR simple_interface->setFieldOrder(
204 "ERROR", 0); // approximation order for error
205
206 // Setup problem. At that point database is constructed, i.e. fields,
207 // finite elements and problem data structures with local and global
208 // indexing.
209 CHKERR simple_interface->setUp();
210 }
211
212 // Get access to PETSC-MoFEM DM manager.
213 // or more derails see
214 // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
215 // Form that point internal MoFEM data structures are linked with PETSc
216 // interface. If DM functions contains string MoFEM is is MoFEM specific DM
217 // interface function, otherwise DM function part of standard PETSc
218 // interface.
219
220 DM dm;
221 // Get dm
222 CHKERR simple_interface->getDM(&dm);
223
224 // Set KSP context for DM. At that point only elements are added to DM
225 // operators. Calculations of matrices and vectors is executed by KSP
226 // solver. This part of the code makes connection between implementation of
227 // finite elements and data operators with finite element declarations in DM
228 // manager. From that point DM takes responsibility for executing elements,
229 // calculations of matrices and vectors, and solution of the problem.
230 {
231 // Set operators for SNES solver
232 CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getDomainFEName(),
233 domain_lhs_fe, null, null);
234 CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getBoundaryFEName(),
235 boundary_lhs_fe, null, null);
236 // Set calculation of the right hand side vector for KSP solver
237 CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getDomainFEName(),
238 domain_rhs_fe, null, null);
239 CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getBoundaryFEName(),
240 boundary_rhs_fe, null, null);
241 }
242
243 // Solve problem, only PETEc interface here.
244 {
245
246 // Create the right hand side vector and vector of unknowns
247 Vec F, D;
248 CHKERR DMCreateGlobalVector(dm, &F);
249 // Create unknown vector by creating duplicate copy of F vector. only
250 // structure is duplicated no values.
251 CHKERR VecDuplicate(F, &D);
252
253 // Create solver and link it to DM
254 SNES solver;
255 CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
256 CHKERR SNESSetFromOptions(solver);
257 CHKERR SNESSetDM(solver, dm);
258 // Set-up solver, is type of solver and pre-conditioners
259 CHKERR SNESSetUp(solver);
260 // At solution process, KSP solver using DM creates matrices, Calculate
261 // values of the left hand side and the right hand side vector. then
262 // solves system of equations. Results are stored in vector D.
263 CHKERR SNESSolve(solver, F, D);
264
265 // Scatter solution on the mesh. Stores unknown vector on field on the
266 // mesh.
267 CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
268
269 // Clean data. Solver and vector are not needed any more.
270 CHKERR SNESDestroy(&solver);
271 CHKERR VecDestroy(&D);
272 CHKERR VecDestroy(&F);
273 }
274
275 // Calculate error
276 {
277 // Loop over all elements in mesh, and run users operators on each
278 // element.
279 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
280 domain_error);
282 global_error);
283 CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
284 if (flg_test == PETSC_TRUE) {
285 CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
286 }
287 }
288
289 {
290 // Loop over all elements in the mesh and for each execute
291 // post_proc_volume element and operators on it.
292 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
293 post_proc_volume);
294 // Write results
295 post_proc_volume->writeFile("out_vol.h5m");
296 }
297
298 // Destroy DM, no longer needed.
299 CHKERR DMDestroy(&dm);
300
301 // Destroy ghost vector
302 CHKERR VecDestroy(&global_error);
303 }
305
306 // finish work cleaning memory, getting statistics, etc.
308
309 return 0;
310}
static Index< 'p', 3 > p
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
constexpr int order
@ F
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:704
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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:572
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
double D
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
double operator()(const double u)
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
double operator()(const double u)
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:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
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:341
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 createFEToAssembleMatrixAndVectorForNonlinearProblem(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::function< double(const double)> a, boost::function< double(const double)> diff_a, 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, ForcesAndSourcesCore::RuleHookFun vol_rule, ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule(), bool trans=true) const
Create finite element to calculate matrix and vectors.
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.
int operator()(int, int, int p) const