--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 1999 - 2021 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.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Wolfgang Bangerth, University of Heidelberg, 1999
+ */
+
+
+// Like the step-04 test, but output data via DataOutFaces
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out_faces.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+
+template <int dim>
+class Step4
+{
+public:
+ Step4();
+ void
+ run();
+
+private:
+ void
+ make_grid();
+ void
+ setup_system();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ output_results() const;
+
+ Triangulation<dim> triangulation;
+ FE_SimplexP<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ MappingFE<dim> mapping;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+};
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double return_value = 0.0;
+ for (unsigned int i = 0; i < dim; ++i)
+ return_value += 4.0 * std::pow(p(i), 4.0);
+
+ return return_value;
+}
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+ : fe(1)
+ , dof_handler(triangulation)
+ , mapping(fe)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+ Triangulation<dim, dim> temp_tria;
+ GridGenerator::subdivided_hyper_cube(temp_tria, 1, -1., 1., false);
+ GridGenerator::convert_hypercube_to_simplex_mesh<dim, dim>(temp_tria,
+ triangulation);
+
+ deallog << " Number of active cells: " << triangulation.n_active_cells()
+ << std::endl
+ << " Total number of cells: " << triangulation.n_cells()
+ << std::endl;
+}
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+ dof_handler.distribute_dofs(fe);
+
+ deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void
+Step4<dim>::assemble_system()
+{
+ QGaussSimplex<dim> quadrature_formula(fe.degree + 1);
+
+ RightHandSide<dim> right_hand_side;
+
+ const UpdateFlags flag = update_JxW_values | update_values |
+ update_gradients | update_quadrature_points;
+
+ FEValues<dim> fe_values(mapping, fe, quadrature_formula, flag);
+
+
+ const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ for (const unsigned int q_index : fe_values.quadrature_point_indices())
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ for (const unsigned int j : fe_values.dof_indices())
+ cell_matrix(i, j) +=
+ (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+ fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+ fe_values.JxW(q_index)); // dx
+
+ const auto x_q = fe_values.quadrature_point(q_index);
+ cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+ right_hand_side.value(x_q) * // f(x_q)
+ fe_values.JxW(q_index)); // dx
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ for (const unsigned int j : fe_values.dof_indices())
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i, j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values(
+ mapping, dof_handler, 0, BoundaryValues<dim>(), boundary_values);
+ MatrixTools::apply_boundary_values(boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::solve()
+{
+ SolverControl solver_control(1000, 1e-12);
+ SolverCG<Vector<double>> solver(solver_control);
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ deallog << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence." << std::endl;
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::output_results() const
+{
+ DataOutFaces<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(mapping, 0);
+
+ data_out.write_vtk(deallog.get_file_stream());
+}
+
+
+template <int dim>
+void
+Step4<dim>::run()
+{
+ deallog << "Solving problem in " << dim << " space dimensions." << std::endl;
+
+ make_grid();
+ setup_system();
+ assemble_system();
+ solve();
+ output_results();
+}
+
+
+
+int
+main()
+{
+ initlog();
+ deallog.depth_console(2);
+ {
+ Step4<2> laplace_problem_2d;
+ laplace_problem_2d.run();
+ }
+
+ {
+ Step4<3> laplace_problem_3d;
+ laplace_problem_3d.run();
+ }
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Solving problem in 2 space dimensions.
+DEAL:: Number of active cells: 8
+DEAL:: Total number of cells: 8
+DEAL:: Number of degrees of freedom: 9
+DEAL:cg::Starting value 4.18267
+DEAL:cg::Convergence step 1 value 0.00000
+DEAL:: 1 CG iterations needed to obtain convergence.
+# vtk DataFile Version 3.0
+#This file was generated
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 16 double
+-1.00000 -1.00000 0
+0.00000 -1.00000 0
+-1.00000 0.00000 0
+-1.00000 -1.00000 0
+1.00000 -1.00000 0
+1.00000 0.00000 0
+0.00000 -1.00000 0
+1.00000 -1.00000 0
+-1.00000 1.00000 0
+-1.00000 0.00000 0
+0.00000 1.00000 0
+-1.00000 1.00000 0
+1.00000 1.00000 0
+0.00000 1.00000 0
+1.00000 0.00000 0
+1.00000 1.00000 0
+
+CELLS 8 24
+2 0 1
+2 2 3
+2 4 5
+2 6 7
+2 8 9
+2 10 11
+2 12 13
+2 14 15
+
+CELL_TYPES 8
+ 3 3 3 3 3 3 3 3
+POINT_DATA 16
+SCALARS solution double 1
+LOOKUP_TABLE default
+2.00000 1.00000 1.00000 2.00000 2.00000 1.00000 1.00000 2.00000 2.00000 1.00000 1.00000 2.00000 2.00000 1.00000 1.00000 2.00000
+DEAL::Solving problem in 3 space dimensions.
+DEAL:: Number of active cells: 24
+DEAL:: Total number of cells: 24
+DEAL:: Number of degrees of freedom: 14
+DEAL:cg::Starting value 0.00000
+DEAL:cg::Convergence step 0 value 0.00000
+DEAL:: 0 CG iterations needed to obtain convergence.
+# vtk DataFile Version 3.0
+#This file was generated
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 96 double
+-1.00000 -1.00000 -1.00000
+1.00000 -1.00000 -1.00000
+0.00000 0.00000 -1.00000
+2.00000 0.00000 -1.00000
+1.00000 -1.00000 -1.00000
+-1.00000 -1.00000 -1.00000
+0.00000 -1.00000 0.00000
+-2.00000 -1.00000 0.00000
+-1.00000 1.00000 -1.00000
+1.00000 1.00000 -1.00000
+0.00000 1.00000 0.00000
+2.00000 1.00000 0.00000
+1.00000 1.00000 -1.00000
+-1.00000 1.00000 -1.00000
+0.00000 0.00000 -1.00000
+-2.00000 0.00000 -1.00000
+1.00000 1.00000 1.00000
+-1.00000 1.00000 1.00000
+0.00000 1.00000 0.00000
+-2.00000 1.00000 0.00000
+-1.00000 1.00000 1.00000
+1.00000 1.00000 1.00000
+0.00000 0.00000 1.00000
+2.00000 0.00000 1.00000
+1.00000 -1.00000 1.00000
+-1.00000 -1.00000 1.00000
+0.00000 0.00000 1.00000
+-2.00000 0.00000 1.00000
+-1.00000 -1.00000 1.00000
+1.00000 -1.00000 1.00000
+0.00000 -1.00000 0.00000
+2.00000 -1.00000 0.00000
+-1.00000 -1.00000 -1.00000
+-1.00000 1.00000 -1.00000
+-1.00000 0.00000 0.00000
+-1.00000 2.00000 0.00000
+-1.00000 1.00000 -1.00000
+-1.00000 -1.00000 -1.00000
+0.00000 0.00000 -1.00000
+0.00000 -2.00000 -1.00000
+-1.00000 -1.00000 1.00000
+-1.00000 1.00000 1.00000
+0.00000 0.00000 1.00000
+0.00000 2.00000 1.00000
+-1.00000 1.00000 1.00000
+-1.00000 -1.00000 1.00000
+-1.00000 0.00000 0.00000
+-1.00000 -2.00000 0.00000
+1.00000 -1.00000 1.00000
+0.00000 0.00000 1.00000
+1.00000 1.00000 1.00000
+0.00000 2.00000 1.00000
+1.00000 -1.00000 1.00000
+1.00000 1.00000 1.00000
+1.00000 0.00000 0.00000
+1.00000 2.00000 0.00000
+1.00000 -1.00000 -1.00000
+1.00000 0.00000 0.00000
+1.00000 1.00000 -1.00000
+1.00000 2.00000 0.00000
+1.00000 -1.00000 -1.00000
+1.00000 1.00000 -1.00000
+0.00000 0.00000 -1.00000
+0.00000 2.00000 -1.00000
+-1.00000 -1.00000 -1.00000
+-1.00000 0.00000 0.00000
+-1.00000 -1.00000 1.00000
+-1.00000 0.00000 2.00000
+-1.00000 -1.00000 -1.00000
+-1.00000 -1.00000 1.00000
+0.00000 -1.00000 0.00000
+0.00000 -1.00000 2.00000
+1.00000 -1.00000 -1.00000
+1.00000 -1.00000 1.00000
+1.00000 0.00000 0.00000
+1.00000 0.00000 2.00000
+1.00000 -1.00000 1.00000
+1.00000 -1.00000 -1.00000
+0.00000 -1.00000 0.00000
+0.00000 -1.00000 -2.00000
+1.00000 1.00000 -1.00000
+1.00000 1.00000 1.00000
+0.00000 1.00000 0.00000
+0.00000 1.00000 2.00000
+1.00000 1.00000 1.00000
+1.00000 1.00000 -1.00000
+1.00000 0.00000 0.00000
+1.00000 0.00000 -2.00000
+-1.00000 1.00000 -1.00000
+-1.00000 1.00000 1.00000
+-1.00000 0.00000 0.00000
+-1.00000 0.00000 2.00000
+-1.00000 1.00000 1.00000
+-1.00000 1.00000 -1.00000
+0.00000 1.00000 0.00000
+0.00000 1.00000 -2.00000
+
+CELLS 24 120
+4 0 1 3 2
+4 4 5 7 6
+4 8 9 11 10
+4 12 13 15 14
+4 16 17 19 18
+4 20 21 23 22
+4 24 25 27 26
+4 28 29 31 30
+4 32 33 35 34
+4 36 37 39 38
+4 40 41 43 42
+4 44 45 47 46
+4 48 49 51 50
+4 52 53 55 54
+4 56 57 59 58
+4 60 61 63 62
+4 64 65 67 66
+4 68 69 71 70
+4 72 73 75 74
+4 76 77 79 78
+4 80 81 83 82
+4 84 85 87 86
+4 88 89 91 90
+4 92 93 95 94
+
+CELL_TYPES 24
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+POINT_DATA 96
+SCALARS solution double 1
+LOOKUP_TABLE default
+3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000