libMesh
vector_fe_ex1.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
20 // <h1>Vector Finite Element Example 1 - Solving an uncoupled Poisson Problem</h1>
21 // \author Paul Bauman
22 // \date 2012
23 //
24 // This is the first vector FE example program. It builds on
25 // the introduction_ex3 example program by showing how to solve a simple
26 // uncoupled Poisson system using vector Lagrange elements.
27
28 // C++ include files that we need
29 #include <iostream>
30 #include <algorithm>
31 #include <math.h>
32
33 // Basic include files needed for the mesh functionality.
34 #include "libmesh/libmesh.h"
35 #include "libmesh/mesh.h"
36 #include "libmesh/mesh_generation.h"
37 #include "libmesh/linear_implicit_system.h"
38 #include "libmesh/equation_systems.h"
39 #include "libmesh/exodusII_io.h"
40 #include "libmesh/gmv_io.h"
41
42 // Define the Finite Element object.
43 #include "libmesh/fe.h"
44
45 // Define Gauss quadrature rules.
47
48 // Define useful datatypes for finite element
49 // matrix and vector components.
50 #include "libmesh/sparse_matrix.h"
51 #include "libmesh/numeric_vector.h"
52 #include "libmesh/dense_matrix.h"
53 #include "libmesh/dense_vector.h"
54 #include "libmesh/elem.h"
55
56 // Define the DofMap, which handles degree of freedom
57 // indexing.
58 #include "libmesh/dof_map.h"
59
60 // Bring in everything from the libMesh namespace
61 using namespace libMesh;
62
63 // Function prototype. This is the function that will assemble
64 // the linear system for our Poisson problem. Note that the
65 // function will take the EquationSystems object and the
66 // name of the system we are assembling as input. From the
68 // other objects we might need.
70  const std::string & system_name);
71
72 // Function prototype for the exact solution.
73 Real exact_solution (const int component,
74  const Real x,
75  const Real y,
76  const Real z = 0.);
77
78 int main (int argc, char ** argv)
79 {
80  // Initialize libraries.
81  LibMeshInit init (argc, argv);
82
83  // This example requires a linear solver package.
84  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
85  "--enable-petsc, --enable-trilinos, or --enable-eigen");
86
87  // Brief message to the user regarding the program name
88  // and command line arguments.
89  libMesh::out << "Running " << argv[0];
90
91  for (int i=1; i<argc; i++)
92  libMesh::out << " " << argv[i];
93
94  libMesh::out << std::endl << std::endl;
95
96  // Skip this 2D example if libMesh was compiled as 1D-only.
97  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
98
99  // Create a mesh, with dimension to be overridden later, on the
100  // default MPI communicator.
101  Mesh mesh(init.comm());
102
103  // Use the MeshTools::Generation mesh generator to create a uniform
104  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
105  // to build a mesh of 15x15 QUAD9 elements.
107  15, 15,
108  -1., 1.,
109  -1., 1.,
111
112  // Print information about the mesh to the screen.
113  mesh.print_info();
114
115  // Create an equation systems object.
116  EquationSystems equation_systems (mesh);
117
118  // Declare the Poisson system and its variables.
119  // The Poisson system is another example of a steady system.
121
122  // Adds the variable "u" to "Poisson". "u"
123  // will be approximated using second-order approximation
124  // using vector Lagrange elements. Since the mesh is 2-D, "u" will
125  // have two components.
127
128  // Give the system a pointer to the matrix assembly
129  // function. This will be called when needed by the
130  // library.
131  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
132
133  // Initialize the data structures for the equation system.
134  equation_systems.init();
135
136  // Prints information about the system to the screen.
137  equation_systems.print_info();
138
139  // Solve the system "Poisson". Note that calling this
140  // member will assemble the linear system and invoke
141  // the default numerical solver. With PETSc the solver can be
142  // controlled from the command line. For example,
143  // you can invoke conjugate gradient with:
144  //
145  // ./vector_fe_ex1 -ksp_type cg
146  //
147  // You can also get a nice X-window that monitors the solver
148  // convergence with:
149  //
150  // ./vector_fe_ex1 -ksp_xmonitor
151  //
152  // if you linked against the appropriate X libraries when you
153  // built PETSc.
154  equation_systems.get_system("Poisson").solve();
155
156 #ifdef LIBMESH_HAVE_EXODUS_API
157  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
158 #endif
159
160 #ifdef LIBMESH_HAVE_GMV
161  GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
162 #endif
163
164  // All done.
165  return 0;
166 }
167
168
169
170 // We now define the matrix assembly function for the
171 // Poisson system. We need to first compute element
172 // matrices and right-hand sides, and then take into
173 // account the boundary conditions, which will be handled
174 // via a penalty method.
176  const std::string & libmesh_dbg_var(system_name))
177 {
178
179  // It is a good idea to make sure we are assembling
180  // the proper system.
181  libmesh_assert_equal_to (system_name, "Poisson");
182
183  // Get a constant reference to the mesh object.
184  const MeshBase & mesh = es.get_mesh();
185
186  // The dimension that we are running
187  const unsigned int dim = mesh.mesh_dimension();
188
189  // Get a reference to the LinearImplicitSystem we are solving
190  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Poisson");
191
192  // A reference to the DofMap object for this system. The DofMap
193  // object handles the index translation from node and element numbers
194  // to degree of freedom numbers. We will talk more about the DofMap
195  // in future examples.
196  const DofMap & dof_map = system.get_dof_map();
197
198  // Get a constant reference to the Finite Element type
199  // for the first (and only) variable in the system.
200  FEType fe_type = dof_map.variable_type(0);
201
202  // Build a Finite Element object of the specified type.
203  // Note that FEVectorBase is a typedef for the templated FE
204  // class.
205  UniquePtr<FEVectorBase> fe (FEVectorBase::build(dim, fe_type));
206
207  // A 5th order Gauss quadrature rule for numerical integration.
208  QGauss qrule (dim, FIFTH);
209
210  // Tell the finite element object to use our quadrature rule.
212
213  // Declare a special finite element object for
214  // boundary integration.
215  UniquePtr<FEVectorBase> fe_face (FEVectorBase::build(dim, fe_type));
216
217  // Boundary integration requires one quadrature rule,
218  // with dimensionality one less than the dimensionality
219  // of the element.
220  QGauss qface(dim-1, FIFTH);
221
222  // Tell the finite element object to use our
225
226  // Here we define some references to cell-specific data that
227  // will be used to assemble the linear system.
228  //
229  // The element Jacobian * quadrature weight at each integration point.
230  const std::vector<Real> & JxW = fe->get_JxW();
231
232  // The physical XY locations of the quadrature points on the element.
233  // These might be useful for evaluating spatially varying material
234  // properties at the quadrature points.
235  const std::vector<Point> & q_point = fe->get_xyz();
236
237  // The element shape functions evaluated at the quadrature points.
238  // Notice the shape functions are a vector rather than a scalar.
239  const std::vector<std::vector<RealGradient>> & phi = fe->get_phi();
240
242  // points. Notice that the shape function gradients are a tensor.
243  const std::vector<std::vector<RealTensor>> & dphi = fe->get_dphi();
244
245  // Define data structures to contain the element matrix
246  // and right-hand-side vector contribution. Following
247  // basic finite element terminology we will denote these
248  // "Ke" and "Fe". These datatypes are templated on
249  // Number, which allows the same code to work for real
250  // or complex numbers.
253
254  // This vector will hold the degree of freedom indices for
255  // the element. These define where in the global system
256  // the element degrees of freedom get mapped.
257  std::vector<dof_id_type> dof_indices;
258
259  // Now we will loop over all the elements in the mesh.
260  // We will compute the element matrix and right-hand-side
261  // contribution.
262  //
263  // Element iterators are a nice way to iterate through all the
264  // elements, or all the elements that have some property. The
265  // iterator el will iterate from the first to the last element on
266  // the local processor. The iterator end_el tells us when to stop.
267  // It is smart to make this one const so that we don't accidentally
268  // mess it up! In case users later modify this program to include
269  // refinement, we will be safe and will only consider the active
270  // elements; hence we use a variant of the active_elem_iterator.
271  for (const auto & elem : mesh.active_local_element_ptr_range())
272  {
273  // Get the degree of freedom indices for the
274  // current element. These define where in the global
275  // matrix and right-hand-side this element will
276  // contribute to.
277  dof_map.dof_indices (elem, dof_indices);
278
279  // Compute the element-specific data for the current
280  // element. This involves computing the location of the
281  // quadrature points (q_point) and the shape functions
282  // (phi, dphi) for the current element.
283  fe->reinit (elem);
284
285  // Zero the element matrix and right-hand side before
286  // summing them. We use the resize member here because
287  // the number of degrees of freedom might have changed from
288  // the last element. Note that this will be the case if the
289  // element type is different (i.e. the last element was a
290  // triangle, now we are on a quadrilateral).
291
292  // The DenseMatrix::resize() and the DenseVector::resize()
293  // members will automatically zero out the matrix and vector.
294  Ke.resize (dof_indices.size(),
295  dof_indices.size());
296
297  Fe.resize (dof_indices.size());
298
299  // Now loop over the quadrature points. This handles
300  // the numeric integration.
301  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
302  {
303  // Now we will build the element matrix. This involves
304  // a double loop to integrate the test functions (i) against
305  // the trial functions (j).
306  for (std::size_t i=0; i<phi.size(); i++)
307  for (std::size_t j=0; j<phi.size(); j++)
308  Ke(i,j) += JxW[qp] * dphi[i][qp].contract(dphi[j][qp]);
309
310  // This is the end of the matrix summation loop
311  // Now we build the element right-hand-side contribution.
312  // This involves a single loop in which we integrate the
313  // "forcing function" in the PDE against the test functions.
314  {
315  const Real x = q_point[qp](0);
316  const Real y = q_point[qp](1);
317  const Real eps = 1.e-3;
318
319  // "f" is the forcing function for the Poisson equation.
320  // In this case we set f to be a finite difference
321  // Laplacian approximation to the (known) exact solution.
322  //
323  // We will use the second-order accurate FD Laplacian
324  // approximation, which in 2D is
325  //
326  // u_xx + u_yy = (u(i,j-1) + u(i,j+1) +
327  // u(i-1,j) + u(i+1,j) +
328  // -4*u(i,j))/h^2
329  //
330  // Since the value of the forcing function depends only
331  // on the location of the quadrature point (q_point[qp])
332  // we will compute it here, outside of the i-loop
333  const Real fx = -(exact_solution(0, x, y-eps) +
334  exact_solution(0, x, y+eps) +
335  exact_solution(0, x-eps, y) +
336  exact_solution(0, x+eps, y) -
337  4.*exact_solution(0, x, y))/eps/eps;
338
339  const Real fy = -(exact_solution(1, x, y-eps) +
340  exact_solution(1, x, y+eps) +
341  exact_solution(1, x-eps, y) +
342  exact_solution(1, x+eps, y) -
343  4.*exact_solution(1, x, y))/eps/eps;
344
346
347  for (std::size_t i=0; i<phi.size(); i++)
348  Fe(i) += JxW[qp]*f*phi[i][qp];
349  }
350  }
351
352  // We have now reached the end of the RHS summation,
353  // and the end of quadrature point loop, so
354  // the interior element integration has
355  // been completed. However, we have not yet addressed
356  // boundary conditions. For this example we will only
357  // consider simple Dirichlet boundary conditions.
358  //
359  // There are several ways Dirichlet boundary conditions
360  // can be imposed. A simple approach, which works for
361  // interpolary bases like the standard Lagrange polynomials,
362  // is to assign function values to the
363  // degrees of freedom living on the domain boundary. This
364  // works well for interpolary bases, but is more difficult
365  // when non-interpolary (e.g Legendre or Hierarchic) bases
366  // are used.
367  //
368  // Dirichlet boundary conditions can also be imposed with a
369  // "penalty" method. In this case essentially the L2 projection
370  // of the boundary values are added to the matrix. The
371  // projection is multiplied by some large factor so that, in
372  // floating point arithmetic, the existing (smaller) entries
373  // in the matrix and right-hand-side are effectively ignored.
374  //
375  // This amounts to adding a term of the form (in latex notation)
376  //
377  // \frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i
378  //
379  // where
380  //
381  // \frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1
382  {
383  // The following loop is over the sides of the element.
384  // If the element has no neighbor on a side then that
385  // side MUST live on a boundary of the domain.
386  for (auto side : elem->side_index_range())
387  if (elem->neighbor_ptr(side) == libmesh_nullptr)
388  {
389  // The value of the shape functions at the quadrature
390  // points.
391  const std::vector<std::vector<RealGradient>> & phi_face = fe_face->get_phi();
392
394  // points on the face.
395  const std::vector<Real> & JxW_face = fe_face->get_JxW();
396
397  // The XYZ locations (in physical space) of the
398  // quadrature points on the face. This is where
399  // we will interpolate the boundary value function.
400  const std::vector<Point> & qface_point = fe_face->get_xyz();
401
402  // Compute the shape function values on the element
403  // face.
404  fe_face->reinit(elem, side);
405
406  // Loop over the face quadrature points for integration.
407  for (unsigned int qp=0; qp<qface.n_points(); qp++)
408  {
409  // The location on the boundary of the current
411  const Real xf = qface_point[qp](0);
412  const Real yf = qface_point[qp](1);
413
414  // The penalty value. \frac{1}{\epsilon}
415  // in the discussion above.
416  const Real penalty = 1.e10;
417
418  // The boundary values.
419  const RealGradient f(exact_solution(0, xf, yf),
420  exact_solution(1, xf, yf));
421
422  // Matrix contribution of the L2 projection.
423  for (std::size_t i=0; i<phi_face.size(); i++)
424  for (std::size_t j=0; j<phi_face.size(); j++)
425  Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
426
427  // Right-hand-side contribution of the L2
428  // projection.
429  for (std::size_t i=0; i<phi_face.size(); i++)
430  Fe(i) += JxW_face[qp]*penalty*f*phi_face[i][qp];
431  }
432  }
433  }
434
435  // We have now finished the quadrature point loop,
436  // and have therefore applied all the boundary conditions.
437
438  // If this assembly program were to be used on an adaptive mesh,
439  // we would have to apply any hanging node constraint equations
440  //dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
441
442  // The element matrix and right-hand-side are now built
443  // for this element. Add them to the global matrix and
444  // right-hand-side vector. The SparseMatrix::add_matrix()
445  // and NumericVector::add_vector() members do this for us.
448  }
449
450  // All done!
451 }
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
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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 writing meshes in the GMV format.
Definition: gmv_io.h:46
This is the MeshBase class.
Definition: mesh_base.h:68
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)
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
int main(int argc, char **argv)
Definition: vector_fe_ex1.C:78
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
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
Real exact_solution(const int component, const Real x, const Real y, const Real z=0.)
This is the exact solution that we are trying to obtain.
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