--- /dev/null
+/* ---------------------------------------------------------------------\r
+ *\r
+ * Copyright (C) 2013 - 2015 by the deal.II authors\r
+ *\r
+ * This file is part of the deal.II library.\r
+ *\r
+ * The deal.II library is free software; you can use it, redistribute\r
+ * it, and/or modify it under the terms of the GNU Lesser General\r
+ * Public License as published by the Free Software Foundation; either\r
+ * version 2.1 of the License, or (at your option) any later version.\r
+ * The full text of the license can be found in the file LICENSE at\r
+ * the top level of the deal.II distribution.\r
+ *\r
+ * ---------------------------------------------------------------------\r
+ */\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
+#include <deal.II/grid/tria.h>\r
+#include <deal.II/grid/grid_generator.h>\r
+#include <deal.II/grid/tria_accessor.h>\r
+#include <deal.II/grid/tria_iterator.h>\r
+#include <deal.II/dofs/dof_handler.h>\r
+#include <deal.II/dofs/dof_accessor.h>\r
+#include <deal.II/dofs/dof_tools.h>\r
+#include <deal.II/fe/fe_q.h>\r
+#include <deal.II/fe/fe_values.h>\r
+#include <deal.II/numerics/vector_tools.h>\r
+#include <deal.II/numerics/matrix_tools.h>\r
+#include <deal.II/numerics/data_out.h>\r
+#include <deal.II/lac/full_matrix.h>\r
+#include <deal.II/lac/trilinos_sparse_matrix.h>\r
+#include <deal.II/lac/trilinos_vector.h>\r
+#include <deal.II/lac/trilinos_solver.h>\r
+#include <deal.II/lac/trilinos_precondition.h>\r
+#include <deal.II/lac/dynamic_sparsity_pattern.h>\r
+\r
+#include <fstream>\r
+#include <iostream>\r
+\r
+using namespace dealii;\r
+\r
+// Test that deal.II is working with Trilinos by solving the Laplace's\r
+// problem in 2d.\r
+class LaplaceProblem\r
+{\r
+public:\r
+ LaplaceProblem ();\r
+ void run ();\r
+\r
+private:\r
+ void setup_system ();\r
+ void assemble_system ();\r
+ void solve ();\r
+\r
+ Triangulation<2> triangulation;\r
+ FE_Q<2> fe;\r
+ DoFHandler<2> dof_handler;\r
+\r
+ TrilinosWrappers::SparseMatrix A;\r
+ TrilinosWrappers::Vector b, x;\r
+ ConstraintMatrix constraints;\r
+\r
+ TableHandler output_table;\r
+};\r
+\r
+LaplaceProblem::LaplaceProblem ()\r
+ :\r
+ fe (1),\r
+ dof_handler (triangulation)\r
+{}\r
+\r
+void LaplaceProblem::setup_system ()\r
+{\r
+ dof_handler.distribute_dofs (fe);\r
+\r
+ constraints.clear ();\r
+ DoFTools::make_zero_boundary_constraints (dof_handler, constraints);\r
+ constraints.close ();\r
+\r
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());\r
+ DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);\r
+\r
+ A.reinit (dsp);\r
+ b.reinit (dof_handler.n_dofs());\r
+ x.reinit (dof_handler.n_dofs());\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 LaplaceProblem::assemble_system ()\r
+{\r
+ QGauss<2> quadrature_formula(2);\r
+\r
+ FEValues<2> fe_values (fe, quadrature_formula,\r
+ update_values |\r
+ update_gradients |\r
+ update_quadrature_points |\r
+ update_JxW_values);\r
+\r
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;\r
+ const unsigned int n_q_points = quadrature_formula.size();\r
+\r
+ FullMatrix<double> cell_A (dofs_per_cell, dofs_per_cell);\r
+ Vector<double> cell_b (dofs_per_cell);\r
+\r
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);\r
+\r
+ DoFHandler<2>::active_cell_iterator\r
+ cell = dof_handler.begin_active (),\r
+ endc = dof_handler.end ();\r
+\r
+ for (; cell!=endc; ++cell)\r
+ {\r
+ fe_values.reinit (cell);\r
+ cell_A = 0;\r
+ cell_b = 0;\r
+\r
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)\r
+ for (unsigned int i=0; i<dofs_per_cell; ++i)\r
+ {\r
+ for (unsigned int j=0; j<dofs_per_cell; ++j)\r
+ {\r
+ cell_A (i, j)\r
+ +=\r
+ fe_values.shape_grad (i, q_point) *\r
+ fe_values.shape_grad (j, q_point)\r
+ *\r
+ fe_values.JxW (q_point);\r
+ }\r
+\r
+ cell_b (i)\r
+ +=\r
+ fe_values.shape_value (i, q_point)\r
+ *\r
+ fe_values.JxW (q_point);\r
+ }\r
+\r
+ cell->get_dof_indices (local_dof_indices);\r
+\r
+ constraints.distribute_local_to_global (cell_A, local_dof_indices, A);\r
+ constraints.distribute_local_to_global (cell_b, local_dof_indices, b);\r
+ }\r
+\r
+ A.compress (VectorOperation::add);\r
+ b.compress (VectorOperation::add);\r
+}\r
+\r
+void LaplaceProblem::solve ()\r
+{\r
+ SolverControl solver_control (1e03, 1e-03);\r
+ TrilinosWrappers::SolverCG cg_solver (solver_control);\r
+ TrilinosWrappers::PreconditionBlockJacobi preconditioner;\r
+ preconditioner.initialize (A);\r
+ cg_solver.solve (A, x, b, preconditioner);\r
+\r
+ // some output\r
+ // ?\r
+}\r
+\r
+void LaplaceProblem::run ()\r
+{\r
+ GridGenerator::hyper_cube (triangulation, -1, 1);\r
+\r
+ for (unsigned int c=0; c<5; ++c)\r
+ {\r
+ triangulation.refine_global (1);\r
+ setup_system ();\r
+ assemble_system ();\r
+ solve ();\r
+ }\r
+\r
+ // finialise output\r
+ output_table.write_text (std::cout);\r
+ deallog << std::endl;\r
+}\r
+\r
+\r
+int main (int argc, char **argv)\r
+{\r
+ try\r
+ {\r
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);\r
+ {\r
+ LaplaceProblem problem;\r
+ problem.run ();\r
+ deallog << "OK" << std::endl;\r
+ }\r
+ }\r
+\r
+ catch (std::exception &exc)\r
+ {\r
+ std::cerr << std::endl << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ std::cerr << "Exception on processing: " << std::endl\r
+ << exc.what() << std::endl\r
+ << "Aborting!" << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+\r
+ return 1;\r
+ }\r
+ catch (...)\r
+ {\r
+ std::cerr << std::endl << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ std::cerr << "Unknown exception!" << std::endl\r
+ << "Aborting!" << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ return 1;\r
+ }\r
+\r
+ return 0;\r
+}\r