libMesh
Functions
eigenproblems_ex1.C File Reference

Go to the source code of this file.

Functions

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

Function Documentation

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

Referenced by main().

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

Definition at line 183 of file eigenproblems_ex1.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix(), dim, libMesh::DofMap::dof_indices(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::libmesh_ignore(), libMesh::EigenSystem::matrix_A, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseMatrix< T >::resize(), and libMesh::DofMap::variable_type().

185 {
186  // It is a good idea to make sure we are assembling
187  // the proper system.
188  libmesh_assert_equal_to (system_name, "Eigensystem");
189 
190 #ifdef LIBMESH_HAVE_SLEPC
191 
192  // Get a constant reference to the mesh object.
193  const MeshBase & mesh = es.get_mesh();
194 
195  // The dimension that we are running.
196  const unsigned int dim = mesh.mesh_dimension();
197 
198  // Get a reference to our system.
199  EigenSystem & eigen_system = es.get_system<EigenSystem> ("Eigensystem");
200 
201  // Get a constant reference to the Finite Element type
202  // for the first (and only) variable in the system.
203  FEType fe_type = eigen_system.get_dof_map().variable_type(0);
204 
205  // A reference to the system matrix
206  SparseMatrix<Number> & matrix_A = *eigen_system.matrix_A;
207 
208  // Build a Finite Element object of the specified type. Since the
209  // FEBase::build() member dynamically creates memory we will
210  // store the object as a UniquePtr<FEBase>. This can be thought
211  // of as a pointer that will clean up after itself.
212  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
213 
214  // A Gauss quadrature rule for numerical integration.
215  // Use the default quadrature order.
216  QGauss qrule (dim, fe_type.default_quadrature_order());
217 
218  // Tell the finite element object to use our quadrature rule.
219  fe->attach_quadrature_rule (&qrule);
220 
221  // The element Jacobian * quadrature weight at each integration point.
222  const std::vector<Real> & JxW = fe->get_JxW();
223 
224  // The element shape functions evaluated at the quadrature points.
225  const std::vector<std::vector<Real>> & phi = fe->get_phi();
226 
227  // A reference to the DofMap object for this system. The DofMap
228  // object handles the index translation from node and element numbers
229  // to degree of freedom numbers.
230  const DofMap & dof_map = eigen_system.get_dof_map();
231 
232  // The element mass matrix.
234 
235  // This vector will hold the degree of freedom indices for
236  // the element. These define where in the global system
237  // the element degrees of freedom get mapped.
238  std::vector<dof_id_type> dof_indices;
239 
240  // Now we will loop over all the elements in the mesh that
241  // live on the local processor. We will compute the element
242  // matrix and right-hand-side contribution. In case users
243  // later modify this program to include refinement, we will
244  // be safe and will only consider the active elements;
245  // hence we use a variant of the active_elem_iterator.
246  for (const auto & elem : mesh.active_local_element_ptr_range())
247  {
248  // Get the degree of freedom indices for the
249  // current element. These define where in the global
250  // matrix and right-hand-side this element will
251  // contribute to.
252  dof_map.dof_indices (elem, dof_indices);
253 
254  // Compute the element-specific data for the current
255  // element. This involves computing the location of the
256  // quadrature points (q_point) and the shape functions
257  // (phi, dphi) for the current element.
258  fe->reinit (elem);
259 
260  // Zero the element matrices and rhs before
261  // summing them. We use the resize member here because
262  // the number of degrees of freedom might have changed from
263  // the last element. Note that this will be the case if the
264  // element type is different (i.e. the last element was a
265  // triangle, now we are on a quadrilateral).
266  Me.resize (dof_indices.size(), dof_indices.size());
267 
268  // Now loop over the quadrature points. This handles
269  // the numeric integration.
270  //
271  // We will build the element matrix. This involves
272  // a double loop to integrate the test functions (i) against
273  // the trial functions (j).
274  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
275  for (std::size_t i=0; i<phi.size(); i++)
276  for (std::size_t j=0; j<phi.size(); j++)
277  Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp];
278 
279  // On an unrefined mesh, constrain_element_matrix does
280  // nothing. If this assembly function is ever repurposed to
281  // run on a refined mesh, getting the hanging node constraints
282  // right will be important. Note that, even with
283  // asymmetric_constraint_rows = false, the constrained dof
284  // diagonals still exist in the matrix, with diagonal entries
285  // that are there to ensure non-singular matrices for linear
286  // solves but which would generate positive non-physical
287  // eigenvalues for eigensolves.
288  dof_map.constrain_element_matrix(Me, dof_indices, false);
289 
290  // Finally, simply add the element contribution to the
291  // overall matrix.
292  matrix_A.add_matrix (Me, dof_indices);
293  } // end of element loop
294 
295 #else
296  // Avoid compiler warnings
297  libmesh_ignore(es);
298 #endif // LIBMESH_HAVE_SLEPC
299 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
MeshBase & mesh
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
void libmesh_ignore(const T &)
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
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.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
This class provides a specific system class.
Definition: eigen_system.h:50
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
int main ( int  argc,
char **  argv 
)

Definition at line 59 of file eigenproblems_ex1.C.

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_mass(), libMesh::System::attach_assemble_function(), libMesh::MeshTools::Generation::build_square(), libMesh::LibMeshInit::comm(), libMesh::FIRST, libMesh::EigenSystem::get_eigenpair(), libMesh::EigenSystem::get_n_converged(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::out, libMesh::EquationSystems::parameters, std::pow(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD4, libMesh::Real, libMesh::Parameters::set(), libMesh::EigenSystem::solve(), libMesh::TOLERANCE, and libMesh::MeshOutput< MT >::write_equation_systems().

60 {
61  // Initialize libMesh and the dependent libraries.
62  LibMeshInit init (argc, argv);
63 
64 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
65  // SLEPc currently gives us a nasty crash with Real==float
66  libmesh_example_requires(false, "--disable-singleprecision");
67 #endif
68 
69 #ifndef LIBMESH_HAVE_SLEPC
70  libmesh_example_requires(false, "--enable-slepc");
71 #else
72  // Check for proper usage.
73  if (argc < 3)
74  libmesh_error_msg("\nUsage: " << argv[0] << " -n <number of eigen values>");
75 
76  // Tell the user what we are doing.
77  else
78  {
79  libMesh::out << "Running " << argv[0];
80 
81  for (int i=1; i<argc; i++)
82  libMesh::out << " " << argv[i];
83 
84  libMesh::out << std::endl << std::endl;
85  }
86 
87  // Get the number of eigen values to be computed from argv[2]
88  const unsigned int nev = std::atoi(argv[2]);
89 
90  // Skip this 2D example if libMesh was compiled as 1D-only.
91  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
92 
93  // Create a mesh, with dimension to be overridden later, on the
94  // default MPI communicator.
95  Mesh mesh(init.comm());
96 
97  // Use the internal mesh generator to create a uniform
98  // 2D grid on a square.
100  20, 20,
101  -1., 1.,
102  -1., 1.,
103  QUAD4);
104 
105  // Print information about the mesh to the screen.
106  mesh.print_info();
107 
108  // Create an equation systems object.
109  EquationSystems equation_systems (mesh);
110 
111  // Create a EigenSystem named "Eigensystem" and (for convenience)
112  // use a reference to the system we create.
113  EigenSystem & eigen_system =
114  equation_systems.add_system<EigenSystem> ("Eigensystem");
115 
116  // Declare the system variables.
117  // Adds the variable "p" to "Eigensystem". "p"
118  // will be approximated using second-order approximation.
119  eigen_system.add_variable("p", FIRST);
120 
121  // Give the system a pointer to the matrix assembly
122  // function defined below.
124 
125  // Set necessary parameters used in EigenSystem::solve(),
126  // i.e. the number of requested eigenpairs nev and the number
127  // of basis vectors ncv used in the solution algorithm. Note that
128  // ncv >= nev must hold and ncv >= 2*nev is recommended.
129  equation_systems.parameters.set<unsigned int>("eigenpairs") = nev;
130  equation_systems.parameters.set<unsigned int>("basis vectors") = nev*3;
131 
132  // You may optionally change the default eigensolver used by SLEPc.
133  // The Krylov-Schur method is mathematically equivalent to implicitly
134  // restarted Arnoldi, the method of Arpack, so there is currently no
135  // point in using SLEPc with Arpack.
136  // ARNOLDI = default in SLEPc 2.3.1 and earlier
137  // KRYLOVSCHUR default in SLEPc 2.3.2 and later
138  // eigen_system.eigen_solver->set_eigensolver_type(KRYLOVSCHUR);
139 
140  // Set the solver tolerance and the maximum number of iterations.
141  equation_systems.parameters.set<Real>
142  ("linear solver tolerance") = pow(TOLERANCE, 5./3.);
143  equation_systems.parameters.set<unsigned int>
144  ("linear solver maximum iterations") = 1000;
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  // Solve the system "Eigensystem".
153  eigen_system.solve();
154 
155  // Get the number of converged eigen pairs.
156  unsigned int nconv = eigen_system.get_n_converged();
157 
158  libMesh::out << "Number of converged eigenpairs: " << nconv
159  << "\n" << std::endl;
160 
161  // Get the last converged eigenpair
162  if (nconv != 0)
163  {
164  eigen_system.get_eigenpair(nconv-1);
165 
166 #ifdef LIBMESH_HAVE_EXODUS_API
167  // Write the eigen vector to file.
168  ExodusII_IO (mesh).write_equation_systems ("out.e", equation_systems);
169 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
170  }
171  else
172  {
173  libMesh::out << "WARNING: Solver did not converge!\n" << nconv << std::endl;
174  }
175 
176  // All done.
177  return 0;
178 #endif // LIBMESH_HAVE_SLEPC
179 }
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:257
This is the EquationSystems class.
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
MeshBase & mesh
unsigned int get_n_converged() const
Definition: eigen_system.h:124
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
virtual void solve() libmesh_override
Assembles & solves the eigen system.
Definition: eigen_system.C:199
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.
double pow(double a, int b)
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
OStreamProxy out
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
void assemble_mass(EquationSystems &es, const std::string &system_name)
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448