1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #include "libmesh/enum_solver_package.h"
76
77 // Bring in everything from the libMesh namespace
78 using namespace libMesh;
79
80 // Function prototype. This function will assemble the system matrix
81 // and right-hand-side.
82 void assemble (EquationSystems & es,
83 const std::string & system_name);
84
85 // Begin the main program. Note that the first
86 // statement in the program throws an error if
87 // you are in complex number mode, since this
88 // example is only intended to work with real
89 // numbers.
main(int argc,char ** argv)90 int main (int argc, char ** argv)
91 {
92 // Initialize libMesh.
93 LibMeshInit init (argc, argv);
94
95 #if !defined(LIBMESH_ENABLE_AMR)
96 libmesh_example_requires(false, "--enable-amr");
97 #else
98 libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
99
100 // Brief message to the user regarding the program name
101 // and command line arguments.
102
103 libMesh::out << "Running: " << argv[0];
104
105 for (int i=1; i<argc; i++)
106 libMesh::out << " " << argv[i];
107
108 libMesh::out << std::endl << std::endl;
109
110 // Skip this 2D example if libMesh was compiled as 1D-only.
111 libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
112
113 // Create a mesh, with dimension to be overridden later, distributed
114 // across the default MPI communicator.
115 Mesh mesh(init.comm());
116
117 // Create an equation systems object.
118 EquationSystems equation_systems (mesh);
119
120 MeshTools::Generation::build_square (mesh,
121 16,
122 16,
123 -1., 1.,
124 -1., 1.,
125 QUAD4);
126
127 LinearImplicitSystem & system =
128 equation_systems.add_system<LinearImplicitSystem>
129 ("System");
130
131 // Adds the variable "u" to "System". "u"
132 // will be approximated using first-order approximation.
133 system.add_variable ("u", FIRST);
134
135 // Also, we need to add two vectors. The tensor matrix v*w^T of
136 // these two vectors will be part of the system matrix.
137 system.add_vector("v", false);
138 system.add_vector("w", false);
139
140 // We need an additional matrix to be used for preconditioning since
141 // a shell matrix is not suitable for that.
142 system.add_matrix("Preconditioner");
143
144 // Give the system a pointer to the matrix assembly function.
145 system.attach_assemble_function (assemble);
146
147 // Initialize the data structures for the equation system.
148 equation_systems.init ();
149
150 // Prints information about the system to the screen.
151 equation_systems.print_info();
152
153 equation_systems.parameters.set<unsigned int>
154 ("linear solver maximum iterations") = 250;
155 equation_systems.parameters.set<Real>
156 ("linear solver tolerance") = TOLERANCE;
157
158 // Refine arbitrarily some elements.
159 for (unsigned int i=0; i<2; i++)
160 {
161 MeshRefinement mesh_refinement(mesh);
162 for (auto & elem : mesh.element_ptr_range())
163 {
164 if (elem->active())
165 {
166 if ((elem->id()%20)>8)
167 elem->set_refinement_flag(Elem::REFINE);
168 else
169 elem->set_refinement_flag(Elem::DO_NOTHING);
170 }
171 else
172 elem->set_refinement_flag(Elem::INACTIVE);
173 }
174 mesh_refinement.refine_elements();
175 equation_systems.reinit();
176 }
177
178 // Prints information about the system to the screen.
179 equation_systems.print_info();
180
181 // Before assembling the matrix, we have to clear the two
182 // vectors that form the tensor matrix (since this is not performed
183 // automatically).
184 system.get_vector("v").init(system.n_dofs(), system.n_local_dofs());
185 system.get_vector("w").init(system.n_dofs(), system.n_local_dofs());
186
187 // We need a shell matrix to solve. There is currently no way to
188 // store the shell matrix in the system. We just create it locally
189 // here (a shell matrix does not occupy much memory).
190 SumShellMatrix<Number> shellMatrix(system.comm());
191 TensorShellMatrix<Number> shellMatrix0(system.get_vector("v"), system.get_vector("w"));
192 shellMatrix.matrices.push_back(&shellMatrix0);
193 SparseShellMatrix<Number> shellMatrix1(*system.matrix);
194 shellMatrix.matrices.push_back(&shellMatrix1);
195
196 // Attach that to the system.
197 system.attach_shell_matrix(&shellMatrix);
198
199 // Reset the preconditioning matrix to zero (for the system matrix,
200 // the same thing is done automatically).
201 system.get_matrix("Preconditioner").zero();
202
203 // Assemble & solve the linear system
204 system.solve();
205
206 // Detach the shell matrix from the system since it will go out of
207 // scope. Nobody should solve the system outside this function.
208 system.detach_shell_matrix();
209
210 // Print a nice message.
211 libMesh::out << "Solved linear system in "
212 << system.n_linear_iterations()
213 << " iterations, residual norm is "
214 << system.final_linear_residual()
215 << "."
216 << std::endl;
217
218 #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
219 // Write result to file.
220 VTKIO(mesh).write_equation_systems ("out.pvtu", equation_systems);
221 #endif // #ifdef LIBMESH_HAVE_VTK
222
223 #endif // #ifndef LIBMESH_ENABLE_AMR
224
225 return 0;
226 }
227
228
229
230 // This function defines the assembly routine. It is responsible for
231 // computing the proper matrix entries for the element stiffness
232 // matrices and right-hand sides.
assemble(EquationSystems & es,const std::string & system_name)233 void assemble (EquationSystems & es,
234 const std::string & system_name)
235 {
236 // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
237 libmesh_ignore(es, 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 std::unique_ptr<FEBase>. This can be thought
261 // of as a pointer that will clean up after itself.
262 std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
263 std::unique_ptr<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".
301 DenseMatrix<Number> Ke;
302 DenseVector<Number> Fe;
303
304 // Analogous data structures for thw two vectors v and w that form
305 // the tensor shell matrix.
306 DenseVector<Number> Ve;
307 DenseVector<Number> We;
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 const unsigned int n_dofs =
340 cast_int<unsigned int>(dof_indices.size());
341
342 Ke.resize (n_dofs, n_dofs);
343
344 Fe.resize (n_dofs);
345 Ve.resize (n_dofs);
346 We.resize (n_dofs);
347
348 // Now we will build the element matrix and right-hand-side.
349 // Constructing the RHS requires the solution and its
350 // gradient from the previous timestep. This myst be
351 // calculated at each quadrature point by summing the
352 // solution degree-of-freedom values by the appropriate
353 // weight functions.
354 for (unsigned int qp=0; qp<qrule.n_points(); qp++)
355 {
356 // Now compute the element matrix and RHS contributions.
357 for (unsigned int i=0; i<n_dofs; i++)
358 {
359 // The RHS contribution
360 Fe(i) += JxW[qp]*phi[i][qp];
361
362 for (unsigned int j=0; j<n_dofs; j++)
363 {
364 // The matrix contribution
365 Ke(i,j) += JxW[qp]*(
366 // Stiffness matrix
367 (dphi[i][qp]*dphi[j][qp])
368 );
369 }
370
371 // V and W are the same for this example.
372 Ve(i) += JxW[qp]*phi[i][qp];
373 We(i) += JxW[qp]*phi[i][qp];
374 }
375 }
376
377 // At this point the interior element integration has
378 // been completed. However, we have not yet addressed
379 // boundary conditions. For this example we will only
380 // consider simple Dirichlet boundary conditions imposed
381 // via the penalty method.
382 //
383 // The following loops 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 {
387 // The penalty value.
388 const Real penalty = 1.e10;
389
390 // The following loops over the sides of the element.
391 // If the element has no neighbor on a side then that
392 // side MUST live on a boundary of the domain.
393 for (auto s : elem->side_index_range())
394 if (elem->neighbor_ptr(s) == nullptr)
395 {
396 fe_face->reinit(elem, s);
397
398 for (unsigned int qp=0; qp<qface.n_points(); qp++)
399 {
400 // Matrix contribution
401 for (unsigned int i=0; i<n_dofs; i++)
402 for (unsigned int j=0; j<n_dofs; j++)
403 Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
404 }
405 }
406 }
407
408
409 // We have now built the element matrix and RHS vector in terms
410 // of the element degrees of freedom. However, it is possible
411 // that some of the element DOFs are constrained to enforce
412 // solution continuity, i.e. they are not really "free". We need
413 // to constrain those DOFs in terms of non-constrained DOFs to
414 // ensure a continuous solution. The
415 // DofMap::constrain_element_matrix_and_vector() method does
416 // just that.
417
418 // However, constraining both the sparse matrix (and right hand
419 // side) plus the rank 1 matrix is tricky. The dof_indices
420 // vector has to be backed up for that because the constraining
421 // functions modify it.
422
423 std::vector<dof_id_type> dof_indices_backup(dof_indices);
424 dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
425 dof_indices = dof_indices_backup;
426 dof_map.constrain_element_dyad_matrix(Ve, We, dof_indices);
427
428 // The element matrix and right-hand-side are now built
429 // for this element. Add them to the global matrix and
430 // right-hand-side vector. The SparseMatrix::add_matrix()
431 // and NumericVector::add_vector() members do this for us.
432 system.matrix->add_matrix (Ke, dof_indices);
433 system.get_matrix("Preconditioner").add_matrix (Ke, dof_indices);
434 system.rhs->add_vector (Fe, dof_indices);
435 system.get_vector("v").add_vector(Ve, dof_indices);
436 system.get_vector("w").add_vector(We, dof_indices);
437 }
438 // Finished computing the system matrix and right-hand side.
439
440 // Matrices and vectors must be closed manually. This is necessary
441 // because the matrix is not directly used as the system matrix (in
442 // which case the solver closes it) but as a part of a shell matrix.
443 system.matrix->close();
444 system.get_matrix("Preconditioner").close();
445 system.rhs->close();
446 system.get_vector("v").close();
447 system.get_vector("w").close();
448
449 #endif // #ifdef LIBMESH_ENABLE_AMR
450 }
451