*/\r
\r
#include <deal.II/base/logstream.h>\r
+#include <deal.II/base/table_handler.h>\r
#include <deal.II/base/quadrature_lib.h>\r
#include <deal.II/base/function.h>\r
#include <deal.II/base/utilities.h>\r
\r
PETScWrappers::SparseMatrix A, B;\r
std::vector<PETScWrappers::Vector> x;\r
- std::vector<double> lambda;\r
- \r
- ConstraintMatrix constraints;\r
+ std::vector<double> lambda; \r
+ ConstraintMatrix constraints;\r
+\r
+ TableHandler output_table;\r
};\r
\r
LaplaceEigenspectrumProblem::LaplaceEigenspectrumProblem ()\r
\r
void LaplaceEigenspectrumProblem::setup_system ()\r
{\r
- GridGenerator::hyper_cube (triangulation, -1, 1);\r
- triangulation.refine_global (3);\r
dof_handler.distribute_dofs (fe);\r
\r
+ constraints.clear ();\r
DoFTools::make_zero_boundary_constraints (dof_handler, constraints);\r
constraints.close ();\r
\r
x[0].reinit (dof_handler.n_dofs ());\r
lambda.resize (1);\r
lambda[0] = 0.;\r
+\r
+ // some output\r
+ output_table.add_value ("cells", triangulation.n_active_cells());\r
+ output_table.add_value ("dofs", dof_handler.n_dofs());\r
}\r
\r
void LaplaceEigenspectrumProblem::assemble_system ()\r
\r
void LaplaceEigenspectrumProblem::solve ()\r
{\r
- SolverControl solver_control (1., 1e-03);\r
- SLEPcWrappers::SolverLAPACK eigensolver (solver_control);\r
+ SolverControl solver_control (1000, 1e-03);\r
+ SLEPcWrappers::SolverArnoldi eigensolver (solver_control);\r
+ eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);\r
eigensolver.solve (A, B, lambda, x, x.size());\r
+\r
+ // some output\r
+ output_table.add_value ("lambda", lambda[0]);\r
+ output_table.add_value ("error", std::fabs(2.-lambda[0]));\r
}\r
\r
void LaplaceEigenspectrumProblem::run ()\r
{\r
- setup_system ();\r
- assemble_system ();\r
- solve ();\r
+ const double radius = dealii::numbers::PI/2.;\r
+ GridGenerator::hyper_cube (triangulation, -radius, radius);\r
+\r
+ // set the old eigenvalue to a silly number.\r
+ double old_lambda = 1000;\r
+\r
+ for (unsigned int c=0; c<5; ++c)\r
+ {\r
+ // obtain numerical result\r
+ triangulation.refine_global (1);\r
+ setup_system ();\r
+ assemble_system ();\r
+ solve ();\r
+\r
+ // check energy convergence with previous result\r
+ AssertThrow (lambda[0]<old_lambda, ExcMessage("solution is not converging"));\r
+ old_lambda = lambda[0];\r
+ }\r
+\r
+ // push back analytic result\r
+ output_table.add_value ("cells", "inf");\r
+ output_table.add_value ("dofs", "inf");\r
+ output_table.add_value ("lambda", 2.);\r
+ output_table.add_value ("error", "-");\r
+\r
+ // finialise output\r
+ output_table.write_text (std::cout);\r
+ deallog << std::endl;\r
}\r
\r
\r
deallog.depth_console (0);\r
LaplaceEigenspectrumProblem problem;\r
problem.run ();\r
- deallog << "OK" << std::endl;\r
}\r
}\r
\r