libMesh
eigenproblems_ex3.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>Eigenproblems Example 3 - Can you "hear the shape" of a drum?</h1>
21 // \author David Knezevic
22 // \date 2012
23 //
24 // The sound that a drum makes is determined by it's resonant frequencies,
25 // which are given by the eigenvalues of the Laplacian. "Can One Hear the
26 // Shape of a Drum?" was the title of an article by Mark Kac
27 // in the American Mathematical Monthly in 1966, where he raised the question:
28 // If we know all the eigenvalues of a drum, can we uniquely determine it's shape?
29 // This question was resolved in 1992, when Gordon, Webb, and Wolpert constructed
30 // a pair of regions in 2D that have different shapes but identical eigenvalues.
31 // So the answer to Kac's question is no: the spectrum of the Laplacian does
32 // not uniquely determine the shape of the domain.
33 
34 // In this example, we compute the first few eigenvalues of the two domains proposed
35 // by Gordon, Webb and Wolpert. This amounts to solving a generalized eigenvalue
36 // problem in each case. The computed eigenvalues are stored in drum1_evals.txt and
37 // drum2_evals.txt. We can compare these to the (highly accurate) values reported in:
38 // T.A. Driscoll, "Eigenmodes of Isospectral Drums", SIAM Review, Vol. 39, No. 1, pp. 1-17, 1997.
39 //
40 // The first five eigenvalues from Driscoll are listed below (the author states that
41 // "all digits shown are believed to be correct"):
42 // 2.53794399980
43 // 3.65550971352
44 // 5.17555935622
45 // 6.53755744376
46 // 7.24807786256
47 
48 
49 // C++ include files
50 #include <fstream>
51 
52 // libMesh include files.
53 #include "libmesh/libmesh.h"
54 #include "libmesh/mesh.h"
55 #include "libmesh/mesh_generation.h"
56 #include "libmesh/exodusII_io.h"
57 #include "libmesh/condensed_eigen_system.h"
58 #include "libmesh/equation_systems.h"
59 #include "libmesh/fe.h"
60 #include "libmesh/quadrature_gauss.h"
61 #include "libmesh/dense_matrix.h"
62 #include "libmesh/sparse_matrix.h"
63 #include "libmesh/numeric_vector.h"
64 #include "libmesh/dof_map.h"
65 #include "libmesh/fe_interface.h"
66 #include "libmesh/getpot.h"
67 #include "libmesh/elem.h"
68 #include "libmesh/zero_function.h"
69 #include "libmesh/dirichlet_boundaries.h"
70 
71 #define BOUNDARY_ID 100
72 
73 // Bring in everything from the libMesh namespace
74 using namespace libMesh;
75 
76 
77 // Function prototype. This is the function that will assemble
78 // the eigen system. Here, we will simply assemble a mass matrix.
80  const std::string & system_name);
81 
82 // We store the Dirichlet dofs in a set in order to impose the boundary conditions
84  const std::string & system_name,
85  std::set<unsigned int> & global_dirichlet_dofs_set);
86 
87 
88 int main (int argc, char ** argv)
89 {
90  // Initialize libMesh and the dependent libraries.
91  LibMeshInit init (argc, argv);
92 
93  // This example uses an ExodusII input file
94 #ifndef LIBMESH_HAVE_EXODUS_API
95  libmesh_example_requires(false, "--enable-exodus");
96 #endif
97 
98  // This example is designed for the SLEPc eigen solver interface.
99 #ifndef LIBMESH_HAVE_SLEPC
100  if (init.comm().rank() == 0)
101  libMesh::err << "ERROR: This example requires libMesh to be\n"
102  << "compiled with SLEPc eigen solvers support!"
103  << std::endl;
104 
105  return 0;
106 #else
107 
108 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
109  // SLEPc currently gives us a nasty crash with Real==float
110  libmesh_example_requires(false, "--disable-singleprecision");
111 #endif
112 
113 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
114  // SLEPc currently gives us an "inner product not well defined" with
115  // Number==complex
116  libmesh_example_requires(false, "--disable-complex");
117 #endif
118 
119  // Tell the user what we are doing.
120  {
121  libMesh::out << "Running " << argv[0];
122 
123  for (int i=1; i<argc; i++)
124  libMesh::out << " " << argv[i];
125 
126  libMesh::out << std::endl << std::endl;
127  }
128 
129  // Skip this 2D example if libMesh was compiled as 1D-only.
130  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
131 
132  // Use GetPot to parse the command line arguments
133  GetPot command_line (argc, argv);
134 
135  // Read the mesh name from the command line
136  std::string mesh_name = "";
137  if (command_line.search(1, "-mesh_name"))
138  mesh_name = command_line.next(mesh_name);
139 
140  // Also, read in the index of the eigenvector that we should plot
141  // (zero-based indexing, as usual!)
142  unsigned int plotting_index = 0;
143  if (command_line.search(1, "-plotting_index"))
144  plotting_index = command_line.next(plotting_index);
145 
146  // Finally, read in the number of eigenpairs we want to compute!
147  unsigned int n_evals = 0;
148  if (command_line.search(1, "-n_evals"))
149  n_evals = command_line.next(n_evals);
150 
151  // Append the .e to mesh_name
152  std::ostringstream mesh_name_exodus;
153  mesh_name_exodus << mesh_name << "_mesh.e";
154 
155  // Create a mesh, with dimension to be overridden by the file, on
156  // the default MPI communicator.
157  Mesh mesh(init.comm());
158 
159  mesh.read(mesh_name_exodus.str());
160 
161  // Add boundary IDs to this mesh so that we can use DirichletBoundary
162  // Each processor should know about each boundary condition it can
163  // see, so we loop over all elements, not just local elements.
164  for (const auto & elem : mesh.element_ptr_range())
165  for (auto side : elem->side_index_range())
166  if (elem->neighbor_ptr (side) == NULL)
167  mesh.get_boundary_info().add_side(elem, side, BOUNDARY_ID);
168 
169  // Print information about the mesh to the screen.
170  mesh.print_info();
171 
172  // Create an equation systems object.
173  EquationSystems equation_systems (mesh);
174 
175  // Create a CondensedEigenSystem named "Eigensystem" and (for convenience)
176  // use a reference to the system we create.
177  CondensedEigenSystem & eigen_system =
178  equation_systems.add_system<CondensedEigenSystem> ("Eigensystem");
179 
180  // Declare the system variables.
181  // Adds the variable "p" to "Eigensystem". "p"
182  // will be approximated using second-order approximation.
183  eigen_system.add_variable("p", SECOND);
184 
185  // Give the system a pointer to the matrix assembly
186  // function defined below.
188 
189  // Set the number of requested eigenpairs n_evals and the number
190  // of basis vectors used in the solution algorithm.
191  equation_systems.parameters.set<unsigned int>("eigenpairs") = n_evals;
192  equation_systems.parameters.set<unsigned int>("basis vectors") = n_evals*3;
193 
194  // Set the solver tolerance and the maximum number of iterations.
195  equation_systems.parameters.set<Real>("linear solver tolerance") = pow(TOLERANCE, 5./3.);
196  equation_systems.parameters.set<unsigned int>
197  ("linear solver maximum iterations") = 1000;
198 
199  // Set the type of the problem, here we deal with
200  // a generalized Hermitian problem.
201  eigen_system.set_eigenproblem_type(GHEP);
202 
203  // Order the eigenvalues "smallest first"
204  eigen_system.eigen_solver->set_position_of_spectrum(SMALLEST_MAGNITUDE);
205 
206  {
207  std::set<boundary_id_type> boundary_ids;
208  boundary_ids.insert(BOUNDARY_ID);
209 
210  std::vector<unsigned int> variables;
211  variables.push_back(0);
212 
213  ZeroFunction<> zf;
214 
215  // Most DirichletBoundary users will want to supply a "locally
216  // indexed" functor
217  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
219 
220  eigen_system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
221  }
222 
223  // Initialize the data structures for the equation system.
224  equation_systems.init();
225 
226  // Prints information about the system to the screen.
227  equation_systems.print_info();
228 
229  eigen_system.initialize_condensed_dofs();
230 
231  // Solve the system "Eigensystem".
232  eigen_system.solve();
233 
234  // Get the number of converged eigen pairs.
235  unsigned int nconv = eigen_system.get_n_converged();
236 
237  libMesh::out << "Number of converged eigenpairs: "
238  << nconv
239  << "\n"
240  << std::endl;
241 
242  if (plotting_index > n_evals)
243  {
244  libMesh::out << "WARNING: Solver did not converge for the requested eigenvector!" << std::endl;
245  }
246 
247  // write out all of the computed eigenvalues and plot the specified eigenvector
248  std::ostringstream eigenvalue_output_name;
249  eigenvalue_output_name << mesh_name << "_evals.txt";
250  std::ofstream evals_file(eigenvalue_output_name.str().c_str());
251 
252  for (unsigned int i=0; i<nconv; i++)
253  {
254  std::pair<Real,Real> eval = eigen_system.get_eigenpair(i);
255 
256  // The eigenvalues should be real!
257  libmesh_assert_less (eval.second, TOLERANCE);
258  evals_file << eval.first << std::endl;
259 
260  // plot the specified eigenvector
261  if (i == plotting_index)
262  {
263 #ifdef LIBMESH_HAVE_EXODUS_API
264  // Write the eigen vector to file.
265  std::ostringstream eigenvector_output_name;
266  eigenvector_output_name << mesh_name << "_evec.e";
267  ExodusII_IO (mesh).write_equation_systems (eigenvector_output_name.str(), equation_systems);
268 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
269  }
270  }
271 
272  evals_file.close();
273 
274 #endif // LIBMESH_HAVE_SLEPC
275 
276  // All done.
277  return 0;
278 }
279 
280 
281 
283  const std::string & libmesh_dbg_var(system_name))
284 {
285 
286  // It is a good idea to make sure we are assembling
287  // the proper system.
288  libmesh_assert_equal_to (system_name, "Eigensystem");
289 
290 #ifdef LIBMESH_HAVE_SLEPC
291 
292  // Get a constant reference to the mesh object.
293  const MeshBase & mesh = es.get_mesh();
294 
295  // The dimension that we are running.
296  const unsigned int dim = mesh.mesh_dimension();
297 
298  // Get a reference to our system.
299  EigenSystem & eigen_system = es.get_system<EigenSystem> ("Eigensystem");
300 
301  // Get a constant reference to the Finite Element type
302  // for the first (and only) variable in the system.
303  FEType fe_type = eigen_system.get_dof_map().variable_type(0);
304 
305  // A reference to the two system matrices
306  SparseMatrix<Number> & matrix_A = *eigen_system.matrix_A;
307  SparseMatrix<Number> & matrix_B = *eigen_system.matrix_B;
308 
309  // Build a Finite Element object of the specified type. Since the
310  // FEBase::build() member dynamically creates memory we will
311  // store the object as a UniquePtr<FEBase>. This can be thought
312  // of as a pointer that will clean up after itself.
313  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
314 
315  // A Gauss quadrature rule for numerical integration.
316  // Use the default quadrature order.
317  QGauss qrule (dim, fe_type.default_quadrature_order());
318 
319  // Tell the finite element object to use our quadrature rule.
320  fe->attach_quadrature_rule (&qrule);
321 
322  // The element Jacobian * quadrature weight at each integration point.
323  const std::vector<Real> & JxW = fe->get_JxW();
324 
325  // The element shape functions evaluated at the quadrature points.
326  const std::vector<std::vector<Real>> & phi = fe->get_phi();
327 
328  // The element shape function gradients evaluated at the quadrature
329  // points.
330  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
331 
332  // A reference to the DofMap object for this system. The DofMap
333  // object handles the index translation from node and element numbers
334  // to degree of freedom numbers.
335  const DofMap & dof_map = eigen_system.get_dof_map();
336 
337  // The element mass and stiffness matrices.
340 
341  // This vector will hold the degree of freedom indices for
342  // the element. These define where in the global system
343  // the element degrees of freedom get mapped.
344  std::vector<dof_id_type> dof_indices;
345 
346 
347  // Now we will loop over all the elements in the mesh that
348  // live on the local processor. We will compute the element
349  // matrix and right-hand-side contribution. In case users
350  // later modify this program to include refinement, we will
351  // be safe and will only consider the active elements;
352  // hence we use a variant of the active_elem_iterator.
353  for (const auto & elem : mesh.active_local_element_ptr_range())
354  {
355  // Get the degree of freedom indices for the
356  // current element. These define where in the global
357  // matrix and right-hand-side this element will
358  // contribute to.
359  dof_map.dof_indices (elem, dof_indices);
360 
361  // Compute the element-specific data for the current
362  // element. This involves computing the location of the
363  // quadrature points (q_point) and the shape functions
364  // (phi, dphi) for the current element.
365  fe->reinit (elem);
366 
367  // Zero the element matrices before
368  // summing them. We use the resize member here because
369  // the number of degrees of freedom might have changed from
370  // the last element. Note that this will be the case if the
371  // element type is different (i.e. the last element was a
372  // triangle, now we are on a quadrilateral).
373  Ke.resize (dof_indices.size(), dof_indices.size());
374  Me.resize (dof_indices.size(), dof_indices.size());
375 
376  // Now loop over the quadrature points. This handles
377  // the numeric integration.
378  //
379  // We will build the element matrix. This involves
380  // a double loop to integrate the test functions (i) against
381  // the trial functions (j).
382  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
383  for (std::size_t i=0; i<phi.size(); i++)
384  for (std::size_t j=0; j<phi.size(); j++)
385  {
386  Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp];
387  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
388  }
389 
390  // The calls to constrain_element_matrix below have no effect in
391  // the current example. However, if users modify this example to
392  // include hanging nodes due to mesh refinement, for example,
393  // then it is essential to call constrain_element_matrix.
394  // As a result we include constrain_element_matrix here to
395  // ensure this example is ready to be used with hanging nodes.
396  // (Note that constrained rows/cols will be eliminated from
397  // the eigenproblem by the CondensedEigenSystem.)
398  dof_map.constrain_element_matrix(Ke, dof_indices, false);
399  dof_map.constrain_element_matrix(Me, dof_indices, false);
400 
401  // Finally, simply add the element contribution to the
402  // overall matrices A and B.
403  matrix_A.add_matrix (Ke, dof_indices);
404  matrix_B.add_matrix (Me, dof_indices);
405  } // end of element loop
406 
407 
408 #else
409  // Avoid compiler warnings
410  libmesh_ignore(es);
411 #endif // LIBMESH_HAVE_SLEPC
412 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
void get_dirichlet_dofs(EquationSystems &es, const std::string &system_name, std::set< unsigned int > &global_dirichlet_dofs_set)
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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
UniquePtr< SparseMatrix< Number > > matrix_B
A second system matrix for generalized eigenvalue problems.
Definition: eigen_system.h:155
MeshBase & mesh
virtual void solve() libmesh_override
Override to solve the condensed eigenproblem with the dofs in local_non_condensed_dofs_vector strippe...
unsigned int get_n_converged() const
Definition: eigen_system.h:124
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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.
UniquePtr< SparseMatrix< Number > > matrix_A
The system matrix for standard eigenvalue problems.
Definition: eigen_system.h:150
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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.
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual SimpleRange< element_iterator > element_ptr_range()=0
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
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
void initialize_condensed_dofs(const std::set< dof_id_type > &global_condensed_dofs_set=std::set< dof_id_type >())
Loop over the dofs on each processor to initialize the list of non-condensed dofs.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
This class extends EigenSystem to allow a simple way of solving (standard or generalized) eigenvalue ...
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)
void libmesh_ignore(const T &)
void assemble_matrices(EquationSystems &es, const std::string &system_name)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:1793
UniquePtr< EigenSolver< Number > > eigen_solver
The EigenSolver, defining which interface, i.e solver package to use.
Definition: eigen_system.h:161
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
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 set_eigenproblem_type(EigenProblemType ept)
Sets the type of the current eigen problem.
Definition: eigen_system.C:78
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
unsigned int rank() const
Definition: parallel.h:724
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
This class provides a specific system class.
Definition: eigen_system.h:50
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i) libmesh_override
Override get_eigenpair() to retrieve the eigenpair for the condensed eigensolve.
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
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.