v0.15.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
14
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 */
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 */
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 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
100 "none");
101 // Set approximation order
102 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
103 PETSC_NULLPTR);
104 // Set testing (used by CTest)
105 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
106 &flg_test, PETSC_NULLPTR);
107 PetscOptionsEnd();
108
109 // Create MoFEM database and link it to MoAB
110 MoFEM::Core mofem_core(moab); // create database
111 MoFEM::Interface &m_field = mofem_core; // create interface to database
112 // Register DM Manager
113 CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
114
115 // Create vector to store approximation global error
116 Vec global_error;
117 CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
118
119 // First we crate elements, implementation of elements is problem
120 // independent, we do not know yet what fields are present in the problem,
121 // or even we do not decided yet what approximation base or spaces we are
122 // going to use. Implementation of element is free from those constrains and
123 // can be used in different context.
124
125 // Elements used by KSP & DM to assemble system of equations
126 boost::shared_ptr<ForcesAndSourcesCore>
127 domain_lhs_fe; ///< Volume element for the matrix
128 boost::shared_ptr<ForcesAndSourcesCore>
129 boundary_lhs_fe; ///< Boundary element for the matrix
130 boost::shared_ptr<ForcesAndSourcesCore>
131 domain_rhs_fe; ///< Volume element to assemble vector
132 boost::shared_ptr<ForcesAndSourcesCore>
133 boundary_rhs_fe; ///< Volume element to assemble vector
134 boost::shared_ptr<ForcesAndSourcesCore>
135 domain_error; ///< Volume element evaluate error
136 boost::shared_ptr<PoissonExample::PostProcFE>
137 post_proc_volume; ///< Volume element to Post-process results
138 boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
139 {
140 // Add problem specific operators the generic finite elements to calculate
141 // matrices and vectors.
145 domain_lhs_fe, boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe,
147 // Add problem specific operators the generic finite elements to calculate
148 // error on elements and global error in H1 norm
151 global_error, domain_error);
152 // Post-process results
154 .creatFEToPostProcessResults(post_proc_volume);
155 }
156
157 // Get simple interface is simplified version enabling quick and
158 // easy construction of problem.
159 Simple *simple_interface;
160 // Query interface and get pointer to Simple interface
161 CHKERR m_field.getInterface(simple_interface);
162
163 // Build problem with simple interface
164 {
165
166 // Get options for simple interface from command line
167 CHKERR simple_interface->getOptions();
168 // Load mesh file to database
169 CHKERR simple_interface->loadFile();
170
171 // Add field on domain and boundary. Field is declared by space and base
172 // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
173 // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
174 // number of coefficients for dof.
175 //
176 // Simple interface assumes that there are four types of field; 1) defined
177 // on body domain, 2) fields defined on body boundary, 3) skeleton field
178 // defined on finite element skeleton. Finally data field is defined on
179 // body domain. Data field is not solved but used for post-process or to
180 // keep material parameters, etc. Here we using it to calculate
181 // approximation error on elements.
182
183 // Add domain field "U" in space H1 and Legendre base, Ainsworth recipe is
184 // used to construct base functions.
185 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
186 1);
187 // Add Lagrange multiplier field on body boundary
188 CHKERR simple_interface->addBoundaryField("L", H1,
190 // Add error (data) field, we need only L2 norm. Later order is set to 0,
191 // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
192 // element. You can use more DOFs and collate moments of error to drive
193 // anisotropic h/p-adaptivity, however this is beyond this example.
194 CHKERR simple_interface->addDataField("ERROR", L2,
196
197 // Set fields order domain and boundary fields.
198 CHKERR simple_interface->setFieldOrder("U",
199 order); // to approximate function
200 CHKERR simple_interface->setFieldOrder("L",
201 order); // to Lagrange multipliers
202 CHKERR simple_interface->setFieldOrder(
203 "ERROR", 0); // approximation order for error
204
205 // Setup problem. At that point database is constructed, i.e. fields,
206 // finite elements and problem data structures with local and global
207 // indexing.
208 CHKERR simple_interface->setUp();
209 }
210
211 // Get access to PETSC-MoFEM DM manager.
212 // or more derails see
213 // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
214 // Form that point internal MoFEM data structures are linked with PETSc
215 // interface. If DM functions contains string MoFEM is is MoFEM specific DM
216 // interface function, otherwise DM function part of standard PETSc
217 // interface.
218
219 DM dm;
220 // Get dm
221 CHKERR simple_interface->getDM(&dm);
222
223 // Set KSP context for DM. At that point only elements are added to DM
224 // operators. Calculations of matrices and vectors is executed by KSP
225 // solver. This part of the code makes connection between implementation of
226 // finite elements and data operators with finite element declarations in DM
227 // manager. From that point DM takes responsibility for executing elements,
228 // calculations of matrices and vectors, and solution of the problem.
229 {
230 // Set operators for SNES solver
231 CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getDomainFEName(),
232 domain_lhs_fe, null, null);
233 CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getBoundaryFEName(),
234 boundary_lhs_fe, null, null);
235 // Set calculation of the right hand side vector for KSP solver
236 CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getDomainFEName(),
237 domain_rhs_fe, null, null);
238 CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getBoundaryFEName(),
239 boundary_rhs_fe, null, null);
240 }
241
242 // Solve problem, only PETEc interface here.
243 {
244
245 // Create the right hand side vector and vector of unknowns
246 Vec F, D;
247 CHKERR DMCreateGlobalVector(dm, &F);
248 // Create unknown vector by creating duplicate copy of F vector. only
249 // structure is duplicated no values.
250 CHKERR VecDuplicate(F, &D);
251
252 // Create solver and link it to DM
253 SNES solver;
254 CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
255 CHKERR SNESSetFromOptions(solver);
256 CHKERR SNESSetDM(solver, dm);
257 // Set-up solver, is type of solver and pre-conditioners
258 CHKERR SNESSetUp(solver);
259 // At solution process, KSP solver using DM creates matrices, Calculate
260 // values of the left hand side and the right hand side vector. then
261 // solves system of equations. Results are stored in vector D.
262 CHKERR SNESSolve(solver, F, D);
263
264 // Scatter solution on the mesh. Stores unknown vector on field on the
265 // mesh.
266 CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
267
268 // Clean data. Solver and vector are not needed any more.
269 CHKERR SNESDestroy(&solver);
270 CHKERR VecDestroy(&D);
271 CHKERR VecDestroy(&F);
272 }
273
274 // Calculate error
275 {
276 // Loop over all elements in mesh, and run users operators on each
277 // element.
278 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
279 domain_error);
281 global_error);
282 CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
283 if (flg_test == PETSC_TRUE) {
284 CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
285 }
286 }
287
288 {
289 // Loop over all elements in the mesh and for each execute
290 // post_proc_volume element and operators on it.
291 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
292 post_proc_volume);
293 // Write results
294 post_proc_volume->writeFile("out_vol.h5m");
295 }
296
297 // Destroy DM, no longer needed.
298 CHKERR DMDestroy(&dm);
299
300 // Destroy ghost vector
301 CHKERR VecDestroy(&global_error);
302 }
304
305 // finish work cleaning memory, getting statistics, etc.
307
308 return 0;
309}
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 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:708
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:749
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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
double D
double operator()(const double u)
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: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 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