v0.14.0
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 #include <BasicFiniteElements.hpp>
15
16 #include <PoissonOperators.hpp>
17
18 #include <AuxPoissonFunctions.hpp>
19
20 static 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  */
34 struct 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  */
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;
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
74 struct FunA {
75  double operator()(const double u) { return 1 + u * u; }
76 };
77
78 struct 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
86 int 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();
108  CHKERRG(ierr)
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,
147  VolRuleNonlinear());
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
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  }
304  CATCH_ERRORS;
305
306  // finish work cleaning memory, getting statistics, etc.
308
309  return 0;
310 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
main
int main(int argc, char *argv[])
Definition: analytical_nonlinear_poisson.cpp:86
VolRuleNonlinear::operator()
int operator()(int, int, int p) const
Definition: analytical_nonlinear_poisson.cpp:83
DiffFunA::operator()
double operator()(const double u)
Definition: analytical_nonlinear_poisson.cpp:79
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FTensor::Tensor1< double, 3 >
AuxPoissonFunctions.hpp
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
Definition: Simple.cpp:194
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
PoissonExample::CreateFiniteElements
Create finite elements instances.
Definition: PoissonOperators.hpp:858
Definition: analytical_nonlinear_poisson.cpp:43
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
PoissonExample::AuxFunctions::createGhostVec
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
Definition: AuxPoissonFunctions.hpp:30
ExactLaplacianFunction::operator()
double operator()(const double x, const double y, const double z) const
Definition: analytical_nonlinear_poisson.cpp:59
ExactFunction::operator()
double operator()(const double x, const double y, const double z) const
Definition: analytical_nonlinear_poisson.cpp:35
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:693
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
PoissonExample::CreateFiniteElements::creatFEToPostProcessResults
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
Definition: PoissonOperators.hpp:944
ExactLaplacianFunction
Laplacian of function.
Definition: analytical_nonlinear_poisson.cpp:58
FunA
Definition: analytical_nonlinear_poisson.cpp:74
double
PoissonExample::AuxFunctions::testError
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
Definition: AuxPoissonFunctions.hpp:73
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
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
DiffFunA
Definition: analytical_nonlinear_poisson.cpp:78
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
MoFEM::DMMoFEMSNESSetJacobian
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:759
PoissonExample::AuxFunctions::printError
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Definition: AuxPoissonFunctions.hpp:60
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:338
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:491
FunA::operator()
double operator()(const double u)
Definition: analytical_nonlinear_poisson.cpp:75
ExactFunction
Function.
Definition: analytical_nonlinear_poisson.cpp:34
MoFEM::CoreTmp< 0 >::Initialize
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
PoissonOperators.hpp
help
static char help[]
Definition: analytical_nonlinear_poisson.cpp:20
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PoissonExample::CreateFiniteElements::createFEToAssembleMatrixAndVectorForNonlinearProblem
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.
Definition: PoissonOperators.hpp:1010
PoissonExample::CreateFiniteElements::createFEToEvaluateError
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.
Definition: PoissonOperators.hpp:908
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::Simple::getBoundaryFEName
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:375
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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:300
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:629
MoFEM::DMoFEMLoopFiniteElements
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:586
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
PoissonExample::AuxFunctions::assembleGhostVector
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
Definition: AuxPoissonFunctions.hpp:43
F
@ F
Definition: free_surface.cpp:394