]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-65
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 9 May 2019 16:32:55 +0000 (18:32 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 13 Jun 2019 09:43:59 +0000 (11:43 +0200)
examples/step-65/CMakeLists.txt [new file with mode: 0644]
examples/step-65/doc/intro.dox [new file with mode: 0644]
examples/step-65/doc/results.dox [new file with mode: 0644]
examples/step-65/step-65.cc [new file with mode: 0644]

diff --git a/examples/step-65/CMakeLists.txt b/examples/step-65/CMakeLists.txt
new file mode 100644 (file)
index 0000000..bcedad9
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-65 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-65")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#    FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#    FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#    SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.2.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-65/doc/intro.dox b/examples/step-65/doc/intro.dox
new file mode 100644 (file)
index 0000000..26b84f2
--- /dev/null
@@ -0,0 +1,14 @@
+
+<br>
+
+<i>
+This program was contributed by Martin Kronbichler.
+</i>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+This tutorial program presents a way to handle expensive manifold classes,
+exemplified for the TransfiniteInterpolationManifold, a special manifold class
+that can blend a curved approximation on a boundary with straight-sided
+descriptions elsewhere.
diff --git a/examples/step-65/doc/results.dox b/examples/step-65/doc/results.dox
new file mode 100644 (file)
index 0000000..9738b1c
--- /dev/null
@@ -0,0 +1,61 @@
+
+<h1>Results</h1>
+
+<h3>Program output</h3>
+
+@code
+$ make run
+Scanning dependencies of target step-65
+[ 33%] Building CXX object CMakeFiles/step-65.dir/step-65.cc.o
+[ 66%] Linking CXX executable step-65
+[ 66%] Built target step-65
+[100%] Run step-65 with Release configuration
+
+====== Running with the basic MappingQGeneric class ====== 
+
+   Number of active cells:       6656
+   Number of degrees of freedom: 181609
+   Number of solver iterations:  285
+   L2 error vs exact solution:   8.99328e-08
+   H1 error vs exact solution:   6.45341e-06
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start    |      56.3s |            |
+|                                             |            |            |
+| Section                         | no. calls |  wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| Assemble linear system          |         1 |      5.39s |       9.6% |
+| Compute constraints             |         1 |     0.109s |      0.19% |
+| Compute error estimator         |         1 |      17.3s |        31% |
+| Compute error norms             |         1 |      9.53s |        17% |
+| Solve linear system             |         1 |      10.6s |        19% |
+| Write output                    |         2 |      10.3s |        18% |
++---------------------------------+-----------+------------+------------+
+
+====== Running with the optimized MappingQCache class ====== 
+
+   Memory consumption cache:     22.9976 MB
+   Number of active cells:       6656
+   Number of degrees of freedom: 181609
+   Number of solver iterations:  285
+   L2 error vs exact solution:   8.99328e-08
+   H1 error vs exact solution:   6.45341e-06
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start    |      19.7s |            |
+|                                             |            |            |
+| Section                         | no. calls |  wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| Assemble linear system          |         1 |     0.876s |       4.4% |
+| Compute constraints             |         1 |   0.00368s |         0% |
+| Compute error estimator         |         1 |     0.488s |       2.5% |
+| Compute error norms             |         1 |     0.513s |       2.6% |
+| Initialize mapping cache        |         1 |       5.2s |        26% |
+| Solve linear system             |         1 |      10.6s |        54% |
+| Write output                    |         2 |      1.86s |       9.5% |
++---------------------------------+-----------+------------+------------+
+
+[100%] Built target run
+@endcode
diff --git a/examples/step-65/step-65.cc b/examples/step-65/step-65.cc
new file mode 100644 (file)
index 0000000..44bdb5a
--- /dev/null
@@ -0,0 +1,418 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+// This tutorial program was contributed by Martin Kronbichler
+
+#include <deal.II/base/timer.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.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/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_cache.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+
+namespace step65
+{
+  using namespace dealii;
+
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      if (p.norm_square() < 0.25)
+        return p.norm_square();
+      else
+        return 0.1 * p.norm_square() + (0.25 - 0.025);
+    }
+
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> &p,
+             const unsigned int /*component*/ = 0) const override
+    {
+      if (p.norm_square() < 0.25)
+        return 2. * p;
+      else
+        return 0.2 * p;
+    }
+  };
+
+
+  template <int dim>
+  class Coefficient : public Function<dim>
+  {
+  public:
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      if (p.norm_square() < 0.25)
+        return 0.5;
+      else
+        return 5.;
+    }
+  };
+
+
+
+  template <int dim>
+  class PoissonProblem
+  {
+  public:
+    PoissonProblem();
+    void run();
+
+  private:
+    void create_grid();
+    void setup_system(const Mapping<dim> &mapping);
+    void assemble_system(const Mapping<dim> &mapping);
+    void solve();
+    void postprocess(const Mapping<dim> &mapping);
+    void output_results(const Mapping<dim> &mapping);
+
+    Triangulation<dim> triangulation;
+    FE_Q<dim>          fe;
+    DoFHandler<dim>    dof_handler;
+
+    AffineConstraints<double> constraints;
+    SparsityPattern           sparsity_pattern;
+    SparseMatrix<double>      system_matrix;
+    Vector<double>            solution;
+    Vector<double>            system_rhs;
+
+    TimerOutput timer;
+  };
+
+
+
+  template <int dim>
+  PoissonProblem<dim>::PoissonProblem()
+    : fe(3)
+    , dof_handler(triangulation)
+    , timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
+  {}
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::create_grid()
+  {
+    Triangulation<dim> tria_outer, tria_inner;
+    GridGenerator::hyper_shell(
+      tria_outer, Point<dim>(), 0.5, std::sqrt(dim), 2 * dim);
+    GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
+    GridGenerator::merge_triangulations(tria_inner, tria_outer, triangulation);
+    triangulation.reset_all_manifolds();
+    triangulation.set_all_manifold_ids(0);
+    for (const auto &cell : triangulation.cell_iterators())
+      {
+        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          {
+            bool face_at_sphere_boundary = true;
+            for (unsigned int v = 0;
+                 v < GeometryInfo<dim - 1>::vertices_per_cell;
+                 ++v)
+              if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
+                  1e-12)
+                face_at_sphere_boundary = false;
+            if (face_at_sphere_boundary)
+              cell->face(f)->set_all_manifold_ids(1);
+          }
+        if (cell->center().norm_square() < 0.25)
+          cell->set_material_id(1);
+        else
+          cell->set_material_id(0);
+      }
+    triangulation.set_manifold(1, SphericalManifold<dim>());
+    TransfiniteInterpolationManifold<dim> transfinite_manifold;
+    transfinite_manifold.initialize(triangulation);
+    triangulation.set_manifold(0, transfinite_manifold);
+
+    triangulation.refine_global(9 - 2 * dim);
+  }
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::setup_system(const Mapping<dim> &mapping)
+  {
+    dof_handler.distribute_dofs(fe);
+    std::cout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl;
+    std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl;
+
+    {
+      TimerOutput::Scope scope(timer, "Compute constraints");
+      constraints.clear();
+      DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+      VectorTools::interpolate_boundary_values(
+        mapping, dof_handler, 0, ExactSolution<dim>(), constraints);
+      constraints.close();
+    }
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+
+    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 PoissonProblem<dim>::assemble_system(const Mapping<dim> &mapping)
+  {
+    TimerOutput::Scope scope(timer, "Assemble linear system");
+    QGauss<dim>        quadrature_formula(fe.degree + 1);
+    FEValues<dim>      fe_values(mapping,
+                            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_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+    FullMatrix<double> partial_matrix(dofs_per_cell, dim * n_q_points);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        cell_rhs = 0.;
+        fe_values.reinit(cell);
+        for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
+          {
+            const double current_coefficient =
+              Coefficient<dim>().value(fe_values.quadrature_point(q_index));
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                for (unsigned int d = 0; d < dim; ++d)
+                  partial_matrix(i, q_index * dim + d) =
+                    std::sqrt(fe_values.JxW(q_index) * current_coefficient) *
+                    fe_values.shape_grad(i, q_index)[d];
+                cell_rhs(i) +=
+                  (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+                   (-dim) *                            // f(x_q)
+                   fe_values.JxW(q_index));            // dx
+              }
+          }
+        partial_matrix.mTmult(cell_matrix, partial_matrix);
+        cell->get_dof_indices(local_dof_indices);
+
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+      }
+  }
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::solve()
+  {
+    TimerOutput::Scope   scope(timer, "Solve linear system");
+    SolverControl        solver_control(1000, 1e-12);
+    SolverCG<>           solver(solver_control);
+    PreconditionJacobi<> preconditioner;
+    preconditioner.initialize(system_matrix);
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+    constraints.distribute(solution);
+    std::cout << "   Number of solver iterations:  "
+              << solver_control.last_step() << std::endl;
+  }
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::postprocess(const Mapping<dim> &mapping)
+  {
+    {
+      TimerOutput::Scope scope(timer, "Write output");
+      Timer              time;
+      DataOut<dim>       data_out;
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "solution");
+
+      DataOutBase::VtkFlags flags;
+      flags.write_higher_order_cells = true;
+      data_out.set_flags(flags);
+
+      Vector<double> material_ids(triangulation.n_active_cells());
+      for (const auto &cell : triangulation.active_cell_iterators())
+        material_ids[cell->active_cell_index()] = cell->material_id();
+      data_out.add_data_vector(material_ids, "material_ids");
+
+      data_out.build_patches(mapping,
+                             fe.degree,
+                             DataOut<dim>::curved_inner_cells);
+
+      std::ofstream file(
+        ("solution-" +
+         std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) +
+         ".vtu")
+          .c_str());
+      data_out.write_vtu(file);
+    }
+
+    {
+      TimerOutput::Scope scope(timer, "Compute error norms");
+      Vector<double>     norm_per_cell_p(triangulation.n_active_cells());
+
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        ExactSolution<dim>(),
+                                        norm_per_cell_p,
+                                        QGauss<dim>(fe.degree + 2),
+                                        VectorTools::L2_norm);
+      std::cout << "   L2 error vs exact solution:   "
+                << norm_per_cell_p.l2_norm() << std::endl;
+
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        ExactSolution<dim>(),
+                                        norm_per_cell_p,
+                                        QGauss<dim>(fe.degree + 2),
+                                        VectorTools::H1_norm);
+      std::cout << "   H1 error vs exact solution:   "
+                << norm_per_cell_p.l2_norm() << std::endl;
+    }
+
+    {
+      TimerOutput::Scope scope(timer, "Compute error estimator");
+      Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+      KellyErrorEstimator<dim>::estimate(
+        mapping,
+        dof_handler,
+        QGauss<dim - 1>(fe.degree + 1),
+        std::map<types::boundary_id, const Function<dim> *>(),
+        solution,
+        estimated_error_per_cell);
+    }
+  }
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::output_results(const Mapping<dim> &mapping)
+  {
+    TimerOutput::Scope scope(timer, "Write output");
+    Timer              time;
+    DataOut<dim>       data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "solution");
+
+    DataOutBase::VtkFlags flags;
+    flags.write_higher_order_cells = true;
+    data_out.set_flags(flags);
+
+    Vector<double> material_ids(triangulation.n_active_cells());
+    for (const auto &cell : triangulation.active_cell_iterators())
+      material_ids[cell->active_cell_index()] = cell->material_id();
+    data_out.add_data_vector(material_ids, "material_ids");
+
+    data_out.build_patches(mapping,
+                           fe.degree,
+                           DataOut<dim>::curved_inner_cells);
+
+    std::ofstream file(
+      ("solution-" +
+       std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) + ".vtu")
+        .c_str());
+    data_out.write_vtu(file);
+  }
+
+
+
+  template <int dim>
+  void PoissonProblem<dim>::run()
+  {
+    create_grid();
+    {
+      std::cout << std::endl
+                << "====== Running with the basic MappingQGeneric class ====== "
+                << std::endl
+                << std::endl;
+
+      MappingQGeneric<dim> mapping(fe.degree + 1);
+      setup_system(mapping);
+      assemble_system(mapping);
+      solve();
+      postprocess(mapping);
+      output_results(mapping);
+
+      timer.print_summary();
+    }
+
+    timer.reset();
+
+    {
+      std::cout
+        << "====== Running with the optimized MappingQCache class ====== "
+        << std::endl
+        << std::endl;
+
+      MappingQCache<dim> mapping(fe.degree + 1);
+      {
+        TimerOutput::Scope scope(timer, "Initialize mapping cache");
+        mapping.initialize(triangulation, MappingQGeneric<dim>(fe.degree + 1));
+      }
+      std::cout << "   Memory consumption cache:     "
+                << 1e-6 * mapping.memory_consumption() << " MB" << std::endl;
+
+      setup_system(mapping);
+      assemble_system(mapping);
+      solve();
+      postprocess(mapping);
+      output_results(mapping);
+
+      timer.print_summary();
+    }
+  }
+
+} // namespace step65
+
+
+int main()
+{
+  step65::PoissonProblem<3> test_program;
+  test_program.run();
+  return 0;
+}

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.