/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
- // @sect3{Include files}
-
+ // @sect3{Include files} We bundle the
+ // "usual" deal.II include files as we
+ // did in step-4:
#include <base/logstream.h>
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <base/parameter_handler.h>
#include <base/utilities.h>
#include <lac/full_matrix.h>
-#include <lac/petsc_sparse_matrix.h>
-#include <lac/petsc_vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
#include <lac/compressed_simple_sparsity_pattern.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
- // We need the interfaces for solvers
- // that SLEPc provides:
+ // PETSc appears here because SLEPc
+ // depends on him:
+#include <lac/petsc_sparse_matrix.h>
+#include <lac/petsc_vector.h>
+
+ // and then we need to actually import
+ // the interfaces for solvers that
+ // SLEPc provides:
#include <lac/slepc_solver.h>
- // and then some standard C++:
+ // and some standard C++:
#include <fstream>
#include <iostream>
- // and the final step, as in previous
- // programs, is to import all the
- // deal.II class and function names
- // into the global namespace:
+ // and the finally, as in previous
+ // programs, we import all the deal.II
+ // class and function names into the
+ // global namespace:
using namespace dealii;
- // @sect3{The
+ // @sect1{The
// <code>EigenvalueProblem</code>
// class template}
void solve ();
void output_results () const;
- Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
// These are data types pertaining to
// the generalized problem.
ParameterHandler parameters;
};
- // @sect3{Implementation of the
+ // @sect2{Implementation of the
// <code>EigenvalueProblem</code>
// class}
- // @sect4{EigenvalueProblem::EigenvalueProblem}
+ // @sect3{EigenvalueProblem::EigenvalueProblem}
template <int dim>
EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
fe (1),
dof_handler (triangulation)
{
+
+ // Declare some of the needed data
+ // from a file; you can always change
+ // this!
parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "10",
Patterns::Integer (0, 100),
"The number of eigenvalues/eigenfunctions "
Patterns::Anything(),
"A functional description of the potential.");
+ // Entries are declared, so now we
+ // read them in...
parameters.read_input (prm_file);
}
triangulation.refine_global (5);
dof_handler.distribute_dofs (fe);
- CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
- dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler, csp);
- stiffness_matrix.reinit (csp);
- mass_matrix.reinit (csp);
+ // CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+ // dof_handler.n_dofs());
+ // DoFTools::make_sparsity_pattern (dof_handler, csp);
+ // stiffness_matrix.reinit (csp);
+ // mass_matrix.reinit (csp);
+ stiffness_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ mass_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+
+ // Set the collective eigenfunction
+ // block to be as big as we wanted!
eigenfunctions
.resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions"));
- for (unsigned int i=0; i<eigenfunctions.size(); ++i)
- eigenfunctions[i].reinit (dof_handler.n_dofs());
+ for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+ eigenfunctions[i].reinit (dof_handler.n_dofs ());
+
+ // and do the same for the eigenvalue
+ // output, which had better be the
+ // same size as the eigenfunction
+ // block
eigenvalues
- .resize (eigenfunctions.size());
+ .resize (eigenfunctions.size ());
}
- // @sect4{EigenvalueProblem::assemble_system}
+ // @sect5{EigenvalueProblem::assemble_system}
template <int dim>
void EigenvalueProblem<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
FEValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
+ update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
"x,y,z"),
parameters.get ("Potential"),
typename FunctionParser<dim>::ConstMap());
+
+ // Here we call the actual
+ // quadrature-point values of the
+ // potential defined in the prm file
+ // which we read in earlier
std::vector<double> potential_values (n_q_points);
// Initialize objects denoting zero
// present grid.
ConstraintMatrix constraints;
constraints.clear ();
- DoFTools::make_zero_boundary_constraints (dof_handler,
- constraints);
+ DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
constraints.close ();
typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ cell = dof_handler.begin_active (),
+ endc = dof_handler.end ();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
potential.value_list (fe_values.get_quadrature_points(),
potential_values);
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int q_point=0; q_point<n_q_points;
+ ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell;
+ ++i)
+ for (unsigned int j=0; j<dofs_per_cell;
+ ++j)
{
- cell_stiffness_matrix(i,j)
- += ((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) *
- fe_values.shape_value (j, q_point)) *
- fe_values.JxW (q_point));
+ cell_stiffness_matrix (i, j)
+ += (fe_values.shape_grad (i, q_point) *
+ 0.5 *
+ fe_values.shape_grad (j, q_point)
+ +
+ fe_values.shape_value (i, q_point) *
+ // potential_values[q_point] *
+ 0 * // infinite potential well
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
- cell_mass_matrix(i,j)
+ cell_mass_matrix (i, j)
+= (fe_values.shape_value (i, q_point) *
- fe_values.shape_value (j, q_point) *
- fe_values.JxW (q_point));
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
}
// Now we have the local system we
stiffness_matrix.compress();
mass_matrix.compress();
- // make sure that the diagonal entries of
- // constrained degrees of freedom are
- // non-zero to ensure that the matrix is
- // actually invertible
- for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
- if (constraints.is_constrained(i))
- {
- stiffness_matrix.set (i, i, 1);
- mass_matrix.set (i, i, 1);
- }
+ // make sure that the diagonal entries
+ // of constrained degrees of freedom
+ // are non-zero to ensure that the
+ // matrix is actually invertible
+ // for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ // if (constraints.is_constrained(i))
+ // {
+ // stiffness_matrix.set (i, i, 1);
+ // mass_matrix.set (i, i, 1);
+ // }
- stiffness_matrix.compress();
- mass_matrix.compress();
+ // finally set the matrices in an
+ // assembled state so SLEPc likes:
+ // stiffness_matrix.compress ();
+ // mass_matrix.compress ();
}
- // @sect4{EigenvalueProblem::solve}
+ // @sect6{EigenvalueProblem::solve}
+ // Now that the system is set up, here
+ // is a good time to actually solve
+ // the problem: As with other examples
+ // this is done using a "solve"
+ // routine
template <int dim>
void EigenvalueProblem<dim>::solve ()
{
- // assign the accuracy to which we
- // would like to solve the system,
+ // We start by assigning 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
+ // and assign our solver of
+ // choice. Here we want to use the
+ // Krylov-Schur solver, which is
+ // pretty darn fast and robust:
+ SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
+ // Lets assign the solver which part
+ // of the spectrum we want to solve
eigensolver.set_which_eigenpairs (EPS_SMALLEST_MAGNITUDE);
- // then actually solve the system,
-
+ // Finally, we actually solve the
+ // generalized eigenproblem:
eigensolver.solve (stiffness_matrix, mass_matrix,
eigenvalues, eigenfunctions,
eigenfunctions.size());
}
- // @sect4{EigenvalueProblem::output_results}
-
+ // @sect7{EigenvalueProblem::output_results}
template <int dim>
void EigenvalueProblem<dim>::output_results () const
{
std::ofstream output ("eigenvectors.vtk");
data_out.write_vtk (output);
+
+ for (unsigned int i=0; i<eigenvalues.size(); ++i)
+ std::cout << std::endl
+ << " eigenvalue " << i
+ << " : " << eigenvalues[i];
+
}
- // @sect4{EigenvalueProblem::run}
+ // @sect8{EigenvalueProblem::run}
// This is the function which has the
- // top-level control over everything. It is
- // the same as for the previous example.
+ // top-level control over
+ // everything. It is very similar as
+ // for the previous examples.
template <int dim>
void EigenvalueProblem<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
- make_grid_and_dofs();
+ make_grid_and_dofs ();
+
+ // While we are here, lets count the
+ // number of active cells and degrees
+ // of freedom like we always do.
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs();
+
assemble_system ();
solve ();
output_results ();
}
- // @sect3{The <code>main</code> function}
+ // @sect9{The <code>main</code> function}
int main (int argc, char **argv)
{
try
{
- SlepcInitialize(&argc,&argv,0,0);
+
+ // Here is another difference from
+ // other steps: We initialize the
+ // SLEPc work space which inherently
+ // initializes the PETSc work space
+ SlepcInitialize (&argc,&argv,0,0);
{
deallog.depth_console (0);
problem.run ();
}
- SlepcFinalize();
+ // and then unitialize the SLEPc
+ // work space when the job is done:
+ SlepcFinalize ();
}
+
+ // or panic if something goes wrong.
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
return 1;
}
+ // or show that we happy and didn't
+ // crap out on the calculation...
+ std::cout << std::endl
+ << "Completed"
+ << std::endl;
+
return 0;
}