]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-36: Apply SLEPcWrappers and give an error message when the matrices are not...
authorToby D. Young <tyoung@ippt.pan.pl>
Wed, 1 Jul 2009 06:13:38 +0000 (06:13 +0000)
committerToby D. Young <tyoung@ippt.pan.pl>
Wed, 1 Jul 2009 06:13:38 +0000 (06:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@19005 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-36/Makefile
deal.II/examples/step-36/doc/intro.dox
deal.II/examples/step-36/step-36.cc

index b753dea12385db29194f3a962ee7457d303dc6ec..9a5b39cbe7d7b4f78748bc9f8c555d92fe63b8fd 100644 (file)
@@ -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
 
 
 
index 928f5587391e882b86b45075dfa2cd1ae93c2384..5aa971f5c6d6a768502b36670182964d49575afa 100644 (file)
@@ -1,10 +1,14 @@
 <br>
 
-<i>This program was contributed by Toby Young and Wolfgang
-Bangerth. 
-</i>
+<i>This program was contributed by Toby D. Young and Wolfgang
+Bangerth.  </i>
 
 
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
+@f[
+       ( \frac{1}{2} \partial_x \partial_y + f(x,y) ) \Phi
+       = 
+       \varepsilon \Psi \quad,
+@f]
\ No newline at end of file
index c6df42a99b99f200f6600d1b5cff8b1d44eeb470..347b9931e9ee9cd87a070630acf3116817c305ef 100644 (file)
@@ -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 <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),
@@ -91,16 +100,16 @@ EigenvalueProblem<dim>::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 <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);
 
@@ -114,27 +123,29 @@ void EigenvalueProblem<dim>::make_grid_and_dofs ()
     .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" :
@@ -142,11 +153,13 @@ void EigenvalueProblem<dim>::assemble_system ()
                        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 ();
   
@@ -157,8 +170,8 @@ void EigenvalueProblem<dim>::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<dim>::assemble_system ()
          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
@@ -215,6 +254,7 @@ 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 ?
@@ -232,7 +272,6 @@ void EigenvalueProblem<dim>::output_results () const
 }
 
 
-
                                  // @sect4{EigenvalueProblem::run}
 
                                  // This is the function which has the
@@ -251,16 +290,20 @@ void EigenvalueProblem<dim>::run ()
 
 
                                  // @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)
     {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.