Loading [MathJax]/jax/output/HTML-CSS/config.js
libMesh
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
Functions
systems_of_equations_ex1.C File Reference

Go to the source code of this file.

Functions

void assemble_stokes (EquationSystems &es, const std::string &system_name)
 
int main (int argc, char **argv)
 
void assemble_stokes (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 

Function Documentation

◆ assemble_stokes() [1/2]

void assemble_stokes ( EquationSystems es,
const std::string &  system_name 
)

Referenced by main().

◆ assemble_stokes() [2/2]

void assemble_stokes ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 151 of file systems_of_equations_ex1.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::Elem::node_id(), libMesh::Elem::node_index_range(), libMesh::Elem::point(), libMesh::Real, libMesh::DenseSubVector< T >::reposition(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::System::variable_number(), and libMesh::System::variable_type().

153 {
154  // It is a good idea to make sure we are assembling
155  // the proper system.
156  libmesh_assert_equal_to (system_name, "Stokes");
157 
158  // Get a constant reference to the mesh object.
159  const MeshBase & mesh = es.get_mesh();
160 
161  // The dimension that we are running
162  const unsigned int dim = mesh.mesh_dimension();
163 
164  // Get a reference to the Convection-Diffusion system object.
165  LinearImplicitSystem & system =
166  es.get_system<LinearImplicitSystem> ("Stokes");
167 
168  // Numeric ids corresponding to each variable in the system
169  const unsigned int u_var = system.variable_number ("u");
170  const unsigned int v_var = system.variable_number ("v");
171  const unsigned int p_var = system.variable_number ("p");
172 
173  // Get the Finite Element type for "u". Note this will be
174  // the same as the type for "v".
175  FEType fe_vel_type = system.variable_type(u_var);
176 
177  // Get the Finite Element type for "p".
178  FEType fe_pres_type = system.variable_type(p_var);
179 
180  // Build a Finite Element object of the specified type for
181  // the velocity variables.
182  std::unique_ptr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
183 
184  // Build a Finite Element object of the specified type for
185  // the pressure variables.
186  std::unique_ptr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
187 
188  // A Gauss quadrature rule for numerical integration.
189  // Let the FEType object decide what order rule is appropriate.
190  QGauss qrule (dim, fe_vel_type.default_quadrature_order());
191 
192  // Tell the finite element objects to use our quadrature rule.
193  fe_vel->attach_quadrature_rule (&qrule);
194  fe_pres->attach_quadrature_rule (&qrule);
195 
196  // Here we define some references to cell-specific data that
197  // will be used to assemble the linear system.
198  //
199  // The element Jacobian * quadrature weight at each integration point.
200  const std::vector<Real> & JxW = fe_vel->get_JxW();
201 
202  // The element shape function gradients for the velocity
203  // variables evaluated at the quadrature points.
204  const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
205 
206  // The element shape functions for the pressure variable
207  // evaluated at the quadrature points.
208  const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
209 
210  // A reference to the DofMap object for this system. The DofMap
211  // object handles the index translation from node and element numbers
212  // to degree of freedom numbers. We will talk more about the DofMap
213  // in future examples.
214  const DofMap & dof_map = system.get_dof_map();
215 
216  // Define data structures to contain the element matrix
217  // and right-hand-side vector contribution. Following
218  // basic finite element terminology we will denote these
219  // "Ke" and "Fe".
222 
224  Kuu(Ke), Kuv(Ke), Kup(Ke),
225  Kvu(Ke), Kvv(Ke), Kvp(Ke),
226  Kpu(Ke), Kpv(Ke), Kpp(Ke);
227 
229  Fu(Fe),
230  Fv(Fe),
231  Fp(Fe);
232 
233  // This vector will hold the degree of freedom indices for
234  // the element. These define where in the global system
235  // the element degrees of freedom get mapped.
236  std::vector<dof_id_type> dof_indices;
237  std::vector<dof_id_type> dof_indices_u;
238  std::vector<dof_id_type> dof_indices_v;
239  std::vector<dof_id_type> dof_indices_p;
240 
241  SparseMatrix<Number> & matrix = system.get_system_matrix();
242 
243  // Now we will loop over all the elements in the mesh that
244  // live on the local processor. We will compute the element
245  // matrix and right-hand-side contribution. In case users later
246  // modify this program to include refinement, we will be safe and
247  // will only consider the active elements; hence we use a variant of
248  // the active_elem_iterator.
249  for (const auto & elem : mesh.active_local_element_ptr_range())
250  {
251  // Get the degree of freedom indices for the
252  // current element. These define where in the global
253  // matrix and right-hand-side this element will
254  // contribute to.
255  dof_map.dof_indices (elem, dof_indices);
256  dof_map.dof_indices (elem, dof_indices_u, u_var);
257  dof_map.dof_indices (elem, dof_indices_v, v_var);
258  dof_map.dof_indices (elem, dof_indices_p, p_var);
259 
260  const unsigned int n_dofs = dof_indices.size();
261  const unsigned int n_u_dofs = dof_indices_u.size();
262  const unsigned int n_v_dofs = dof_indices_v.size();
263  const unsigned int n_p_dofs = dof_indices_p.size();
264 
265  // Compute the element-specific data for the current
266  // element. This involves computing the location of the
267  // quadrature points (q_point) and the shape functions
268  // (phi, dphi) for the current element.
269  fe_vel->reinit (elem);
270  fe_pres->reinit (elem);
271 
272  // Zero the element matrix and right-hand side before
273  // summing them. We use the resize member here because
274  // the number of degrees of freedom might have changed from
275  // the last element. Note that this will be the case if the
276  // element type is different (i.e. the last element was a
277  // triangle, now we are on a quadrilateral).
278  Ke.resize (n_dofs, n_dofs);
279  Fe.resize (n_dofs);
280 
281  // Reposition the submatrices... The idea is this:
282  //
283  // - - - -
284  // | Kuu Kuv Kup | | Fu |
285  // Ke = | Kvu Kvv Kvp |; Fe = | Fv |
286  // | Kpu Kpv Kpp | | Fp |
287  // - - - -
288  //
289  // The DenseSubMatrix.reposition () member takes the
290  // (row_offset, column_offset, row_size, column_size).
291  //
292  // Similarly, the DenseSubVector.reposition () member
293  // takes the (row_offset, row_size)
294  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
295  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
296  Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
297 
298  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
299  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
300  Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
301 
302  Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
303  Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
304  Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
305 
306  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
307  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
308  Fp.reposition (p_var*n_u_dofs, n_p_dofs);
309 
310  // Now we will build the element matrix.
311  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
312  {
313  // Assemble the u-velocity row
314  // uu coupling
315  for (unsigned int i=0; i<n_u_dofs; i++)
316  for (unsigned int j=0; j<n_u_dofs; j++)
317  Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
318 
319  // up coupling
320  for (unsigned int i=0; i<n_u_dofs; i++)
321  for (unsigned int j=0; j<n_p_dofs; j++)
322  Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0);
323 
324 
325  // Assemble the v-velocity row
326  // vv coupling
327  for (unsigned int i=0; i<n_v_dofs; i++)
328  for (unsigned int j=0; j<n_v_dofs; j++)
329  Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
330 
331  // vp coupling
332  for (unsigned int i=0; i<n_v_dofs; i++)
333  for (unsigned int j=0; j<n_p_dofs; j++)
334  Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1);
335 
336 
337  // Assemble the pressure row
338  // pu coupling
339  for (unsigned int i=0; i<n_p_dofs; i++)
340  for (unsigned int j=0; j<n_u_dofs; j++)
341  Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0);
342 
343  // pv coupling
344  for (unsigned int i=0; i<n_p_dofs; i++)
345  for (unsigned int j=0; j<n_v_dofs; j++)
346  Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1);
347 
348  } // end of the quadrature point qp-loop
349 
350  // At this point the interior element integration has
351  // been completed. However, we have not yet addressed
352  // boundary conditions. For this example we will only
353  // consider simple Dirichlet boundary conditions imposed
354  // via the penalty method. The penalty method used here
355  // is equivalent (for Lagrange basis functions) to lumping
356  // the matrix resulting from the L2 projection penalty
357  // approach introduced in example 3.
358  {
359  // To avoid extraneous memory allocation when building element sides
360  ElemSideBuilder side_builder;
361 
362  // The following loops over the sides of the element.
363  // If the element has no neighbor on a side then that
364  // side MUST live on a boundary of the domain.
365  for (auto s : elem->side_index_range())
366  if (elem->neighbor_ptr(s) == nullptr)
367  {
368  const Elem & side = side_builder(*elem, s);
369 
370  // Loop over the nodes on the side.
371  for (auto ns : side.node_index_range())
372  {
373  // The location on the boundary of the current
374  // node.
375 
376  // const Real xf = side->point(ns)(0);
377  const Real yf = side.point(ns)(1);
378 
379  // The penalty value. \f$ \frac{1}{\epsilon \f$
380  const Real penalty = 1.e10;
381 
382  // The boundary values.
383 
384  // Set u = 1 on the top boundary, 0 everywhere else
385  const Real u_value = (yf > .99) ? 1. : 0.;
386 
387  // Set v = 0 everywhere
388  const Real v_value = 0.;
389 
390  // Find the node on the element matching this node on
391  // the side. That defined where in the element matrix
392  // the boundary condition will be applied.
393  for (auto n : elem->node_index_range())
394  if (elem->node_id(n) == side.node_id(ns))
395  {
396  // Matrix contribution.
397  Kuu(n,n) += penalty;
398  Kvv(n,n) += penalty;
399 
400  // Right-hand-side contribution.
401  Fu(n) += penalty*u_value;
402  Fv(n) += penalty*v_value;
403  }
404  } // end face node loop
405  } // end if (elem->neighbor(side) == nullptr)
406  } // end boundary condition section
407 
408  // If this assembly program were to be used on an adaptive mesh,
409  // we would have to apply any hanging node constraint equations.
410  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
411 
412  // The element matrix and right-hand-side are now built
413  // for this element. Add them to the global matrix and
414  // right-hand-side vector. The NumericMatrix::add_matrix()
415  // and NumericVector::add_vector() members do this for us.
416  matrix.add_matrix (Ke, dof_indices);
417  system.rhs->add_vector (Fe, dof_indices);
418  } // end of element loop
419 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
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:2330
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
unsigned int dim
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
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]...
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
Order default_quadrature_order() const
Definition: fe_type.h:371
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:75
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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.
Defines a dense submatrix for use in Finite Element-type computations.
Helper for building element sides that minimizes the construction of new elements.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2508
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
const DofMap & get_dof_map() const
Definition: system.h:2374
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
const Point & point(const unsigned int i) const
Definition: elem.h:2453
const SparseMatrix< Number > & get_system_matrix() const

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 70 of file systems_of_equations_ex1.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_stokes(), libMesh::System::attach_assemble_function(), libMesh::MeshTools::Generation::build_square(), libMesh::command_line_next(), libMesh::default_solver_package(), libMesh::FIRST, libMesh::EquationSystems::get_system(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::MeshTools::n_elem(), libMesh::EquationSystems::parameters, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::Real, libMesh::SECOND, libMesh::Parameters::set(), libMesh::TOLERANCE, libMesh::TRILINOS_SOLVERS, and libMesh::ExodusII_IO::write_equation_systems().

71 {
72  // Initialize libMesh.
73  LibMeshInit init (argc, argv);
74 
75  // This example requires a linear solver package.
76  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
77  "--enable-petsc, --enable-trilinos, or --enable-eigen");
78 
79  // Skip this 2D example if libMesh was compiled as 1D-only.
80  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
81 
82  // This example NaNs with the Trilinos solvers
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  // Get the mesh size from the command line.
90  const int n_elem =
91  libMesh::command_line_next("-n_elem", 15);
92 
93  // Use the MeshTools::Generation mesh generator to create a uniform
94  // 2D grid on the square [-1,1]^2. We instruct the mesh generator
95  // to build a mesh of 8x8 Quad9 elements. Building these
96  // higher-order elements allows us to use higher-order
97  // approximation, as in example 3.
99  n_elem, n_elem,
100  0., 1.,
101  0., 1.,
102  QUAD9);
103 
104  // Print information about the mesh to the screen.
105  mesh.print_info();
106 
107  // Create an equation systems object.
108  EquationSystems equation_systems (mesh);
109 
110  // Declare the system and its variables.
111  // Create a transient system named "Stokes"
112  LinearImplicitSystem & system =
113  equation_systems.add_system<LinearImplicitSystem> ("Stokes");
114 
115  // Add the variables "u" & "v" to "Stokes". They
116  // will be approximated using second-order approximation.
117  system.add_variable ("u", SECOND);
118  system.add_variable ("v", SECOND);
119 
120  // Add the variable "p" to "Stokes". This will
121  // be approximated with a first-order basis,
122  // providing an LBB-stable pressure-velocity pair.
123  system.add_variable ("p", FIRST);
124 
125  // Give the system a pointer to the matrix assembly
126  // function.
128 
129  // Initialize the data structures for the equation system.
130  equation_systems.init ();
131 
132  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
133  equation_systems.parameters.set<Real> ("linear solver tolerance") = TOLERANCE;
134 
135  // Prints information about the system to the screen.
136  equation_systems.print_info();
137 
138  // Assemble & solve the linear system,
139  // then write the solution.
140  equation_systems.get_system("Stokes").solve();
141 
142 #ifdef LIBMESH_HAVE_EXODUS_API
144  equation_systems);
145 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
146 
147  // All done.
148  return 0;
149 }
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1078
This is the EquationSystems class.
void assemble_stokes(EquationSystems &es, const std::string &system_name)
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
static constexpr Real TOLERANCE
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
MeshBase & mesh
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:90
SolverPackage default_solver_package()
Definition: libmesh.C:1117
virtual void write_equation_systems(const std::string &fname, const EquationSystems &es, const std::set< std::string > *system_names=nullptr) override
Writes out the solution for no specific time or timestep.
Definition: exodusII_io.C:2033
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1562
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1357
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:2161
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50