libMesh
systems_of_equations_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
6 // License as published by the Free Software Foundation; either
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>Systems Example 1 - Stokes Equations</h1>
21 // \author Benjamin S. Kirk
22 // \date 2003
23 //
24 // This example shows how a simple, linear system of equations
25 // can be solved in parallel. The system of equations are the familiar
26 // Stokes equations for low-speed incompressible fluid flow.
27 
28 // C++ include files that we need
29 #include <iostream>
30 #include <algorithm>
31 #include <math.h>
32 
33 // Basic include file needed for the mesh functionality.
34 #include "libmesh/libmesh.h"
35 #include "libmesh/mesh.h"
36 #include "libmesh/mesh_generation.h"
37 #include "libmesh/exodusII_io.h"
38 #include "libmesh/equation_systems.h"
39 #include "libmesh/fe.h"
40 #include "libmesh/quadrature_gauss.h"
41 #include "libmesh/dof_map.h"
42 #include "libmesh/sparse_matrix.h"
43 #include "libmesh/numeric_vector.h"
44 #include "libmesh/dense_matrix.h"
45 #include "libmesh/dense_vector.h"
46 #include "libmesh/linear_implicit_system.h"
47 
48 // For systems of equations the DenseSubMatrix
49 // and DenseSubVector provide convenient ways for
50 // assembling the element matrix and vector on a
51 // component-by-component basis.
52 #include "libmesh/dense_submatrix.h"
53 #include "libmesh/dense_subvector.h"
54 
55 // The definition of a geometric element
56 #include "libmesh/elem.h"
57 
58 // Bring in everything from the libMesh namespace
59 using namespace libMesh;
60 
61 // Function prototype. This function will assemble the system
62 // matrix and right-hand-side.
64  const std::string & system_name);
65 
66 // The main program.
67 int main (int argc, char ** argv)
68 {
69  // Initialize libMesh.
70  LibMeshInit init (argc, argv);
71 
72  // This example requires a linear solver package.
73  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
74  "--enable-petsc, --enable-trilinos, or --enable-eigen");
75 
76  // Skip this 2D example if libMesh was compiled as 1D-only.
77  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
78 
79  // This example NaNs with the Eigen sparse linear solvers and
80  // Trilinos solvers, but should work OK with either PETSc or
81  // Laspack.
82  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");
83  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");
84 
85  // Create a mesh, with dimension to be overridden later, distributed
86  // across the default MPI communicator.
87  Mesh mesh(init.comm());
88 
89  // Use the MeshTools::Generation mesh generator to create a uniform
90  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
91  // to build a mesh of 8x8 Quad9 elements. Building these
92  // higher-order elements allows us to use higher-order
93  // approximation, as in example 3.
95  15, 15,
96  0., 1.,
97  0., 1.,
98  QUAD9);
99 
100  // Print information about the mesh to the screen.
101  mesh.print_info();
102 
103  // Create an equation systems object.
104  EquationSystems equation_systems (mesh);
105 
106  // Declare the system and its variables.
107  // Create a transient system named "Stokes"
108  LinearImplicitSystem & system =
109  equation_systems.add_system<LinearImplicitSystem> ("Stokes");
110 
111  // Add the variables "u" & "v" to "Stokes". They
112  // will be approximated using second-order approximation.
113  system.add_variable ("u", SECOND);
114  system.add_variable ("v", SECOND);
115 
116  // Add the variable "p" to "Stokes". This will
117  // be approximated with a first-order basis,
118  // providing an LBB-stable pressure-velocity pair.
119  system.add_variable ("p", FIRST);
120 
121  // Give the system a pointer to the matrix assembly
122  // function.
124 
125  // Initialize the data structures for the equation system.
126  equation_systems.init ();
127 
128  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
129  equation_systems.parameters.set<Real> ("linear solver tolerance") = TOLERANCE;
130 
131  // Prints information about the system to the screen.
132  equation_systems.print_info();
133 
134  // Assemble & solve the linear system,
135  // then write the solution.
136  equation_systems.get_system("Stokes").solve();
137 
138 #ifdef LIBMESH_HAVE_EXODUS_API
140  equation_systems);
141 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
142 
143  // All done.
144  return 0;
145 }
146 
148  const std::string & libmesh_dbg_var(system_name))
149 {
150  // It is a good idea to make sure we are assembling
151  // the proper system.
152  libmesh_assert_equal_to (system_name, "Stokes");
153 
154  // Get a constant reference to the mesh object.
155  const MeshBase & mesh = es.get_mesh();
156 
157  // The dimension that we are running
158  const unsigned int dim = mesh.mesh_dimension();
159 
160  // Get a reference to the Convection-Diffusion system object.
161  LinearImplicitSystem & system =
162  es.get_system<LinearImplicitSystem> ("Stokes");
163 
164  // Numeric ids corresponding to each variable in the system
165  const unsigned int u_var = system.variable_number ("u");
166  const unsigned int v_var = system.variable_number ("v");
167  const unsigned int p_var = system.variable_number ("p");
168 
169  // Get the Finite Element type for "u". Note this will be
170  // the same as the type for "v".
171  FEType fe_vel_type = system.variable_type(u_var);
172 
173  // Get the Finite Element type for "p".
174  FEType fe_pres_type = system.variable_type(p_var);
175 
176  // Build a Finite Element object of the specified type for
177  // the velocity variables.
178  UniquePtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
179 
180  // Build a Finite Element object of the specified type for
181  // the pressure variables.
182  UniquePtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
183 
184  // A Gauss quadrature rule for numerical integration.
185  // Let the FEType object decide what order rule is appropriate.
186  QGauss qrule (dim, fe_vel_type.default_quadrature_order());
187 
188  // Tell the finite element objects to use our quadrature rule.
189  fe_vel->attach_quadrature_rule (&qrule);
190  fe_pres->attach_quadrature_rule (&qrule);
191 
192  // Here we define some references to cell-specific data that
193  // will be used to assemble the linear system.
194  //
195  // The element Jacobian * quadrature weight at each integration point.
196  const std::vector<Real> & JxW = fe_vel->get_JxW();
197 
198  // The element shape function gradients for the velocity
199  // variables evaluated at the quadrature points.
200  const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
201 
202  // The element shape functions for the pressure variable
203  // evaluated at the quadrature points.
204  const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
205 
206  // A reference to the DofMap object for this system. The DofMap
207  // object handles the index translation from node and element numbers
208  // to degree of freedom numbers. We will talk more about the DofMap
209  // in future examples.
210  const DofMap & dof_map = system.get_dof_map();
211 
212  // Define data structures to contain the element matrix
213  // and right-hand-side vector contribution. Following
214  // basic finite element terminology we will denote these
215  // "Ke" and "Fe".
218 
220  Kuu(Ke), Kuv(Ke), Kup(Ke),
221  Kvu(Ke), Kvv(Ke), Kvp(Ke),
222  Kpu(Ke), Kpv(Ke), Kpp(Ke);
223 
225  Fu(Fe),
226  Fv(Fe),
227  Fp(Fe);
228 
229  // This vector will hold the degree of freedom indices for
230  // the element. These define where in the global system
231  // the element degrees of freedom get mapped.
232  std::vector<dof_id_type> dof_indices;
233  std::vector<dof_id_type> dof_indices_u;
234  std::vector<dof_id_type> dof_indices_v;
235  std::vector<dof_id_type> dof_indices_p;
236 
237  // Now we will loop over all the elements in the mesh that
238  // live on the local processor. We will compute the element
239  // matrix and right-hand-side contribution. In case users later
240  // modify this program to include refinement, we will be safe and
241  // will only consider the active elements; hence we use a variant of
242  // the active_elem_iterator.
243  for (const auto & elem : mesh.active_local_element_ptr_range())
244  {
245  // Get the degree of freedom indices for the
246  // current element. These define where in the global
247  // matrix and right-hand-side this element will
248  // contribute to.
249  dof_map.dof_indices (elem, dof_indices);
250  dof_map.dof_indices (elem, dof_indices_u, u_var);
251  dof_map.dof_indices (elem, dof_indices_v, v_var);
252  dof_map.dof_indices (elem, dof_indices_p, p_var);
253 
254  const unsigned int n_dofs = dof_indices.size();
255  const unsigned int n_u_dofs = dof_indices_u.size();
256  const unsigned int n_v_dofs = dof_indices_v.size();
257  const unsigned int n_p_dofs = dof_indices_p.size();
258 
259  // Compute the element-specific data for the current
260  // element. This involves computing the location of the
261  // quadrature points (q_point) and the shape functions
262  // (phi, dphi) for the current element.
263  fe_vel->reinit (elem);
264  fe_pres->reinit (elem);
265 
266  // Zero the element matrix and right-hand side before
267  // summing them. We use the resize member here because
268  // the number of degrees of freedom might have changed from
269  // the last element. Note that this will be the case if the
270  // element type is different (i.e. the last element was a
271  // triangle, now we are on a quadrilateral).
272  Ke.resize (n_dofs, n_dofs);
273  Fe.resize (n_dofs);
274 
275  // Reposition the submatrices... The idea is this:
276  //
277  // - - - -
278  // | Kuu Kuv Kup | | Fu |
279  // Ke = | Kvu Kvv Kvp |; Fe = | Fv |
280  // | Kpu Kpv Kpp | | Fp |
281  // - - - -
282  //
283  // The DenseSubMatrix.reposition () member takes the
284  // (row_offset, column_offset, row_size, column_size).
285  //
286  // Similarly, the DenseSubVector.reposition () member
287  // takes the (row_offset, row_size)
288  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
289  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
290  Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
291 
292  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
293  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
294  Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
295 
296  Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
297  Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
298  Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
299 
300  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
301  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
302  Fp.reposition (p_var*n_u_dofs, n_p_dofs);
303 
304  // Now we will build the element matrix.
305  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
306  {
307  // Assemble the u-velocity row
308  // uu coupling
309  for (unsigned int i=0; i<n_u_dofs; i++)
310  for (unsigned int j=0; j<n_u_dofs; j++)
311  Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
312 
313  // up coupling
314  for (unsigned int i=0; i<n_u_dofs; i++)
315  for (unsigned int j=0; j<n_p_dofs; j++)
316  Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0);
317 
318 
319  // Assemble the v-velocity row
320  // vv coupling
321  for (unsigned int i=0; i<n_v_dofs; i++)
322  for (unsigned int j=0; j<n_v_dofs; j++)
323  Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
324 
325  // vp coupling
326  for (unsigned int i=0; i<n_v_dofs; i++)
327  for (unsigned int j=0; j<n_p_dofs; j++)
328  Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1);
329 
330 
331  // Assemble the pressure row
332  // pu coupling
333  for (unsigned int i=0; i<n_p_dofs; i++)
334  for (unsigned int j=0; j<n_u_dofs; j++)
335  Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0);
336 
337  // pv coupling
338  for (unsigned int i=0; i<n_p_dofs; i++)
339  for (unsigned int j=0; j<n_v_dofs; j++)
340  Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1);
341 
342  } // end of the quadrature point qp-loop
343 
344  // At this point the interior element integration has
345  // been completed. However, we have not yet addressed
346  // boundary conditions. For this example we will only
347  // consider simple Dirichlet boundary conditions imposed
348  // via the penalty method. The penalty method used here
349  // is equivalent (for Lagrange basis functions) to lumping
350  // the matrix resulting from the L2 projection penalty
351  // approach introduced in example 3.
352  {
353  // The following loops over the sides of the element.
354  // If the element has no neighbor on a side then that
355  // side MUST live on a boundary of the domain.
356  for (auto s : elem->side_index_range())
357  if (elem->neighbor_ptr(s) == libmesh_nullptr)
358  {
359  UniquePtr<const Elem> side (elem->build_side_ptr(s));
360 
361  // Loop over the nodes on the side.
362  for (auto ns : side->node_index_range())
363  {
364  // The location on the boundary of the current
365  // node.
366 
367  // const Real xf = side->point(ns)(0);
368  const Real yf = side->point(ns)(1);
369 
370  // The penalty value. \f$ \frac{1}{\epsilon \f$
371  const Real penalty = 1.e10;
372 
373  // The boundary values.
374 
375  // Set u = 1 on the top boundary, 0 everywhere else
376  const Real u_value = (yf > .99) ? 1. : 0.;
377 
378  // Set v = 0 everywhere
379  const Real v_value = 0.;
380 
381  // Find the node on the element matching this node on
382  // the side. That defined where in the element matrix
383  // the boundary condition will be applied.
384  for (auto n : elem->node_index_range())
385  if (elem->node_id(n) == side->node_id(ns))
386  {
387  // Matrix contribution.
388  Kuu(n,n) += penalty;
389  Kvv(n,n) += penalty;
390 
391  // Right-hand-side contribution.
392  Fu(n) += penalty*u_value;
393  Fv(n) += penalty*v_value;
394  }
395  } // end face node loop
396  } // end if (elem->neighbor(side) == libmesh_nullptr)
397  } // end boundary condition section
398 
399  // If this assembly program were to be used on an adaptive mesh,
400  // we would have to apply any hanging node constraint equations.
401  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
402 
403  // The element matrix and right-hand-side are now built
404  // for this element. Add them to the global matrix and
405  // right-hand-side vector. The NumericMatrix::add_matrix()
406  // and NumericVector::add_vector() members do this for us.
407  system.matrix->add_matrix (Ke, dof_indices);
408  system.rhs->add_vector (Fe, dof_indices);
409  } // end of element loop
410 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
void assemble_stokes(EquationSystems &es, const std::string &system_name)
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.
int main(int argc, char **argv)
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 int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
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.
static const Real TOLERANCE
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.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
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
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
void reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
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.
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 attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
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.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void reposition(const unsigned int ioff, const unsigned int joff, const unsigned int new_m, const unsigned int new_n)
Changes the location of the submatrix in the parent matrix.
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
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.