From 85caea0a100c0259466e78768fa1013fbc724c16 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Sat, 22 Nov 2014 07:44:54 +0100
Subject: [PATCH] added a step-36 test to arpack category

---
 tests/CMakeLists.txt           |   2 +-
 tests/arpack/step-36_ar.cc     | 374 +++++++++++++++++++++++++++++++++
 tests/arpack/step-36_ar.output |   5 +
 3 files changed, 380 insertions(+), 1 deletion(-)
 create mode 100644 tests/arpack/step-36_ar.cc
 create mode 100644 tests/arpack/step-36_ar.output

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 35bfe8b82d..a59fcf8108 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -62,7 +62,7 @@ SET(_categories
   a-framework algorithms all-headers aniso base bits build_tests codim_one
   deal.II distributed_grids fe gla grid hp integrators lac lapack manifold
   matrix_free mesh_converter metis mpi multigrid mumps opencascade petsc serialization
-  slepc trilinos umfpack
+  slepc trilinos umfpack arpack
   )
 
 #
diff --git a/tests/arpack/step-36_ar.cc b/tests/arpack/step-36_ar.cc
new file mode 100644
index 0000000000..db792017de
--- /dev/null
+++ b/tests/arpack/step-36_ar.cc
@@ -0,0 +1,374 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2014 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Authors: Toby D. Young, Polish Academy of Sciences,
+ *          Wolfgang Bangerth, Texas A&M University
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/lac/arpack_solver.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+
+#include <complex>
+
+#include <fstream>
+#include <iostream>
+
+namespace Step36
+{
+  using namespace dealii;
+
+
+  template <int dim>
+  class EigenvalueProblem
+  {
+  public:
+    EigenvalueProblem (const std::string &prm_file);
+    void run ();
+
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    std::pair<unsigned int, double> solve ();
+    void output_results () const;
+
+    Triangulation<dim> triangulation;
+    FE_Q<dim>          fe;
+    DoFHandler<dim>    dof_handler;
+
+    SparsityPattern                     sparsity_pattern;
+    SparseMatrix<double>                stiffness_matrix, mass_matrix;
+    std::vector<Vector<double> >        eigenfunctions;
+    std::vector<std::complex<double>>        eigenvalues;
+
+    //ParameterHandler parameters;
+
+    ConstraintMatrix constraints;
+  };
+
+
+
+  template <int dim>
+  EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
+    :
+    fe (1),
+    dof_handler (triangulation)
+  {
+//     parameters.declare_entry ("Global mesh refinement steps", "5",
+//                               Patterns::Integer (0, 20),
+//                               "The number of times the 1-cell coarse mesh should "
+//                               "be refined globally for our computations.");
+//     parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "5",
+//                               Patterns::Integer (0, 100),
+//                               "The number of eigenvalues/eigenfunctions "
+//                               "to be computed.");
+//     parameters.declare_entry ("Potential", "0",
+//                               Patterns::Anything(),
+//                               "A functional description of the potential.");
+
+//     parameters.read_input (prm_file);
+  }
+
+
+
+  template <int dim>
+  void EigenvalueProblem<dim>::make_grid_and_dofs ()
+  {
+    GridGenerator::hyper_cube (triangulation, -1, 1);
+    triangulation.refine_global (5/*parameters.get_integer ("Global mesh refinement steps")*/);
+    dof_handler.distribute_dofs (fe);
+
+    DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
+    constraints.close ();
+
+    sparsity_pattern.reinit (dof_handler.n_dofs(),
+                             dof_handler.n_dofs(),
+                             dof_handler.max_couplings_between_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+    constraints.condense (sparsity_pattern);
+    sparsity_pattern.compress();
+    stiffness_matrix.reinit (sparsity_pattern);
+    mass_matrix.reinit (sparsity_pattern);
+
+    eigenfunctions
+    .resize (5/*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 ());
+  }
+
+
+
+  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_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<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    FunctionParser<dim> potential;
+    potential.initialize (FunctionParser<dim>::default_variable_names (),
+                          "0"/*parameters.get ("Potential")*/,
+                          typename FunctionParser<dim>::ConstMap());
+
+    std::vector<double> potential_values (n_q_points);
+
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active (),
+    endc = dof_handler.end ();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+        cell_stiffness_matrix = 0;
+        cell_mass_matrix      = 0;
+
+        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)
+              {
+                cell_stiffness_matrix (i, j)
+                += (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_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);
+      }
+
+    stiffness_matrix.compress (VectorOperation::add);
+    mass_matrix.compress (VectorOperation::add);
+
+
+    double min_spurious_eigenvalue = std::numeric_limits<double>::max(),
+           max_spurious_eigenvalue = -std::numeric_limits<double>::max();
+
+    for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+      if (constraints.is_constrained(i))
+        {
+          const double ev = stiffness_matrix(i,i)/mass_matrix(i,i);
+          min_spurious_eigenvalue = std::min (min_spurious_eigenvalue, ev);
+          max_spurious_eigenvalue = std::max (max_spurious_eigenvalue, ev);
+        }
+
+//     std::cout << "   Spurious eigenvalues are all in the interval "
+//               << "[" << min_spurious_eigenvalue << "," << max_spurious_eigenvalue << "]"
+//               << std::endl;
+
+  }
+
+
+
+  template <int dim>
+  std::pair<unsigned int, double> EigenvalueProblem<dim>::solve ()
+  {
+    SolverControl solver_control (dof_handler.n_dofs(), 1e-9);
+    SparseDirectUMFPACK inverse;
+    inverse.initialize (stiffness_matrix);
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2;
+    ArpackSolver::AdditionalData additional_data(num_arnoldi_vectors,
+                                                 ArpackSolver::largest_real_part);
+    ArpackSolver eigensolver (solver_control, additional_data);
+    eigensolver.solve (stiffness_matrix,
+                       mass_matrix,
+                       inverse,
+                       eigenvalues,
+                       eigenfunctions,
+                       eigenvalues.size());
+
+    for (unsigned int i=0; i<eigenfunctions.size(); ++i)
+      eigenfunctions[i] /= eigenfunctions[i].linfty_norm ();
+
+    return std::make_pair (solver_control.last_step (),
+                           solver_control.last_value ());
+  }
+
+
+
+  template <int dim>
+  void EigenvalueProblem<dim>::output_results () const
+  {
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler (dof_handler);
+
+    for (unsigned int i=0; i<eigenfunctions.size(); ++i)
+      data_out.add_data_vector (eigenfunctions[i],
+                                std::string("eigenfunction_") +
+                                Utilities::int_to_string(i));
+
+    Vector<double> projected_potential (dof_handler.n_dofs());
+    {
+      FunctionParser<dim> potential;
+      potential.initialize (FunctionParser<dim>::default_variable_names (),
+                            "0"/*parameters.get ("Potential")*/,
+                            typename FunctionParser<dim>::ConstMap());
+      VectorTools::interpolate (dof_handler, potential, projected_potential);
+    }
+    data_out.add_data_vector (projected_potential, "interpolated_potential");
+
+    data_out.build_patches ();
+
+    std::ofstream output ("eigenvectors.vtk");
+    data_out.write_vtk (output);
+  }
+
+
+
+  template <int dim>
+  void EigenvalueProblem<dim>::run ()
+  {
+    make_grid_and_dofs ();
+
+//     std::cout << "   Number of active cells:       "
+//               << triangulation.n_active_cells ()
+//               << std::endl
+//               << "   Number of degrees of freedom: "
+//               << dof_handler.n_dofs ()
+//               << std::endl;
+
+    assemble_system ();
+
+    const std::pair<unsigned int, double> res = solve ();
+//     std::cout << "   Solver converged in "   << res.first
+//               << " iterations. Residual is " << res.second << std::endl;
+
+//     output_results ();
+// 
+//     std::cout << std::endl;
+//     for (unsigned int i=0; i<eigenvalues.size(); ++i)
+//       std::cout << "      Eigenvalue " << i
+//                 << " : " << eigenvalues[i]
+//                 << std::endl;
+
+    for (unsigned int i=0; i<eigenvalues.size(); ++i)
+      deallog << "      Eigenvalue " << i
+              << " : " << eigenvalues[i]
+              << std::endl;
+  }
+}
+
+int main (int argc, char **argv)
+{
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step36;
+      
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        deallog.depth_console (0);
+
+        EigenvalueProblem<2> problem ("");
+        problem.run ();
+      }
+    }
+
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+//   std::cout << std::endl
+//             << "   Job done."
+//             << std::endl;
+
+  return 0;
+}
diff --git a/tests/arpack/step-36_ar.output b/tests/arpack/step-36_ar.output
new file mode 100644
index 0000000000..6c321a9872
--- /dev/null
+++ b/tests/arpack/step-36_ar.output
@@ -0,0 +1,5 @@
+DEAL::      Eigenvalue 0 : (4.93877,0.00000)
+DEAL::      Eigenvalue 1 : (19.8027,0.00000)
+DEAL::      Eigenvalue 2 : (12.3707,0.00000)
+DEAL::      Eigenvalue 3 : (12.3707,0.00000)
+DEAL::      Eigenvalue 4 : (24.8370,0.00000)
-- 
2.39.5