libMesh
introduction_ex3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
7 // version 2.1 of the License, or (at your option) any later version.
8
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19 // <h1>Introduction Example 3 - Solving a Poisson Problem</h1>
20 // \author Benjamin S. Kirk
21 // \date 2003
22 //
23 // This is the third example program. It builds on
24 // the second example program by showing how to solve a simple
25 // Poisson system. This example also introduces the notion
26 // of customized matrix assembly functions, working with an
27 // exact solution, and using element iterators.
28 // We will not comment on things that
29 // were already explained in the second example.
30
31 // C++ include files that we need
32 #include <iostream>
33 #include <algorithm>
34 #include <math.h>
35
36 // Basic include files needed for the mesh functionality.
37 #include "libmesh/libmesh.h"
38 #include "libmesh/mesh.h"
39 #include "libmesh/mesh_generation.h"
40 #include "libmesh/vtk_io.h"
41 #include "libmesh/linear_implicit_system.h"
42 #include "libmesh/equation_systems.h"
43
44 // Define the Finite Element object.
45 #include "libmesh/fe.h"
46
47 // Define Gauss quadrature rules.
49
50 // Define useful datatypes for finite element
51 // matrix and vector components.
52 #include "libmesh/sparse_matrix.h"
53 #include "libmesh/numeric_vector.h"
54 #include "libmesh/dense_matrix.h"
55 #include "libmesh/dense_vector.h"
56 #include "libmesh/elem.h"
57
58 // Define the DofMap, which handles degree of freedom
59 // indexing.
60 #include "libmesh/dof_map.h"
61
62 // Bring in everything from the libMesh namespace
63 using namespace libMesh;
64
65 // Function prototype. This is the function that will assemble
66 // the linear system for our Poisson problem. Note that the
67 // function will take the EquationSystems object and the
68 // name of the system we are assembling as input. From the
70 // other objects we might need.
72  const std::string & system_name);
73
74 // Function prototype for the exact solution.
75 Real exact_solution (const Real x,
76  const Real y,
77  const Real z = 0.);
78
79 int main (int argc, char ** argv)
80 {
81  // Initialize libraries, like in example 2.
82  LibMeshInit init (argc, argv);
83
84  // This example requires a linear solver package.
85  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
86  "--enable-petsc, --enable-trilinos, or --enable-eigen");
87
88  // Brief message to the user regarding the program name
89  // and command line arguments.
90  libMesh::out << "Running " << argv[0];
91
92  for (int i=1; i<argc; i++)
93  libMesh::out << " " << argv[i];
94
95  libMesh::out << std::endl << std::endl;
96
97  // Skip this 2D example if libMesh was compiled as 1D-only.
98  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
99
100  // Create a mesh, with dimension to be overridden later, distributed
101  // across the default MPI communicator.
102  Mesh mesh(init.comm());
103
104  // Use the MeshTools::Generation mesh generator to create a uniform
105  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
106  // to build a mesh of 15x15 QUAD9 elements. Building QUAD9
107  // elements instead of the default QUAD4's we used in example 2
108  // allow us to use higher-order approximation.
110  15, 15,
111  -1., 1.,
112  -1., 1.,
114
115  // Print information about the mesh to the screen.
116  // Note that 5x5 QUAD9 elements actually has 11x11 nodes,
117  // so this mesh is significantly larger than the one in example 2.
118  mesh.print_info();
119
120  // Create an equation systems object.
121  EquationSystems equation_systems (mesh);
122
123  // Declare the Poisson system and its variables.
124  // The Poisson system is another example of a steady system.
126
127  // Adds the variable "u" to "Poisson". "u"
128  // will be approximated using second-order approximation.
130
131  // Give the system a pointer to the matrix assembly
132  // function. This will be called when needed by the
133  // library.
134  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
135
136  // Initialize the data structures for the equation system.
137  equation_systems.init();
138
139  // Prints information about the system to the screen.
140  equation_systems.print_info();
141
142  // Solve the system "Poisson". Note that calling this
143  // member will assemble the linear system and invoke
144  // the default numerical solver. With PETSc the solver can be
145  // controlled from the command line. For example,
146  // you can invoke conjugate gradient with:
147  //
148  // ./introduction_ex3 -ksp_type cg
149  //
150  // You can also get a nice X-window that monitors the solver
151  // convergence with:
152  //
153  // ./introduction-ex3 -ksp_xmonitor
154  //
155  // if you linked against the appropriate X libraries when you
156  // built PETSc.
157  equation_systems.get_system("Poisson").solve();
158
159 #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
160
161  // After solving the system write the solution
162  // to a VTK-formatted plot file.
163  VTKIO (mesh).write_equation_systems ("out.pvtu", equation_systems);
164
165 #endif // #ifdef LIBMESH_HAVE_VTK
166
167  // All done.
168  return 0;
169 }
170
171
172
173 // We now define the matrix assembly function for the
174 // Poisson system. We need to first compute element
175 // matrices and right-hand sides, and then take into
176 // account the boundary conditions, which will be handled
177 // via a penalty method.
179  const std::string & libmesh_dbg_var(system_name))
180 {
181
182  // It is a good idea to make sure we are assembling
183  // the proper system.
184  libmesh_assert_equal_to (system_name, "Poisson");
185
186  // Get a constant reference to the mesh object.
187  const MeshBase & mesh = es.get_mesh();
188
189  // The dimension that we are running
190  const unsigned int dim = mesh.mesh_dimension();
191
192  // Get a reference to the LinearImplicitSystem we are solving
193  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Poisson");
194
195  // A reference to the DofMap object for this system. The DofMap
196  // object handles the index translation from node and element numbers
197  // to degree of freedom numbers. We will talk more about the DofMap
198  // in future examples.
199  const DofMap & dof_map = system.get_dof_map();
200
201  // Get a constant reference to the Finite Element type
202  // for the first (and only) variable in the system.
203  FEType fe_type = dof_map.variable_type(0);
204
205  // Build a Finite Element object of the specified type. Since the
206  // FEBase::build() member dynamically creates memory we will
207  // store the object as a UniquePtr<FEBase>. This can be thought
208  // of as a pointer that will clean up after itself. Introduction Example 4
209  // describes some advantages of UniquePtr's in the context of
211  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
212
213  // A 5th order Gauss quadrature rule for numerical integration.
214  QGauss qrule (dim, FIFTH);
215
216  // Tell the finite element object to use our quadrature rule.
218
219  // Declare a special finite element object for
220  // boundary integration.
221  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
222
223  // Boundary integration requires one quadrature rule,
224  // with dimensionality one less than the dimensionality
225  // of the element.
226  QGauss qface(dim-1, FIFTH);
227
228  // Tell the finite element object to use our
231
232  // Here we define some references to cell-specific data that
233  // will be used to assemble the linear system.
234  //
235  // The element Jacobian * quadrature weight at each integration point.
236  const std::vector<Real> & JxW = fe->get_JxW();
237
238  // The physical XY locations of the quadrature points on the element.
239  // These might be useful for evaluating spatially varying material
240  // properties at the quadrature points.
241  const std::vector<Point> & q_point = fe->get_xyz();
242
243  // The element shape functions evaluated at the quadrature points.
244  const std::vector<std::vector<Real>> & phi = fe->get_phi();
245
247  // points.
248  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
249
250  // Define data structures to contain the element matrix
251  // and right-hand-side vector contribution. Following
252  // basic finite element terminology we will denote these
253  // "Ke" and "Fe". These datatypes are templated on
254  // Number, which allows the same code to work for real
255  // or complex numbers.
258
259  // This vector will hold the degree of freedom indices for
260  // the element. These define where in the global system
261  // the element degrees of freedom get mapped.
262  std::vector<dof_id_type> dof_indices;
263
264  // Now we will loop over all the elements in the mesh.
265  // We will compute the element matrix and right-hand-side
266  // contribution.
267  //
268  // Element iterators are a nice way to iterate through all the
269  // elements, or all the elements that have some property. The
270  // iterator el will iterate from the first to the last element on
271  // the local processor. The iterator end_el tells us when to stop.
272  // It is smart to make this one const so that we don't accidentally
273  // mess it up! In case users later modify this program to include
274  // refinement, we will be safe and will only consider the active
275  // elements; hence we use a variant of the active_elem_iterator.
276  for (const auto & elem : mesh.active_local_element_ptr_range())
277  {
278  // Get the degree of freedom indices for the
279  // current element. These define where in the global
280  // matrix and right-hand-side this element will
281  // contribute to.
282  dof_map.dof_indices (elem, dof_indices);
283
284  // Compute the element-specific data for the current
285  // element. This involves computing the location of the
286  // quadrature points (q_point) and the shape functions
287  // (phi, dphi) for the current element.
288  fe->reinit (elem);
289
290  // Zero the element matrix and right-hand side before
291  // summing them. We use the resize member here because
292  // the number of degrees of freedom might have changed from
293  // the last element. Note that this will be the case if the
294  // element type is different (i.e. the last element was a
295  // triangle, now we are on a quadrilateral).
296
297  // The DenseMatrix::resize() and the DenseVector::resize()
298  // members will automatically zero out the matrix and vector.
299  Ke.resize (dof_indices.size(),
300  dof_indices.size());
301
302  Fe.resize (dof_indices.size());
303
304  // Now loop over the quadrature points. This handles
305  // the numeric integration.
306  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
307  {
308
309  // Now we will build the element matrix. This involves
310  // a double loop to integrate the test functions (i) against
311  // the trial functions (j).
312  for (std::size_t i=0; i<phi.size(); i++)
313  for (std::size_t j=0; j<phi.size(); j++)
314  {
315  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
316  }
317
318  // This is the end of the matrix summation loop
319  // Now we build the element right-hand-side contribution.
320  // This involves a single loop in which we integrate the
321  // "forcing function" in the PDE against the test functions.
322  {
323  const Real x = q_point[qp](0);
324  const Real y = q_point[qp](1);
325  const Real eps = 1.e-3;
326
327
328  // "fxy" is the forcing function for the Poisson equation.
329  // In this case we set fxy to be a finite difference
330  // Laplacian approximation to the (known) exact solution.
331  //
332  // We will use the second-order accurate FD Laplacian
333  // approximation, which in 2D is
334  //
335  // u_xx + u_yy = (u(i,j-1) + u(i,j+1) +
336  // u(i-1,j) + u(i+1,j) +
337  // -4*u(i,j))/h^2
338  //
339  // Since the value of the forcing function depends only
340  // on the location of the quadrature point (q_point[qp])
341  // we will compute it here, outside of the i-loop
342  const Real fxy = -(exact_solution(x, y-eps) +
343  exact_solution(x, y+eps) +
344  exact_solution(x-eps, y) +
345  exact_solution(x+eps, y) -
346  4.*exact_solution(x, y))/eps/eps;
347
348  for (std::size_t i=0; i<phi.size(); i++)
349  Fe(i) += JxW[qp]*fxy*phi[i][qp];
350  }
351  }
352
353  // We have now reached the end of the RHS summation,
354  // and the end of quadrature point loop, so
355  // the interior element integration has
356  // been completed. However, we have not yet addressed
357  // boundary conditions. For this example we will only
358  // consider simple Dirichlet boundary conditions.
359  //
360  // There are several ways Dirichlet boundary conditions
361  // can be imposed. A simple approach, which works for
362  // interpolary bases like the standard Lagrange polynomials,
363  // is to assign function values to the
364  // degrees of freedom living on the domain boundary. This
365  // works well for interpolary bases, but is more difficult
366  // when non-interpolary (e.g Legendre or Hierarchic) bases
367  // are used.
368  //
369  // Dirichlet boundary conditions can also be imposed with a
370  // "penalty" method. In this case essentially the L2 projection
371  // of the boundary values are added to the matrix. The
372  // projection is multiplied by some large factor so that, in
373  // floating point arithmetic, the existing (smaller) entries
374  // in the matrix and right-hand-side are effectively ignored.
375  //
376  // This amounts to adding a term of the form (in latex notation)
377  //
378  // \frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i
379  //
380  // where
381  //
382  // \frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1
383  {
384
385  // The following loop is over the sides of the element.
386  // If the element has no neighbor on a side then that
387  // side MUST live on a boundary of the domain.
388  for (auto side : elem->side_index_range())
389  if (elem->neighbor_ptr(side) == libmesh_nullptr)
390  {
391  // The value of the shape functions at the quadrature
392  // points.
393  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
394
396  // points on the face.
397  const std::vector<Real> & JxW_face = fe_face->get_JxW();
398
399  // The XYZ locations (in physical space) of the
400  // quadrature points on the face. This is where
401  // we will interpolate the boundary value function.
402  const std::vector<Point> & qface_point = fe_face->get_xyz();
403
404  // Compute the shape function values on the element
405  // face.
406  fe_face->reinit(elem, side);
407
408  // Loop over the face quadrature points for integration.
409  for (unsigned int qp=0; qp<qface.n_points(); qp++)
410  {
411  // The location on the boundary of the current
413  const Real xf = qface_point[qp](0);
414  const Real yf = qface_point[qp](1);
415
416  // The penalty value. \frac{1}{\epsilon}
417  // in the discussion above.
418  const Real penalty = 1.e10;
419
420  // The boundary value.
421  const Real value = exact_solution(xf, yf);
422
423  // Matrix contribution of the L2 projection.
424  for (std::size_t i=0; i<phi_face.size(); i++)
425  for (std::size_t j=0; j<phi_face.size(); j++)
426  Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
427
428  // Right-hand-side contribution of the L2
429  // projection.
430  for (std::size_t i=0; i<phi_face.size(); i++)
431  Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
432  }
433  }
434  }
435
436  // We have now finished the quadrature point loop,
437  // and have therefore applied all the boundary conditions.
438
439  // If this assembly program were to be used on an adaptive mesh,
440  // we would have to apply any hanging node constraint equations
441  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
442
443  // The element matrix and right-hand-side are now built
444  // for this element. Add them to the global matrix and
445  // right-hand-side vector. The SparseMatrix::add_matrix()
446  // and NumericVector::add_vector() members do this for us.
449  }
450
451  // All done!
452 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
unsigned int dim
This class provides a specific system class.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
This is the MeshBase class.
Definition: mesh_base.h:68
int main(int argc, char **argv)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2030
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
Real exact_solution(const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void assemble_poisson(EquationSystems &es, const std::string &system_name)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
static const bool value
Definition: xdr_io.C:108
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
unsigned int n_points() const