]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert step-7 to a simplex test 11517/head
authorSebastian Proell <proell@lnm.mw.tum.de>
Wed, 13 Jan 2021 18:01:32 +0000 (19:01 +0100)
committerSebastian Proell <proell@lnm.mw.tum.de>
Wed, 13 Jan 2021 18:46:12 +0000 (19:46 +0100)
tests/simplex/step-07.cc [new file with mode: 0644]
tests/simplex/step-07.with_simplex_support=on.output [new file with mode: 0644]

diff --git a/tests/simplex/step-07.cc b/tests/simplex/step-07.cc
new file mode 100644 (file)
index 0000000..4cbbd6a
--- /dev/null
@@ -0,0 +1,805 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2000 - 2020 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+
+// Step-07 on a simplex mesh. Following incompatible modifications had to be
+// made:
+//  - The function GridGenerator::hyper_cube() is only working for hypercube
+//    meshes. Use this function to create a temporary mesh and convert it to
+//    simplex mesh.
+//  - adaptive_refinement cases are not run
+// Other notable adaptations:
+//  - reduced refinement level to keep number of cells manageable
+
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/smartpointer.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.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/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.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/affine_constraints.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.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <array>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#define USE_SIMPLEX
+
+namespace Step7
+{
+  using namespace dealii;
+
+
+  template <int dim>
+  class SolutionBase
+  {
+  protected:
+    static const std::array<Point<dim>, 3> source_centers;
+    static const double                    width;
+  };
+
+  template <>
+  const std::array<Point<1>, 3> SolutionBase<1>::source_centers = {
+    {Point<1>(-1.0 / 3.0), Point<1>(0.0), Point<1>(+1.0 / 3.0)}};
+
+  template <>
+  const std::array<Point<2>, 3> SolutionBase<2>::source_centers = {
+    {Point<2>(-0.5, +0.5), Point<2>(-0.5, -0.5), Point<2>(+0.5, -0.5)}};
+
+  template <int dim>
+  const double SolutionBase<dim>::width = 1. / 8.;
+
+
+
+  template <int dim>
+  class Solution : public Function<dim>, protected SolutionBase<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+  };
+
+  template <int dim>
+  double
+  Solution<dim>::value(const Point<dim> &p, const unsigned int) const
+  {
+    double return_value = 0;
+    for (const auto &center : this->source_centers)
+      {
+        const Tensor<1, dim> x_minus_xi = p - center;
+        return_value +=
+          std::exp(-x_minus_xi.norm_square() / (this->width * this->width));
+      }
+
+    return return_value;
+  }
+
+  template <int dim>
+  Tensor<1, dim>
+  Solution<dim>::gradient(const Point<dim> &p, const unsigned int) const
+  {
+    Tensor<1, dim> return_value;
+
+    for (const auto &center : this->source_centers)
+      {
+        const Tensor<1, dim> x_minus_xi = p - center;
+
+        // For the gradient, note that its direction is along (x-x_i), so we
+        // add up multiples of this distance vector, where the factor is given
+        // by the exponentials.
+        return_value +=
+          (-2. / (this->width * this->width) *
+           std::exp(-x_minus_xi.norm_square() / (this->width * this->width)) *
+           x_minus_xi);
+      }
+
+    return return_value;
+  }
+
+
+
+  template <int dim>
+  class RightHandSide : public Function<dim>, protected SolutionBase<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) const
+  {
+    double return_value = 0;
+    for (const auto &center : this->source_centers)
+      {
+        const Tensor<1, dim> x_minus_xi = p - center;
+
+        // The first contribution is the Laplacian:
+        return_value +=
+          ((2. * dim -
+            4. * x_minus_xi.norm_square() / (this->width * this->width)) /
+           (this->width * this->width) *
+           std::exp(-x_minus_xi.norm_square() / (this->width * this->width)));
+        // And the second is the solution itself:
+        return_value +=
+          std::exp(-x_minus_xi.norm_square() / (this->width * this->width));
+      }
+
+    return return_value;
+  }
+
+
+
+  template <int dim>
+  class HelmholtzProblem
+  {
+  public:
+    enum RefinementMode
+    {
+      global_refinement,
+      adaptive_refinement
+    };
+
+    HelmholtzProblem(const FiniteElement<dim> &fe,
+                     const RefinementMode      refinement_mode);
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+    void
+    process_solution(const unsigned int cycle);
+
+    Triangulation<dim> triangulation;
+    DoFHandler<dim>    dof_handler;
+
+    SmartPointer<const FiniteElement<dim>> fe;
+
+    AffineConstraints<double> hanging_node_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+    const RefinementMode refinement_mode;
+
+    ConvergenceTable convergence_table;
+  };
+
+  template <int dim>
+  HelmholtzProblem<dim>::HelmholtzProblem(const FiniteElement<dim> &fe,
+                                          const RefinementMode refinement_mode)
+    : dof_handler(triangulation)
+    , fe(&fe)
+    , refinement_mode(refinement_mode)
+  {}
+
+  template <int dim>
+  void
+  HelmholtzProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(*fe);
+    DoFRenumbering::Cuthill_McKee(dof_handler);
+
+    hanging_node_constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler,
+                                            hanging_node_constraints);
+    hanging_node_constraints.close();
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+    hanging_node_constraints.condense(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
+  HelmholtzProblem<dim>::assemble_system()
+  {
+#ifdef USE_SIMPLEX
+    Simplex::QGauss<dim>     quadrature_formula(fe->degree + 1);
+    Simplex::QGauss<dim - 1> face_quadrature_formula(fe->degree + 1);
+#else
+    QGauss<dim>     quadrature_formula(fe->degree + 1);
+    QGauss<dim - 1> face_quadrature_formula(fe->degree + 1);
+#endif
+
+    const unsigned int n_q_points      = quadrature_formula.size();
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+    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);
+
+#ifdef USE_SIMPLEX
+    MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+#else
+    MappingFE<dim>  mapping(FE_Q<dim>(1));
+#endif
+
+    FEValues<dim> fe_values(mapping,
+                            *fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    FEFaceValues<dim> fe_face_values(mapping,
+                                     *fe,
+                                     face_quadrature_formula,
+                                     update_values | update_quadrature_points |
+                                       update_normal_vectors |
+                                       update_JxW_values);
+
+    const RightHandSide<dim> right_hand_side;
+    std::vector<double>      rhs_values(n_q_points);
+
+    Solution<dim> exact_solution;
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        cell_matrix = 0.;
+        cell_rhs    = 0.;
+
+        fe_values.reinit(cell);
+
+        right_hand_side.value_list(fe_values.get_quadrature_points(),
+                                   rhs_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_matrix(i, j) +=
+                  ((fe_values.shape_grad(i, q_point) *     // grad phi_i(x_q)
+                      fe_values.shape_grad(j, q_point)     // grad phi_j(x_q)
+                    +                                      //
+                    fe_values.shape_value(i, q_point) *    // phi_i(x_q)
+                      fe_values.shape_value(j, q_point)) * // phi_j(x_q)
+                   fe_values.JxW(q_point));                // dx
+
+
+              cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
+                              rhs_values[q_point] *               // f(x_q)
+                              fe_values.JxW(q_point));            // dx
+            }
+
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && (face->boundary_id() == 1))
+            {
+              fe_face_values.reinit(cell, face);
+
+
+              for (unsigned int q_point = 0; q_point < n_face_q_points;
+                   ++q_point)
+                {
+                  const double neumann_value =
+                    (exact_solution.gradient(
+                       fe_face_values.quadrature_point(q_point)) *
+                     fe_face_values.normal_vector(q_point));
+
+                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                    cell_rhs(i) +=
+                      (neumann_value *                          // g(x_q)
+                       fe_face_values.shape_value(i, q_point) * // phi_i(x_q)
+                       fe_face_values.JxW(q_point));            // dx
+                }
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              system_matrix.add(local_dof_indices[i],
+                                local_dof_indices[j],
+                                cell_matrix(i, j));
+
+            system_rhs(local_dof_indices[i]) += cell_rhs(i);
+          }
+      }
+
+    hanging_node_constraints.condense(system_matrix);
+    hanging_node_constraints.condense(system_rhs);
+
+    std::map<types::global_dof_index, double> boundary_values;
+    VectorTools::interpolate_boundary_values(
+      mapping, dof_handler, 0, Solution<dim>(), boundary_values);
+    MatrixTools::apply_boundary_values(boundary_values,
+                                       system_matrix,
+                                       solution,
+                                       system_rhs);
+  }
+
+  template <int dim>
+  void
+  HelmholtzProblem<dim>::solve()
+  {
+    SolverControl            solver_control(1000, 1e-12);
+    SolverCG<Vector<double>> cg(solver_control);
+
+    PreconditionSSOR<SparseMatrix<double>> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    cg.solve(system_matrix, solution, system_rhs, preconditioner);
+
+    hanging_node_constraints.distribute(solution);
+  }
+
+  template <int dim>
+  void
+  HelmholtzProblem<dim>::refine_grid()
+  {
+    switch (refinement_mode)
+      {
+        case global_refinement:
+          {
+            triangulation.refine_global(1);
+            break;
+          }
+
+        case adaptive_refinement:
+          {
+            Vector<float> estimated_error_per_cell(
+              triangulation.n_active_cells());
+
+            KellyErrorEstimator<dim>::estimate(
+              dof_handler,
+              QGauss<dim - 1>(fe->degree + 1),
+              std::map<types::boundary_id, const Function<dim> *>(),
+              solution,
+              estimated_error_per_cell);
+
+            GridRefinement::refine_and_coarsen_fixed_number(
+              triangulation, estimated_error_per_cell, 0.3, 0.03);
+
+            triangulation.execute_coarsening_and_refinement();
+
+            break;
+          }
+
+        default:
+          {
+            Assert(false, ExcNotImplemented());
+          }
+      }
+  }
+
+
+  template <int dim>
+  void
+  HelmholtzProblem<dim>::process_solution(const unsigned int cycle)
+  {
+    // error norms
+    Vector<float> difference_per_cell(triangulation.n_active_cells());
+
+#ifdef USE_SIMPLEX
+    MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+#else
+    MappingFE<dim>  mapping(FE_Q<dim>(1));
+#endif
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      QGauss<dim>(fe->degree + 1),
+                                      VectorTools::L2_norm);
+    const double L2_error =
+      VectorTools::compute_global_error(triangulation,
+                                        difference_per_cell,
+                                        VectorTools::L2_norm);
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      QGauss<dim>(fe->degree + 1),
+                                      VectorTools::H1_seminorm);
+    const double H1_error =
+      VectorTools::compute_global_error(triangulation,
+                                        difference_per_cell,
+                                        VectorTools::H1_seminorm);
+
+    const QTrapezoid<1>  q_trapez;
+    const QIterated<dim> q_iterated(q_trapez, fe->degree * 2 + 1);
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      q_iterated,
+                                      VectorTools::Linfty_norm);
+    const double Linfty_error =
+      VectorTools::compute_global_error(triangulation,
+                                        difference_per_cell,
+                                        VectorTools::Linfty_norm);
+
+    const unsigned int n_active_cells = triangulation.n_active_cells();
+    const unsigned int n_dofs         = dof_handler.n_dofs();
+
+    deallog << "Cycle " << cycle << ':' << std::endl
+            << "   Number of active cells:       " << n_active_cells
+            << std::endl
+            << "   Number of degrees of freedom: " << n_dofs << std::endl;
+
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", n_active_cells);
+    convergence_table.add_value("dofs", n_dofs);
+    convergence_table.add_value("L2", L2_error);
+    convergence_table.add_value("H1", H1_error);
+    convergence_table.add_value("Linfty", Linfty_error);
+  }
+
+  template <int dim>
+  void
+  HelmholtzProblem<dim>::run()
+  {
+    // note: reduced refinement to keep problem size manageable
+    const unsigned int n_cycles =
+      (refinement_mode == global_refinement) ? 3 : 9;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        if (cycle == 0)
+          {
+#ifdef USE_SIMPLEX
+            Triangulation<dim> temp;
+            GridGenerator::hyper_cube(temp, -1., 1.);
+            GridGenerator::convert_hypercube_to_simplex_mesh(temp,
+                                                             triangulation);
+#else
+            GridGenerator::hyper_cube(triangulation, -1., 1.);
+#endif
+            triangulation.refine_global(3);
+
+            for (const auto &cell : triangulation.cell_iterators())
+              for (const auto &face : cell->face_iterators())
+                {
+                  const auto center = face->center();
+                  if ((std::fabs(center(0) - (-1.0)) < 1e-12) ||
+                      (std::fabs(center(1) - (-1.0)) < 1e-12))
+                    face->set_boundary_id(1);
+                }
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        assemble_system();
+        solve();
+
+        process_solution(cycle);
+      }
+
+    std::string vtk_filename;
+    switch (refinement_mode)
+      {
+        case global_refinement:
+          vtk_filename = "solution-global";
+          break;
+        case adaptive_refinement:
+          vtk_filename = "solution-adaptive";
+          break;
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    switch (fe->degree)
+      {
+        case 1:
+          vtk_filename += "-q1";
+          break;
+        case 2:
+          vtk_filename += "-q2";
+          break;
+
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    vtk_filename += ".vtk";
+    std::ofstream output(vtk_filename);
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "solution");
+
+#ifdef USE_SIMPLEX
+    MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+#else
+    MappingFE<dim> mapping(FE_Q<dim>(1));
+#endif
+
+    data_out.build_patches(mapping, fe->degree);
+    // data_out.write_vtk(output);
+
+    convergence_table.set_precision("L2", 3);
+    convergence_table.set_precision("H1", 3);
+    convergence_table.set_precision("Linfty", 3);
+
+    convergence_table.set_scientific("L2", true);
+    convergence_table.set_scientific("H1", true);
+    convergence_table.set_scientific("Linfty", true);
+
+    convergence_table.set_tex_caption("cells", "\\# cells");
+    convergence_table.set_tex_caption("dofs", "\\# dofs");
+    convergence_table.set_tex_caption("L2", "$L^2$-error");
+    convergence_table.set_tex_caption("H1", "$H^1$-error");
+    convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error");
+
+    convergence_table.set_tex_format("cells", "r");
+    convergence_table.set_tex_format("dofs", "r");
+
+    deallog << std::endl;
+    convergence_table.write_text(deallog.get_file_stream());
+
+    std::string error_filename = "error";
+    switch (refinement_mode)
+      {
+        case global_refinement:
+          error_filename += "-global";
+          break;
+        case adaptive_refinement:
+          error_filename += "-adaptive";
+          break;
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    switch (fe->degree)
+      {
+        case 1:
+          error_filename += "-q1";
+          break;
+        case 2:
+          error_filename += "-q2";
+          break;
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    error_filename += ".tex";
+    std::ofstream error_table_file(error_filename);
+
+    convergence_table.write_tex(deallog.get_file_stream());
+
+    if (refinement_mode == global_refinement)
+      {
+        convergence_table.add_column_to_supercolumn("cycle", "n cells");
+        convergence_table.add_column_to_supercolumn("cells", "n cells");
+
+        std::vector<std::string> new_order;
+        new_order.emplace_back("n cells");
+        new_order.emplace_back("H1");
+        new_order.emplace_back("L2");
+        convergence_table.set_column_order(new_order);
+
+        convergence_table.evaluate_convergence_rates(
+          "L2", ConvergenceTable::reduction_rate);
+        convergence_table.evaluate_convergence_rates(
+          "L2", ConvergenceTable::reduction_rate_log2);
+        convergence_table.evaluate_convergence_rates(
+          "H1", ConvergenceTable::reduction_rate);
+        convergence_table.evaluate_convergence_rates(
+          "H1", ConvergenceTable::reduction_rate_log2);
+
+        deallog << std::endl;
+        convergence_table.write_text(deallog.get_file_stream());
+
+        std::string conv_filename = "convergence";
+        switch (refinement_mode)
+          {
+            case global_refinement:
+              conv_filename += "-global";
+              break;
+            case adaptive_refinement:
+              conv_filename += "-adaptive";
+              break;
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+        switch (fe->degree)
+          {
+            case 1:
+              conv_filename += "-q1";
+              break;
+            case 2:
+              conv_filename += "-q2";
+              break;
+            default:
+              Assert(false, ExcNotImplemented());
+          }
+        conv_filename += ".tex";
+
+        std::ofstream table_file(conv_filename);
+
+        // convergence_table.write_tex(deallog.get_file_stream());
+      }
+  }
+} // namespace Step7
+
+int
+main()
+{
+  initlog();
+
+  const unsigned int dim = 2;
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step7;
+
+      {
+        deallog << "Solving with P1 elements, adaptive refinement" << std::endl
+                << "=============================================" << std::endl
+                << std::endl;
+
+#ifdef USE_SIMPLEX
+        Simplex::FE_P<dim> fe(1);
+#else
+        FE_Q<dim> fe(1);
+#endif
+        HelmholtzProblem<dim> helmholtz_problem_2d(
+          fe, HelmholtzProblem<dim>::adaptive_refinement);
+
+        // TODO needs hanging node constraints
+#ifndef USE_SIMPLEX
+        helmholtz_problem_2d.run();
+#endif
+
+        deallog << std::endl;
+      }
+
+      {
+        deallog << "Solving with P1 elements, global refinement" << std::endl
+                << "===========================================" << std::endl
+                << std::endl;
+
+#ifdef USE_SIMPLEX
+        Simplex::FE_P<dim> fe(1);
+#else
+        FE_Q<dim> fe(1);
+#endif
+        HelmholtzProblem<dim> helmholtz_problem_2d(
+          fe, HelmholtzProblem<dim>::global_refinement);
+
+        helmholtz_problem_2d.run();
+
+        deallog << std::endl;
+      }
+
+      {
+        deallog << "Solving with P2 elements, global refinement" << std::endl
+                << "===========================================" << std::endl
+                << std::endl;
+
+#ifdef USE_SIMPLEX
+        Simplex::FE_P<dim> fe(2);
+#else
+        FE_Q<dim> fe(2);
+#endif
+        HelmholtzProblem<dim> helmholtz_problem_2d(
+          fe, HelmholtzProblem<dim>::global_refinement);
+
+        helmholtz_problem_2d.run();
+
+        deallog << std::endl;
+      }
+      {
+        deallog << "Solving with P2 elements, adaptive refinement" << std::endl
+                << "===========================================" << std::endl
+                << std::endl;
+
+#ifdef USE_SIMPLEX
+        Simplex::FE_P<dim> fe(2);
+#else
+        FE_Q<dim> fe(2);
+#endif
+        HelmholtzProblem<dim> helmholtz_problem_2d(
+          fe, HelmholtzProblem<dim>::adaptive_refinement);
+
+        // TODO needs hanging node constraints
+#ifndef USE_SIMPLEX
+        helmholtz_problem_2d.run();
+#endif
+
+        deallog << std::endl;
+      }
+    }
+  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;
+    }
+
+  return 0;
+}
diff --git a/tests/simplex/step-07.with_simplex_support=on.output b/tests/simplex/step-07.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..4fe8ec1
--- /dev/null
@@ -0,0 +1,99 @@
+
+DEAL::Solving with P1 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::
+DEAL::Solving with P1 elements, global refinement
+DEAL::===========================================
+DEAL::
+DEAL:cg::Starting value 4.36155
+DEAL:cg::Convergence step 32 value 4.41596e-13
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 289
+DEAL:cg::Starting value 2.73911
+DEAL:cg::Convergence step 58 value 8.84507e-13
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       2048
+DEAL::   Number of degrees of freedom: 1089
+DEAL:cg::Starting value 1.49437
+DEAL:cg::Convergence step 113 value 6.96045e-13
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       8192
+DEAL::   Number of degrees of freedom: 4225
+DEAL::
+cycle cells dofs    L2        H1      Linfty   
+0     512   289  1.342e-01 3.074e+00 8.185e-01 
+1     2048  1089 4.203e-02 1.878e+00 3.749e-01 
+2     8192  4225 1.131e-02 9.927e-01 1.125e-01 
+\documentclass[10pt]{report}
+\usepackage{float}
+
+
+\begin{document}
+\begin{table}[H]
+\begin{center}
+\begin{tabular}{|c|r|r|c|c|c|} \hline
+cycle & \# cells & \# dofs & $L^2$-error & $H^1$-error & $L^\infty$-error\\ \hline
+0 & 512 & 289 & 1.342e-01 & 3.074e+00 & 8.185e-01\\ \hline
+1 & 2048 & 1089 & 4.203e-02 & 1.878e+00 & 3.749e-01\\ \hline
+2 & 8192 & 4225 & 1.131e-02 & 9.927e-01 & 1.125e-01\\ \hline
+\end{tabular}
+\end{center}
+\end{table}
+\end{document}
+DEAL::
+n cells         H1                  L2          
+0 512   3.074e+00 -    -    1.342e-01 -    -    
+1 2048  1.878e+00 1.64 0.71 4.203e-02 3.19 1.68 
+2 8192  9.927e-01 1.89 0.92 1.131e-02 3.72 1.89 
+DEAL::
+DEAL::Solving with P2 elements, global refinement
+DEAL::===========================================
+DEAL::
+DEAL:cg::Starting value 3.03300
+DEAL:cg::Convergence step 76 value 7.67913e-13
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 1089
+DEAL:cg::Starting value 1.71371
+DEAL:cg::Convergence step 150 value 8.04906e-13
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       2048
+DEAL::   Number of degrees of freedom: 4225
+DEAL:cg::Starting value 0.880749
+DEAL:cg::Convergence step 296 value 9.07150e-13
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       8192
+DEAL::   Number of degrees of freedom: 16641
+DEAL::
+cycle cells dofs     L2        H1      Linfty   
+0     512   1089  6.841e-02 1.940e+00 8.066e-01 
+1     2048  4225  1.537e-02 7.758e-01 2.471e-01 
+2     8192  16641 2.177e-03 2.177e-01 3.643e-02 
+\documentclass[10pt]{report}
+\usepackage{float}
+
+
+\begin{document}
+\begin{table}[H]
+\begin{center}
+\begin{tabular}{|c|r|r|c|c|c|} \hline
+cycle & \# cells & \# dofs & $L^2$-error & $H^1$-error & $L^\infty$-error\\ \hline
+0 & 512 & 1089 & 6.841e-02 & 1.940e+00 & 8.066e-01\\ \hline
+1 & 2048 & 4225 & 1.537e-02 & 7.758e-01 & 2.471e-01\\ \hline
+2 & 8192 & 16641 & 2.177e-03 & 2.177e-01 & 3.643e-02\\ \hline
+\end{tabular}
+\end{center}
+\end{table}
+\end{document}
+DEAL::
+n cells         H1                  L2          
+0 512   1.940e+00 -    -    6.841e-02 -    -    
+1 2048  7.758e-01 2.50 1.32 1.537e-02 4.45 2.15 
+2 8192  2.177e-01 3.56 1.83 2.177e-03 7.06 2.82 
+DEAL::
+DEAL::Solving with P2 elements, adaptive refinement
+DEAL::===========================================
+DEAL::
+DEAL::

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.