]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add mesh_loop() test: step-11 4961/head
authorTimo Heister <timo.heister@gmail.com>
Fri, 25 Aug 2017 18:56:57 +0000 (14:56 -0400)
committerTimo Heister <timo.heister@gmail.com>
Fri, 25 Aug 2017 19:18:18 +0000 (15:18 -0400)
tests/meshworker/step-11-mesh_loop.cc [new file with mode: 0644]
tests/meshworker/step-11-mesh_loop.output [new file with mode: 0644]

diff --git a/tests/meshworker/step-11-mesh_loop.cc b/tests/meshworker/step-11-mesh_loop.cc
new file mode 100644 (file)
index 0000000..0a4dd4e
--- /dev/null
@@ -0,0 +1,404 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// step-11 but rewritten using mesh_loop()
+
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/table_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.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_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/meshworker/mesh_loop.h>
+
+// Just this one is new: it declares a class
+// DynamicSparsityPattern, which we will use and explain
+// further down below.
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+// We will make use of the std::find algorithm of the C++ standard library, so
+// we have to include the following file for its declaration:
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+// The last step is as in all previous programs:
+namespace Step11
+{
+  using namespace dealii;
+
+  // Then we declare a class which represents the solution of a Laplace
+  // problem. As this example program is based on step-5, the class looks
+  // rather the same, with the sole structural difference that the functions
+  // <code>assemble_system</code> now calls <code>solve</code> itself, and is
+  // thus called <code>assemble_and_solve</code>, and that the output function
+  // was dropped since the solution function is so boring that it is not worth
+  // being viewed.
+  //
+  // The only other noteworthy change is that the constructor takes a value
+  // representing the polynomial degree of the mapping to be used later on,
+  // and that it has another member variable representing exactly this
+  // mapping. In general, this variable will occur in real applications at the
+  // same places where the finite element is declared or used.
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem (const unsigned int mapping_degree);
+    void run ();
+
+  private:
+    void setup_system ();
+    void assemble_and_solve ();
+    void solve ();
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+    MappingQ<dim>        mapping;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    ConstraintMatrix     constraints;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    TableHandler         output_table;
+  };
+
+
+
+  // Construct such an object, by initializing the variables. Here, we use
+  // linear finite elements (the argument to the <code>fe</code> variable
+  // denotes the polynomial degree), and mappings of given order. Print to
+  // screen what we are about to do.
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem (const unsigned int mapping_degree) :
+    fe (1),
+    dof_handler (triangulation),
+    mapping (mapping_degree)
+  {
+    deallog << "Using mapping with degree " << mapping_degree << ":\n"
+            << "============================"
+            << std::endl;
+  }
+
+
+
+  // The first task is to set up the variables for this problem. This includes
+  // generating a valid <code>DoFHandler</code> object, as well as the
+  // sparsity patterns for the matrix, and the object representing the
+  // constraints that the mean value of the degrees of freedom on the boundary
+  // be zero.
+  template <int dim>
+  void LaplaceProblem<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+
+    std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
+    DoFTools::extract_boundary_dofs (dof_handler,
+                                     ComponentMask(),
+                                     boundary_dofs);
+
+    const unsigned int first_boundary_dof
+      = std::distance (boundary_dofs.begin(),
+                       std::find (boundary_dofs.begin(),
+                                  boundary_dofs.end(),
+                                  true));
+
+    // Then generate a constraints object with just this one constraint. First
+    // clear all previous content (which might reside there from the previous
+    // computation on a once coarser grid), then add this one line
+    // constraining the <code>first_boundary_dof</code> to the sum of other
+    // boundary DoFs each with weight -1. Finally, close the constraints
+    // object, i.e. do some internal bookkeeping on it for faster processing
+    // of what is to come later:
+    constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+    constraints.add_line (first_boundary_dof);
+    for (unsigned int i=first_boundary_dof+1; i<dof_handler.n_dofs(); ++i)
+      if (boundary_dofs[i] == true)
+        constraints.add_entry (first_boundary_dof,
+                               i, -1);
+    constraints.close ();
+
+    DynamicSparsityPattern dsp (dof_handler.n_dofs(),
+                                dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints);
+
+
+    sparsity_pattern.copy_from (dsp);
+    system_matrix.reinit (sparsity_pattern);
+  }
+
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData (const Mapping<dim> &mapping,
+                 const FiniteElement<dim> &fe,
+                 const unsigned int quadrature_degree)
+      :
+      fe_values (mapping,
+                 fe,
+                 QGauss<dim>(quadrature_degree),
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values),
+      fe_face_values (mapping,
+                      fe,
+                      QGauss<dim-1>(quadrature_degree),
+                      update_values     | update_quadrature_points   |
+                      update_JxW_values | update_normal_vectors)
+    {}
+
+    ScratchData (const ScratchData<dim> &scratch_data)
+      :
+      fe_values (scratch_data.fe_values.get_mapping(),
+                 scratch_data.fe_values.get_fe(),
+                 scratch_data.fe_values.get_quadrature(),
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values),
+      fe_face_values (scratch_data.fe_face_values.get_mapping(),
+                      scratch_data.fe_face_values.get_fe(),
+                      scratch_data.fe_face_values.get_quadrature(),
+                      update_values     | update_quadrature_points   |
+                      update_JxW_values | update_normal_vectors)
+    {}
+
+    FEValues<dim>     fe_values;
+    FEFaceValues<dim> fe_face_values;
+
+  };
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+  };
+
+  // The next function then assembles the linear system of equations, solves
+  // it, and evaluates the solution.
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_and_solve ()
+  {
+
+
+    typedef decltype(dof_handler.begin_active()) Iterator;
+
+    auto cell_worker = [] (const Iterator &cell, ScratchData<dim> &scratch_data, CopyData &copy_data)
+    {
+      const unsigned int dofs_per_cell   = scratch_data.fe_values.get_fe().dofs_per_cell;
+      const unsigned int n_q_points      = scratch_data.fe_values.get_quadrature().size();
+
+      copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      copy_data.cell_rhs.reinit (dofs_per_cell);
+
+      copy_data.local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices (copy_data.local_dof_indices);
+
+      scratch_data.fe_values.reinit (cell);
+
+      std::vector<double>         rhs_values (n_q_points);
+      ConstantFunction<dim> right_hand_side (-2.0);
+      right_hand_side.value_list (scratch_data.fe_values.get_quadrature_points(),
+                                  rhs_values);
+
+      // ... and assemble the local contributions to the system matrix and
+      // right hand side as also discussed above:
+      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)
+              copy_data.cell_matrix(i,j) += (
+                                              scratch_data.fe_values.shape_grad (i, q_point) *
+                                              scratch_data.fe_values.shape_grad (j, q_point) *
+                                              scratch_data.fe_values.JxW (q_point));
+
+            copy_data.cell_rhs(i) += (scratch_data.fe_values.shape_value(i,q_point)
+                                      * rhs_values[q_point] *
+                                      scratch_data.fe_values.JxW (q_point));
+          }
+    };
+
+    auto boundary_worker = [] (const Iterator &cell, const unsigned int &f, ScratchData<dim> &scratch_data, CopyData &copy_data)
+    {
+      const unsigned int dofs_per_cell   = scratch_data.fe_values.get_fe().dofs_per_cell;
+      const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
+
+      std::vector<double>         face_boundary_values (n_face_q_points);
+      ConstantFunction<dim> boundary_values (1.0);
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        if (cell->face(face)->at_boundary())
+          {
+            scratch_data.fe_face_values.reinit (cell, face);
+            boundary_values.value_list (scratch_data.fe_face_values.get_quadrature_points(),
+                                        face_boundary_values);
+
+            for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                copy_data.cell_rhs(i) -= (face_boundary_values[q_point] *
+                                          scratch_data.fe_face_values.shape_value(i,q_point) *
+                                          scratch_data.fe_face_values.JxW(q_point));
+          }
+    };
+
+    auto copier = [&](const CopyData &c)
+    {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs
+                                            );
+    };
+
+    const unsigned int gauss_degree
+      = std::max (static_cast<unsigned int>(std::ceil(1.*(mapping.get_degree()+1)/2)),
+                  2U);
+
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          ScratchData<dim>(mapping, fe, gauss_degree),
+                          CopyData(),
+                          MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::cells_first,
+                          boundary_worker);
+
+
+    solve ();
+
+
+    // Finally, evaluate what we got as solution. As stated in the
+    // introduction, we are interested in the H1 semi-norm of the
+    // solution. Here, as well, we have a function in the library that does
+    // this, although in a slightly non-obvious way: the
+    // <code>VectorTools::integrate_difference</code> function integrates the
+    // norm of the difference between a finite element function and a
+    // continuous function. If we therefore want the norm of a finite element
+    // field, we just put the continuous function to zero. Note that this
+    // function, just as so many other ones in the library as well, has at
+    // least two versions, one which takes a mapping as argument (which we
+    // make us of here), and the one which we have used in previous examples
+    // which implicitly uses <code>MappingQ1</code>.  Also note that we take a
+    // quadrature formula of one degree higher, in order to avoid
+    // superconvergence effects where the solution happens to be especially
+    // close to the exact solution at certain points (we don't know whether
+    // this might be the case here, but there are cases known of this, and we
+    // just want to make sure):
+    Vector<float> norm_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (mapping, dof_handler,
+                                       solution,
+                                       ZeroFunction<dim>(),
+                                       norm_per_cell,
+                                       QGauss<dim>(gauss_degree+1),
+                                       VectorTools::H1_seminorm);
+    // Then, the function just called returns its results as a vector of
+    // values each of which denotes the norm on one cell. To get the global
+    // norm, we do the following:
+    const double norm = VectorTools::compute_global_error(triangulation,
+                                                          norm_per_cell,
+                                                          VectorTools::H1_seminorm);
+
+    // Last task -- generate output:
+    output_table.add_value ("cells", triangulation.n_active_cells());
+    output_table.add_value ("|u|_1", norm);
+    output_table.add_value ("error", std::fabs(norm-std::sqrt(3.14159265358/2)));
+  }
+
+
+
+  // The following function solving the linear system of equations is copied
+  // from step-5 and is explained there in some detail:
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    SolverControl           solver_control (1000, 1e-12, false, false);
+    SolverCG<>              cg (solver_control);
+
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    cg.solve (system_matrix, solution, system_rhs,
+              preconditioner);
+
+    constraints.distribute (solution);
+  }
+
+
+
+  // Finally the main function controlling the different steps to be
+  // performed.
+  template <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    GridGenerator::hyper_ball (triangulation);
+    static const SphericalManifold<dim> boundary;
+    triangulation.set_all_manifold_ids_on_boundary(0);
+    triangulation.set_manifold (0, boundary);
+
+    for (unsigned int cycle=0; cycle<5; ++cycle)
+      {
+        if (cycle>0)
+          triangulation.refine_global(1);
+
+        setup_system ();
+        assemble_and_solve ();
+      };
+
+    // After all the data is generated, write a table of results to the
+    // screen:
+    output_table.set_precision("|u|_1", 6);
+    output_table.set_precision("error", 6);
+    output_table.write_text (deallog.get_file_stream());
+    deallog << std::endl;
+  }
+}
+
+int main ()
+{
+  initlog();
+
+  for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree)
+    Step11::LaplaceProblem<2>(mapping_degree).run ();
+
+  return 0;
+}
diff --git a/tests/meshworker/step-11-mesh_loop.output b/tests/meshworker/step-11-mesh_loop.output
new file mode 100644 (file)
index 0000000..b1c061c
--- /dev/null
@@ -0,0 +1,28 @@
+
+DEAL::Using mapping with degree 1:
+============================
+cells  |u|_1    error   
+5     0.680402 0.572912 
+20    1.088141 0.165173 
+80    1.210142 0.043172 
+320   1.242375 0.010939 
+1280  1.250569 0.002745 
+DEAL::
+DEAL::Using mapping with degree 2:
+============================
+cells  |u|_1    error   
+5     1.050963 0.202351 
+20    1.203014 0.050300 
+80    1.241120 0.012194 
+320   1.250319 0.002996 
+1280  1.252572 0.000742 
+DEAL::
+DEAL::Using mapping with degree 3:
+============================
+cells  |u|_1    error   
+5     1.069381 0.183934 
+20    1.206606 0.046708 
+80    1.241611 0.011703 
+320   1.250382 0.002932 
+1280  1.252580 0.000734 
+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.