From 784220b721bc5132de7f097b62ff38733bab4151 Mon Sep 17 00:00:00 2001 From: "Toby D. Young" Date: Wed, 1 Jul 2009 06:13:38 +0000 Subject: [PATCH] Step-36: Apply SLEPcWrappers and give an error message when the matrices are not assembled correctly git-svn-id: https://svn.dealii.org/trunk@19005 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-36/Makefile | 2 +- deal.II/examples/step-36/doc/intro.dox | 10 +- deal.II/examples/step-36/step-36.cc | 187 +++++++++++++++---------- 3 files changed, 123 insertions(+), 76 deletions(-) diff --git a/deal.II/examples/step-36/Makefile b/deal.II/examples/step-36/Makefile index b753dea123..9a5b39cbe7 100644 --- a/deal.II/examples/step-36/Makefile +++ b/deal.II/examples/step-36/Makefile @@ -30,7 +30,7 @@ D = ../../ # shall be deleted when calling `make clean'. Object and backup files, # executables and the like are removed anyway. Here, we give a list of # files in the various output formats that deal.II supports. -clean-up-files = *gmv *gnuplot *gpl *eps *pov +clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk diff --git a/deal.II/examples/step-36/doc/intro.dox b/deal.II/examples/step-36/doc/intro.dox index 928f558739..5aa971f5c6 100644 --- a/deal.II/examples/step-36/doc/intro.dox +++ b/deal.II/examples/step-36/doc/intro.dox @@ -1,10 +1,14 @@
-This program was contributed by Toby Young and Wolfgang -Bangerth. - +This program was contributed by Toby D. Young and Wolfgang +Bangerth.

Introduction

+@f[ + ( \frac{1}{2} \partial_x \partial_y + f(x,y) ) \Phi + = + \varepsilon \Psi \quad, +@f] \ No newline at end of file diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc index c6df42a99b..347b9931e9 100644 --- a/deal.II/examples/step-36/step-36.cc +++ b/deal.II/examples/step-36/step-36.cc @@ -1,5 +1,5 @@ /* $Id$ */ -/* Author: Toby Young, Polish Academy of Sciences, */ +/* Author: Toby D. Young, Polish Academy of Sciences, */ /* Wolfgang Bangerth, Texas A&M University */ /* $Id$ */ /* */ @@ -10,7 +10,7 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ - // @sect3{Include files} + // @sect3{Include files} #include #include @@ -37,52 +37,61 @@ #include #include + // We need the interfaces for solvers + // that SLEPc provides: +#include + + // and then some standard C++: #include #include - - // 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 EigenvalueProblem class template} + // @sect3{The + // EigenvalueProblem + // class template} template 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 triangulation; - FE_Q fe; - DoFHandler dof_handler; - - PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; - std::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 triangulation; + FE_Q fe; + DoFHandler dof_handler; + + // These are data types pertaining to + // the generalized problem. + PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; + std::vector eigenfunctions; + std::vector eigenvalues; + + ParameterHandler parameters; }; + // @sect3{Implementation of the + // EigenvalueProblem + // class} - - // @sect3{Implementation of the EigenvalueProblem class} - - // @sect4{EigenvalueProblem::EigenvalueProblem} + // @sect4{EigenvalueProblem::EigenvalueProblem} template EigenvalueProblem::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), @@ -91,16 +100,16 @@ EigenvalueProblem::EigenvalueProblem (const std::string &prm_file) 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 void EigenvalueProblem::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); @@ -114,27 +123,29 @@ void EigenvalueProblem::make_grid_and_dofs () .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions")); for (unsigned int i=0; i void EigenvalueProblem::assemble_system () { - QGauss quadrature_formula(2); - + QGauss quadrature_formula(2); + FEValues 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 cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); - + std::vector local_dof_indices (dofs_per_cell); - + FunctionParser potential; potential.initialize ((dim == 2 ? "x,y" : @@ -142,11 +153,13 @@ void EigenvalueProblem::assemble_system () parameters.get ("Potential"), typename FunctionParser::ConstMap()); std::vector 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(), + constraints.clear (); + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); constraints.close (); @@ -157,8 +170,8 @@ void EigenvalueProblem::assemble_system () { 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); @@ -167,41 +180,67 @@ void EigenvalueProblem::assemble_system () for (unsigned int j=0; jget_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 void EigenvalueProblem::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 void EigenvalueProblem::output_results () const @@ -215,6 +254,7 @@ void EigenvalueProblem::output_results () const std::string("solution") + Utilities::int_to_string(i)); + // How does this work? Vector projected_potential (dof_handler.n_dofs()); FunctionParser potential; potential.initialize ((dim == 2 ? @@ -232,7 +272,6 @@ void EigenvalueProblem::output_results () const } - // @sect4{EigenvalueProblem::run} // This is the function which has the @@ -251,16 +290,20 @@ void EigenvalueProblem::run () // @sect3{The main 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) { -- 2.39.5