v0.9.1
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 /* This file is part of MoFEM.
13 * MoFEM is free software: you can redistribute it and/or modify it under
14 * the terms of the GNU Lesser General Public License as published by the
15 * Free Software Foundation, either version 3 of the License, or (at your
16 * option) any later version.
17 *
18 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21 * License for more details.
22 *
23 * You should have received a copy of the GNU Lesser General Public
24 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #include <BasicFiniteElements.hpp>
27 #include <PoissonOperators.hpp>
28 #include <AuxPoissonFunctions.hpp>
29 
30 static char help[] = "...\n\n";
31 
32 /**
33  * \brief Function
34  *
35  * This is prescribed exact function. If this function is given by polynomial
36  * order of "p" and order of approximation is "p" or higher, solution of
37  * finite element method is exact (with machine precision).
38  *
39  * \f[
40  * u = 1+x^2+y^2+z^3
41  * \f]
42  *
43  */
44 struct ExactFunction {
45  double operator()(const double x, const double y, const double z) const {
46  return 1 + x * x + y * y + z * z * z;
47  }
48 };
49 
50 /**
51  * \brief Exact gradient
52  */
53 struct ExactFunctionGrad {
54  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
55  const double z) const {
57  grad(0) = 2 * x;
58  grad(1) = 2 * y;
59  grad(2) = 3 * z * z;
60  return grad;
61  }
62 };
63 
64 /**
65  * \brief Laplacian of function.
66  *
67  * This is Laplacian of \f$u\f$, it is calculated using formula
68  * \f[
69  * \nabla^2 u(x,y,z) = \nabla \cdot \nabla u
70  * \frac{\partial^2 u}{\partial x^2}+
71  * \frac{\partial^2 u}{\partial y^2}+
72  * \frac{\partial^2 u}{\partial z^2}
73  * \f]
74  *
75  */
77  double operator()(const double x, const double y, const double z) const {
78  return 4 + 6 * z;
79  }
80 };
81 
82 int main(int argc, char *argv[]) {
83 
84  // Initialize PETSc
85  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
86 
87  try {
88 
89  // Create MoAB database
90  moab::Core moab_core; // create database
91  moab::Interface &moab = moab_core; // create interface to database
92 
93  // Get command line options
94  int order = 3; // default approximation order
95  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
96  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
97  "none");
98  // Set approximation order
99  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
100  PETSC_NULL);
101  // Set testing (used by CTest)
102  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
103  &flg_test, PETSC_NULL);
104  ierr = PetscOptionsEnd();
105  CHKERRG(ierr);
106 
107  // Create MoFEM database and link it to MoAB
108  MoFEM::Core mofem_core(moab); // create database
109  MoFEM::Interface &m_field = mofem_core; // create interface to database
110  // Register DM Manager
111  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
112 
113  // Create vector to store approximation global error
114  Vec global_error;
115  CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
116 
117  // First we crate elements, implementation of elements is problem independent,
118  // we do not know yet what fields are present in the problem, or
119  // even we do not decided yet what approximation base or spaces we
120  // are going to use. Implementation of element is free from
121  // those constrains and can be used in different context.
122 
123  // Elements used by KSP & DM to assemble system of equations
124  boost::shared_ptr<ForcesAndSourcesCore> domain_lhs_fe; ///< Volume element for the matrix
125  boost::shared_ptr<ForcesAndSourcesCore> boundary_lhs_fe; ///< Boundary element for the matrix
126  boost::shared_ptr<ForcesAndSourcesCore> domain_rhs_fe; ///< Volume element to assemble vector
127  boost::shared_ptr<ForcesAndSourcesCore> boundary_rhs_fe; ///< Volume element to assemble vector
128  boost::shared_ptr<ForcesAndSourcesCore> domain_error; ///< Volume element evaluate error
129  boost::shared_ptr<ForcesAndSourcesCore> post_proc_volume; ///< Volume element to Post-process results
130  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
131  {
132  // Add problem specific operators the generic finite elements to calculate matrices and vectors.
135  ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
136  boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe);
137  // Add problem specific operators the generic finite elements to calculate error on elements and global error
138  // in H1 norm
141  global_error, domain_error);
142  // Post-process results
144  .creatFEToPostProcessResults(post_proc_volume);
145  }
146 
147  // Get simple interface is simplified version enabling quick and
148  // easy construction of problem.
149  Simple *simple_interface;
150  // Query interface and get pointer to Simple interface
151  CHKERR m_field.getInterface(simple_interface);
152 
153  // Build problem with simple interface
154  {
155 
156  // Get options for simple interface from command line
157  CHKERR simple_interface->getOptions();
158  // Load mesh file to database
159  CHKERR simple_interface->loadFile();
160 
161  // Add field on domain and boundary. Field is declared by space and base and rank. space
162  // can be H1. Hcurl, Hdiv and L2, base can be AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more,
163  // where rank is number of coefficients for dof.
164  //
165  // Simple interface assumes that there are four types of field; 1) defined
166  // on body domain, 2) fields defined on body boundary, 3) skeleton field defined
167  // on finite element skeleton. Finally data field is defined on body domain. Data field
168  // is not solved but used for post-process or to keep material parameters, etc. Here
169  // we using it to calculate approximation error on elements.
170 
171  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is used
172  // to construct base functions.
173  CHKERR simple_interface->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,1);
174  // Add Lagrange multiplier field on body boundary
175  CHKERR simple_interface->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,1);
176  // Add error (data) field, we need only L2 norm. Later order is set to 0, so this
177  // is piecewise discontinuous constant approx., i.e. 1 DOF for element. You can use
178  // more DOFs and collate moments of error to drive anisotropic h/p-adaptivity, however
179  // this is beyond this example.
180  CHKERR simple_interface->addDataField("ERROR",L2,AINSWORTH_LEGENDRE_BASE,1);
181 
182  // Set fields order domain and boundary fields.
183  CHKERR simple_interface->setFieldOrder("U",order); // to approximate function
184  CHKERR simple_interface->setFieldOrder("L",order); // to Lagrange multipliers
185  CHKERR simple_interface->setFieldOrder("ERROR",0); // approximation order for error
186 
187  // Setup problem. At that point database is constructed, i.e. fields, finite elements and
188  // problem data structures with local and global indexing.
189  CHKERR simple_interface->setUp();
190 
191  }
192 
193  // Get access to PETSC-MoFEM DM manager.
194  // or more derails see <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
195  // Form that point internal MoFEM data structures are linked with PETSc interface. If
196  // DM functions contains string MoFEM is is MoFEM specific DM interface function,
197  // otherwise DM function part of standard PETSc interface.
198 
199  DM dm;
200  // Get dm
201  CHKERR simple_interface->getDM(&dm);
202 
203  // Set KSP context for DM. At that point only elements are added to DM operators.
204  // Calculations of matrices and vectors is executed by KSP solver. This part
205  // of the code makes connection between implementation of finite elements and
206  // data operators with finite element declarations in DM manager. From that
207  // point DM takes responsibility for executing elements, calculations of
208  // matrices and vectors, and solution of the problem.
209  {
210  // Set operators for KSP solver
212  dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
214  dm, simple_interface->getBoundaryFEName(), boundary_lhs_fe, null,
215  null);
216  // Set calculation of the right hand side vector for KSP solver
217  CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getDomainFEName(),
218  domain_rhs_fe, null, null);
219  CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getBoundaryFEName(),
220  boundary_rhs_fe, null, null);
221  }
222 
223  // Solve problem, only PETEc interface here.
224  {
225 
226  // Create the right hand side vector and vector of unknowns
227  Vec F,D;
228  CHKERR DMCreateGlobalVector(dm,&F);
229  // Create unknown vector by creating duplicate copy of F vector. only
230  // structure is duplicated no values.
231  CHKERR VecDuplicate(F,&D);
232 
233  // Create solver and link it to DM
234  KSP solver;
235  CHKERR KSPCreate(PETSC_COMM_WORLD,&solver);
236  CHKERR KSPSetFromOptions(solver);
237  CHKERR KSPSetDM(solver,dm);
238  // Set-up solver, is type of solver and pre-conditioners
239  CHKERR KSPSetUp(solver);
240  // At solution process, KSP solver using DM creates matrices, Calculate
241  // values of the left hand side and the right hand side vector. then
242  // solves system of equations. Results are stored in vector D.
243  CHKERR KSPSolve(solver,F,D);
244 
245  // Scatter solution on the mesh. Stores unknown vector on field on the mesh.
246  CHKERR DMoFEMMeshToGlobalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE);
247 
248  // Clean data. Solver and vector are not needed any more.
249  CHKERR KSPDestroy(&solver);
250  CHKERR VecDestroy(&D);
251  CHKERR VecDestroy(&F);
252  }
253 
254  // Calculate error
255  {
256  // Loop over all elements in mesh, and run users operators on each element.
257  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
258  domain_error);
260  global_error);
261  CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
262  if (flg_test == PETSC_TRUE) {
263  CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
264  }
265  }
266 
267  {
268  // Loop over all elements in the mesh and for each execute post_proc_volume
269  // element and operators on it.
270  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
271  post_proc_volume);
272  // Write results
273  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
274  post_proc_volume)
275  ->writeFile("out_vol.h5m");
276  }
277 
278  // Destroy DM, no longer needed.
279  CHKERR DMDestroy(&dm);
280 
281  // Destroy ghost vector
282  CHKERR VecDestroy(&global_error);
283  }
284  CATCH_ERRORS;
285 
286  // finish work cleaning memory, getting statistics, etc.
288 
289  return 0;
290 
291 }
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:672
Deprecated interface functions.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:288
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:706
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: DMMMoFEM.cpp:556
double operator()(const double x, const double y, const double z) const
static char help[]
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
Core (interface) class.
Definition: Core.hpp:50
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: DMMMoFEM.cpp:597
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
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.
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:324
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
int main(int argc, char *argv[])
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:287
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:457
Create finite elements instances.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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:269
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.
#define CHKERR
Inline error check.
Definition: definitions.h:602
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:281
constexpr int order
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:482
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:67
field with C-1 continuity
Definition: definitions.h:180