/* $Id$ */
-/* Author: Toby Young, Polish Academy of Sciences, */
+/* Author: Toby D. Young, Polish Academy of Sciences, */
/* Wolfgang Bangerth, Texas A&M University */
/* $Id$ */
/* */
/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
- // @sect3{Include files}
+ // @sect3{Include files}
#include <base/logstream.h>
#include <base/quadrature_lib.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
+ // We need the interfaces for solvers
+ // that SLEPc provides:
+#include <lac/slepc_solver.h>
+
+ // and then some standard C++:
#include <fstream>
#include <iostream>
-
- // The final step, as in previous
- // programs, is to import all the
- // deal.II class and function names
- // into the global namespace:
+ // and the final step, as in previous
+ // programs, is to import all the
+ // deal.II class and function names
+ // into the global namespace:
using namespace dealii;
- // @sect3{The <code>EigenvalueProblem</code> class template}
+ // @sect3{The
+ // <code>EigenvalueProblem</code>
+ // class template}
template <int dim>
class EigenvalueProblem
{
- public:
- EigenvalueProblem (const std::string &prm_file);
- void run ();
-
- private:
- void make_grid_and_dofs ();
- void assemble_system ();
- void solve ();
- void output_results () const;
-
- Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
-
- PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
- std::vector<PETScWrappers::Vector> eigenfunctions;
-
- ParameterHandler parameters;
+public:
+ EigenvalueProblem (const std::string &prm_file);
+ void run ();
+
+private:
+ void make_grid_and_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results () const;
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ // These are data types pertaining to
+ // the generalized problem.
+ PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
+ std::vector<PETScWrappers::Vector> eigenfunctions;
+ std::vector<double> eigenvalues;
+
+ ParameterHandler parameters;
};
+ // @sect3{Implementation of the
+ // <code>EigenvalueProblem</code>
+ // class}
-
- // @sect3{Implementation of the <code>EigenvalueProblem</code> class}
-
- // @sect4{EigenvalueProblem::EigenvalueProblem}
+ // @sect4{EigenvalueProblem::EigenvalueProblem}
template <int dim>
EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
- :
- fe (1),
- dof_handler (triangulation)
+ :
+ fe (1),
+ dof_handler (triangulation)
{
parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "10",
Patterns::Integer (0, 100),
parameters.declare_entry ("Potential", "0",
Patterns::Anything(),
"A functional description of the potential.");
-
+
parameters.read_input (prm_file);
}
- // @sect4{EigenvalueProblem::make_grid_and_dofs}
+ // @sect4{EigenvalueProblem::make_grid_and_dofs}
template <int dim>
void EigenvalueProblem<dim>::make_grid_and_dofs ()
{
- GridGenerator::hyper_cube (triangulation, -1, 1);
+ GridGenerator::hyper_cube (triangulation, -0.5, 0.5);
triangulation.refine_global (5);
dof_handler.distribute_dofs (fe);
.resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions"));
for (unsigned int i=0; i<eigenfunctions.size(); ++i)
eigenfunctions[i].reinit (dof_handler.n_dofs());
+ eigenvalues
+ .resize (eigenfunctions.size());
}
- // @sect4{EigenvalueProblem::assemble_system}
+ // @sect4{EigenvalueProblem::assemble_system}
template <int dim>
void EigenvalueProblem<dim>::assemble_system ()
{
- QGauss<dim> quadrature_formula(2);
-
+ QGauss<dim> quadrature_formula(2);
+
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
-
+
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
-
+
FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
-
+
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
+
FunctionParser<dim> potential;
potential.initialize ((dim == 2 ?
"x,y" :
parameters.get ("Potential"),
typename FunctionParser<dim>::ConstMap());
std::vector<double> potential_values (n_q_points);
-
+
+ // Initialize objects denoting zero
+ // boundary constraints for the
+ // present grid.
ConstraintMatrix constraints;
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(),
+ constraints.clear ();
+ DoFTools::make_zero_boundary_constraints (dof_handler,
constraints);
constraints.close ();
{
fe_values.reinit (cell);
cell_stiffness_matrix = 0;
- cell_mass_matrix = 0;
-
+ cell_mass_matrix = 0;
+
potential.value_list (fe_values.get_quadrature_points(),
potential_values);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_stiffness_matrix(i,j)
- += ((fe_values.shape_grad (i, q_point) *
- fe_values.shape_grad (j, q_point)
+ += ((0.5 *
+ fe_values.shape_grad (i, q_point) *
+ fe_values.shape_grad (j, q_point)
+
- potential_values[q_point] *
- fe_values.shape_value (i, q_point) *
+ // potential_values[q_point] *
+ fe_values.shape_value (i, q_point) *
fe_values.shape_value (j, q_point)) *
fe_values.JxW (q_point));
+
cell_mass_matrix(i,j)
+= (fe_values.shape_value (i, q_point) *
fe_values.shape_value (j, q_point) *
fe_values.JxW (q_point));
}
-
- cell->get_dof_indices (local_dof_indices);
- constraints.distribute_local_to_global (cell_stiffness_matrix,
- local_dof_indices,
- stiffness_matrix);
- constraints.distribute_local_to_global (cell_mass_matrix,
- local_dof_indices,
- mass_matrix);
+ // Now we have the local system we
+ // transfer it into the global objects
+ // and take care of zero boundary
+ // constraints,
+ cell->get_dof_indices (local_dof_indices);
+
+ constraints
+ .distribute_local_to_global (cell_stiffness_matrix,
+ local_dof_indices,
+ stiffness_matrix);
+ constraints
+ .distribute_local_to_global (cell_mass_matrix,
+ local_dof_indices,
+ mass_matrix);
}
+
+ // and finally, set the matrices and
+ // vectors in an assembled state.
+ stiffness_matrix.compress();
+ mass_matrix.compress();
}
- // @sect4{EigenvalueProblem::solve}
+ // @sect4{EigenvalueProblem::solve}
template <int dim>
void EigenvalueProblem<dim>::solve ()
{
-// do whatever is necessary here. use stiffness_matrix, mass_matrix, and the
-// eigenfunctions[] array
+ // assign the accuracy to which we
+ // would like to solve the system,
+ SolverControl solver_control (1000, 1e-6);
+ SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
+
+ // choose which part of the spectrum to
+ // solve for
+
+ // eigensolver.set_which_eigenpairs ();
+
+ // then actually solve the system,
+
+ eigensolver.solve (stiffness_matrix, mass_matrix,
+ eigenvalues, eigenfunctions,
+ eigenfunctions.size());
}
- // @sect4{EigenvalueProblem::output_results}
+ // @sect4{EigenvalueProblem::output_results}
template <int dim>
void EigenvalueProblem<dim>::output_results () const
std::string("solution") +
Utilities::int_to_string(i));
+ // How does this work?
Vector<double> projected_potential (dof_handler.n_dofs());
FunctionParser<dim> potential;
potential.initialize ((dim == 2 ?
}
-
// @sect4{EigenvalueProblem::run}
// This is the function which has the
// @sect3{The <code>main</code> function}
-
int main (int argc, char **argv)
{
try
{
- PetscInitialize(&argc,&argv,0,0);
- deallog.depth_console (0);
-
- EigenvalueProblem<2> problem ("step-36.prm");
- problem.run ();
+ SlepcInitialize(&argc,&argv,0,0);
+
+ {
+ deallog.depth_console (0);
+
+ EigenvalueProblem<2> problem ("step-36.prm");
+ problem.run ();
+ }
+
+ SlepcFinalize();
}
catch (std::exception &exc)
{