libMesh
miscellaneous_ex4.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>Miscellaneous Example 4 - Using a shell matrix</h1>
21 // \author Tim Kroger
22 // \date 2008
23 //
24 // This example solves the equation
25 //
26 // \f$-\Delta u+\int u = 1\f$
27 //
28 // with homogeneous Dirichlet boundary conditions. This system has
29 // a full system matrix which can be written as the sum of of sparse
30 // matrix and a rank 1 matrix. The shell matrix concept is used to
31 // solve this problem.
32 //
33 // The problem is solved in parallel on a non-uniform grid in order
34 // to demonstrate all the techniques that are required for this.
35 // The grid is fixed, however, i.e. no adaptive mesh refinement is
36 // used, so that the example remains simple.
37 //
38 // The example is 2d; extension to 3d is straight forward.
39 
40 // C++ include files that we need
41 #include <iostream>
42 #include <algorithm>
43 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
44 #include <cmath>
45 
46 // Basic include file needed for the mesh functionality.
47 #include "libmesh/libmesh.h"
48 #include "libmesh/mesh.h"
49 #include "libmesh/mesh_refinement.h"
50 #include "libmesh/vtk_io.h"
51 #include "libmesh/equation_systems.h"
52 #include "libmesh/fe.h"
53 #include "libmesh/quadrature_gauss.h"
54 #include "libmesh/dof_map.h"
55 #include "libmesh/sparse_matrix.h"
56 #include "libmesh/numeric_vector.h"
57 #include "libmesh/dense_matrix.h"
58 #include "libmesh/dense_vector.h"
59 #include "libmesh/mesh_generation.h"
60 #include "libmesh/sum_shell_matrix.h"
61 #include "libmesh/tensor_shell_matrix.h"
62 #include "libmesh/sparse_shell_matrix.h"
63 #include "libmesh/mesh_refinement.h"
64 
65 #include "libmesh/getpot.h"
66 
67 // This example will solve a linear transient system,
68 // so we need to include the TransientLinearImplicitSystem definition.
69 #include "libmesh/transient_system.h"
70 #include "libmesh/linear_implicit_system.h"
71 #include "libmesh/vector_value.h"
72 
73 // The definition of a geometric element
74 #include "libmesh/elem.h"
75 
76 // Bring in everything from the libMesh namespace
77 using namespace libMesh;
78 
79 // Function prototype. This function will assemble the system matrix
80 // and right-hand-side.
81 void assemble (EquationSystems & es,
82  const std::string & system_name);
83 
84 // Begin the main program. Note that the first
85 // statement in the program throws an error if
86 // you are in complex number mode, since this
87 // example is only intended to work with real
88 // numbers.
89 int main (int argc, char ** argv)
90 {
91  // Initialize libMesh.
92  LibMeshInit init (argc, argv);
93 
94 #if !defined(LIBMESH_ENABLE_AMR)
95  libmesh_example_requires(false, "--enable-amr");
96 #else
97  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
98 
99  // Brief message to the user regarding the program name
100  // and command line arguments.
101 
102  libMesh::out << "Running: " << argv[0];
103 
104  for (int i=1; i<argc; i++)
105  libMesh::out << " " << argv[i];
106 
107  libMesh::out << std::endl << std::endl;
108 
109  // Skip this 2D example if libMesh was compiled as 1D-only.
110  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
111 
112  // Create a mesh, with dimension to be overridden later, distributed
113  // across the default MPI communicator.
114  Mesh mesh(init.comm());
115 
116  // Create an equation systems object.
117  EquationSystems equation_systems (mesh);
118 
120  16,
121  16,
122  -1., 1.,
123  -1., 1.,
124  QUAD4);
125 
126  LinearImplicitSystem & system =
127  equation_systems.add_system<LinearImplicitSystem>
128  ("System");
129 
130  // Adds the variable "u" to "System". "u"
131  // will be approximated using first-order approximation.
132  system.add_variable ("u", FIRST);
133 
134  // Also, we need to add two vectors. The tensor matrix v*w^T of
135  // these two vectors will be part of the system matrix.
136  system.add_vector("v", false);
137  system.add_vector("w", false);
138 
139  // We need an additional matrix to be used for preconditioning since
140  // a shell matrix is not suitable for that.
141  system.add_matrix("Preconditioner");
142 
143  // Give the system a pointer to the matrix assembly function.
145 
146  // Initialize the data structures for the equation system.
147  equation_systems.init ();
148 
149  // Prints information about the system to the screen.
150  equation_systems.print_info();
151 
152  equation_systems.parameters.set<unsigned int>
153  ("linear solver maximum iterations") = 250;
154  equation_systems.parameters.set<Real>
155  ("linear solver tolerance") = TOLERANCE;
156 
157  // Refine arbitrarily some elements.
158  for (unsigned int i=0; i<2; i++)
159  {
160  MeshRefinement mesh_refinement(mesh);
161  for (auto & elem : mesh.element_ptr_range())
162  {
163  if (elem->active())
164  {
165  if ((elem->id()%20)>8)
166  elem->set_refinement_flag(Elem::REFINE);
167  else
168  elem->set_refinement_flag(Elem::DO_NOTHING);
169  }
170  else
171  elem->set_refinement_flag(Elem::INACTIVE);
172  }
173  mesh_refinement.refine_elements();
174  equation_systems.reinit();
175  }
176 
177  // Prints information about the system to the screen.
178  equation_systems.print_info();
179 
180  // Before assembling the matrix, we have to clear the two
181  // vectors that form the tensor matrix (since this is not performed
182  // automatically).
183  system.get_vector("v").init(system.n_dofs(), system.n_local_dofs());
184  system.get_vector("w").init(system.n_dofs(), system.n_local_dofs());
185 
186  // We need a shell matrix to solve. There is currently no way to
187  // store the shell matrix in the system. We just create it locally
188  // here (a shell matrix does not occupy much memory).
189  SumShellMatrix<Number> shellMatrix(system.comm());
190  TensorShellMatrix<Number> shellMatrix0(system.get_vector("v"), system.get_vector("w"));
191  shellMatrix.matrices.push_back(&shellMatrix0);
192  SparseShellMatrix<Number> shellMatrix1(*system.matrix);
193  shellMatrix.matrices.push_back(&shellMatrix1);
194 
195  // Attach that to the system.
196  system.attach_shell_matrix(&shellMatrix);
197 
198  // Reset the preconditioning matrix to zero (for the system matrix,
199  // the same thing is done automatically).
200  system.get_matrix("Preconditioner").zero();
201 
202  // Assemble & solve the linear system
203  system.solve();
204 
205  // Detach the shell matrix from the system since it will go out of
206  // scope. Nobody should solve the system outside this function.
207  system.detach_shell_matrix();
208 
209  // Print a nice message.
210  libMesh::out << "Solved linear system in "
211  << system.n_linear_iterations()
212  << " iterations, residual norm is "
213  << system.final_linear_residual()
214  << "."
215  << std::endl;
216 
217 #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
218  // Write result to file.
219  VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
220 #endif // #ifdef LIBMESH_HAVE_VTK
221 
222 #endif // #ifndef LIBMESH_ENABLE_AMR
223 
224  return 0;
225 }
226 
227 
228 
229 // This function defines the assembly routine. It is responsible for
230 // computing the proper matrix entries for the element stiffness
231 // matrices and right-hand sides.
233  const std::string & system_name)
234 {
235  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
236  libmesh_ignore(es);
237  libmesh_ignore(system_name);
238 
239 #ifdef LIBMESH_ENABLE_AMR
240  // It is a good idea to make sure we are assembling
241  // the proper system.
242  libmesh_assert_equal_to (system_name, "System");
243 
244  // Get a constant reference to the mesh object.
245  const MeshBase & mesh = es.get_mesh();
246 
247  // The dimension that we are running
248  const unsigned int dim = mesh.mesh_dimension();
249 
250  // Get a reference to the Convection-Diffusion system object.
251  LinearImplicitSystem & system =
252  es.get_system<LinearImplicitSystem> ("System");
253 
254  // Get the Finite Element type for the first (and only)
255  // variable in the system.
256  FEType fe_type = system.variable_type(0);
257 
258  // Build a Finite Element object of the specified type. Since the
259  // FEBase::build() member dynamically creates memory we will
260  // store the object as a UniquePtr<FEBase>. This can be thought
261  // of as a pointer that will clean up after itself.
262  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
263  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
264 
265  // A Gauss quadrature rule for numerical integration.
266  // Let the FEType object decide what order rule is appropriate.
267  QGauss qrule (dim, fe_type.default_quadrature_order());
268  QGauss qface (dim-1, fe_type.default_quadrature_order());
269 
270  // Tell the finite element object to use our quadrature rule.
271  fe->attach_quadrature_rule (&qrule);
272  fe_face->attach_quadrature_rule (&qface);
273 
274  // Here we define some references to cell-specific data that
275  // will be used to assemble the linear system. We will start
276  // with the element Jacobian * quadrature weight at each integration point.
277  const std::vector<Real> & JxW = fe->get_JxW();
278  const std::vector<Real> & JxW_face = fe_face->get_JxW();
279 
280  // The element shape functions evaluated at the quadrature points.
281  const std::vector<std::vector<Real>> & phi = fe->get_phi();
282  const std::vector<std::vector<Real>> & psi = fe_face->get_phi();
283 
284  // The element shape function gradients evaluated at the quadrature
285  // points.
286  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
287 
288  // The XY locations of the quadrature points used for face integration
289  //const std::vector<Point>& qface_points = fe_face->get_xyz();
290 
291  // A reference to the DofMap object for this system. The DofMap
292  // object handles the index translation from node and element numbers
293  // to degree of freedom numbers. We will talk more about the DofMap
294  // in future examples.
295  const DofMap & dof_map = system.get_dof_map();
296 
297  // Define data structures to contain the element matrix
298  // and right-hand-side vector contribution. Following
299  // basic finite element terminology we will denote these
300  // "Ke" and "Fe".
303 
304  // Analogous data structures for thw two vectors v and w that form
305  // the tensor shell matrix.
308 
309  // This vector will hold the degree of freedom indices for
310  // the element. These define where in the global system
311  // the element degrees of freedom get mapped.
312  std::vector<dof_id_type> dof_indices;
313 
314  // Now we will loop over all the elements in the mesh that
315  // live on the local processor. We will compute the element
316  // matrix and right-hand-side contribution. Since the mesh
317  // will be refined we want to only consider the ACTIVE elements,
318  // hence we use a variant of the active_elem_iterator.
319  for (const auto & elem : mesh.active_local_element_ptr_range())
320  {
321  // Get the degree of freedom indices for the
322  // current element. These define where in the global
323  // matrix and right-hand-side this element will
324  // contribute to.
325  dof_map.dof_indices (elem, dof_indices);
326 
327  // Compute the element-specific data for the current
328  // element. This involves computing the location of the
329  // quadrature points (q_point) and the shape functions
330  // (phi, dphi) for the current element.
331  fe->reinit (elem);
332 
333  // Zero the element matrix and right-hand side before
334  // summing them. We use the resize member here because
335  // the number of degrees of freedom might have changed from
336  // the last element. Note that this will be the case if the
337  // element type is different (i.e. the last element was a
338  // triangle, now we are on a quadrilateral).
339  Ke.resize (dof_indices.size(),
340  dof_indices.size());
341 
342  Fe.resize (dof_indices.size());
343  Ve.resize (dof_indices.size());
344  We.resize (dof_indices.size());
345 
346  // Now we will build the element matrix and right-hand-side.
347  // Constructing the RHS requires the solution and its
348  // gradient from the previous timestep. This myst be
349  // calculated at each quadrature point by summing the
350  // solution degree-of-freedom values by the appropriate
351  // weight functions.
352  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
353  {
354  // Now compute the element matrix and RHS contributions.
355  for (std::size_t i=0; i<phi.size(); i++)
356  {
357  // The RHS contribution
358  Fe(i) += JxW[qp]*phi[i][qp];
359 
360  for (std::size_t j=0; j<phi.size(); j++)
361  {
362  // The matrix contribution
363  Ke(i,j) += JxW[qp]*(
364  // Stiffness matrix
365  (dphi[i][qp]*dphi[j][qp])
366  );
367  }
368 
369  // V and W are the same for this example.
370  Ve(i) += JxW[qp]*phi[i][qp];
371  We(i) += JxW[qp]*phi[i][qp];
372  }
373  }
374 
375  // At this point the interior element integration has
376  // been completed. However, we have not yet addressed
377  // boundary conditions. For this example we will only
378  // consider simple Dirichlet boundary conditions imposed
379  // via the penalty method.
380  //
381  // The following loops over the sides of the element.
382  // If the element has no neighbor on a side then that
383  // side MUST live on a boundary of the domain.
384  {
385  // The penalty value.
386  const Real penalty = 1.e10;
387 
388  // The following loops over the sides of the element.
389  // If the element has no neighbor on a side then that
390  // side MUST live on a boundary of the domain.
391  for (auto s : elem->side_index_range())
392  if (elem->neighbor_ptr(s) == libmesh_nullptr)
393  {
394  fe_face->reinit(elem, s);
395 
396  for (unsigned int qp=0; qp<qface.n_points(); qp++)
397  {
398  // Matrix contribution
399  for (std::size_t i=0; i<psi.size(); i++)
400  for (std::size_t j=0; j<psi.size(); j++)
401  Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
402  }
403  }
404  }
405 
406 
407  // We have now built the element matrix and RHS vector in terms
408  // of the element degrees of freedom. However, it is possible
409  // that some of the element DOFs are constrained to enforce
410  // solution continuity, i.e. they are not really "free". We need
411  // to constrain those DOFs in terms of non-constrained DOFs to
412  // ensure a continuous solution. The
413  // DofMap::constrain_element_matrix_and_vector() method does
414  // just that.
415 
416  // However, constraining both the sparse matrix (and right hand
417  // side) plus the rank 1 matrix is tricky. The dof_indices
418  // vector has to be backed up for that because the constraining
419  // functions modify it.
420 
421  std::vector<dof_id_type> dof_indices_backup(dof_indices);
422  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
423  dof_indices = dof_indices_backup;
424  dof_map.constrain_element_dyad_matrix(Ve, We, dof_indices);
425 
426  // The element matrix and right-hand-side are now built
427  // for this element. Add them to the global matrix and
428  // right-hand-side vector. The SparseMatrix::add_matrix()
429  // and NumericVector::add_vector() members do this for us.
430  system.matrix->add_matrix (Ke, dof_indices);
431  system.get_matrix("Preconditioner").add_matrix (Ke, dof_indices);
432  system.rhs->add_vector (Fe, dof_indices);
433  system.get_vector("v").add_vector(Ve, dof_indices);
434  system.get_vector("w").add_vector(We, dof_indices);
435  }
436  // Finished computing the system matrix and right-hand side.
437 
438  // Matrices and vectors must be closed manually. This is necessary
439  // because the matrix is not directly used as the system matrix (in
440  // which case the solver closes it) but as a part of a shell matrix.
441  system.matrix->close();
442  system.get_matrix("Preconditioner").close();
443  system.rhs->close();
444  system.get_vector("v").close();
445  system.get_vector("w").close();
446 
447 #endif // #ifdef LIBMESH_ENABLE_AMR
448 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
unsigned int n_linear_iterations() const
void attach_shell_matrix(ShellMatrix< Number > *shell_matrix)
This function enables the user to provide a shell matrix, i.e.
This is the EquationSystems class.
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
unsigned int dim
bool refine_elements()
Only refines the user-requested elements.
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 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
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
This class combines any number of shell matrices to a single shell matrix by summing them together...
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.
This class implements reading and writing meshes in the VTK format.
Definition: vtk_io.h:59
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
void assemble(EquationSystems &es, const std::string &system_name)
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
Shell matrix that is given by a tensor product of two vectors, i.e.
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.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
virtual void zero()=0
Set all entries to 0.
virtual SimpleRange< element_iterator > element_ptr_range()=0
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
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
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
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
dof_id_type n_local_dofs() const
Definition: system.C:185
void libmesh_ignore(const T &)
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Parallel::Communicator & comm() const
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
This class allows to use any SparseMatrix object as a shell matrix.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
int main(int argc, char **argv)
This class implements specific orders of Gauss quadrature.
const MeshBase & get_mesh() const
void detach_shell_matrix()
Detaches a shell matrix.
dof_id_type n_dofs() const
Definition: system.C:148
const T_sys & get_system(const std::string &name) const
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void constrain_element_dyad_matrix(DenseVector< Number > &v, DenseVector< Number > &w, std::vector< dof_id_type > &row_dofs, bool asymmetric_constraint_rows=true) const
Constrains a dyadic element matrix B = v w&#39;.
Definition: dof_map.h:1811
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
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.