]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added FEInterfaceValues to ScratchData.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 23 Jun 2020 20:24:58 +0000 (22:24 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 23 Jun 2020 20:45:07 +0000 (22:45 +0200)
include/deal.II/fe/fe_interface_values.h
include/deal.II/meshworker/scratch_data.h
source/meshworker/scratch_data.cc
tests/meshworker/scratch_data_08.cc [new file with mode: 0644]
tests/meshworker/scratch_data_08.output [new file with mode: 0644]

index c2c20978fb86cf6fd1ea6492d30cddd39483d816..7412b0a6089ec47571b6d859c3503eace64b2dd6 100644 (file)
@@ -340,7 +340,7 @@ public:
    * If this is a boundary face (at_boundary() returns true), then
    * $\average{\nabla u}=\nabla u_{\text{cell0}}$.
    */
-  Tensor<1, dim>
+  Tensor<1, spacedim>
   average_gradient(const unsigned int interface_dof_index,
                    const unsigned int q_point,
                    const unsigned int component = 0) const;
@@ -355,7 +355,7 @@ public:
    * If this is a boundary face (at_boundary() returns true), then
    * $\average{\nabla^2 u}=\nabla^2 u_{\text{cell0}}$.
    */
-  Tensor<2, dim>
+  Tensor<2, spacedim>
   average_hessian(const unsigned int interface_dof_index,
                   const unsigned int q_point,
                   const unsigned int component = 0) const;
@@ -369,7 +369,7 @@ public:
    * If this is a boundary face (at_boundary() returns true), then
    * $\jump{\nabla u}=\nabla u_{\text{cell0}}$.
    */
-  Tensor<1, dim>
+  Tensor<1, spacedim>
   jump_gradient(const unsigned int interface_dof_index,
                 const unsigned int q_point,
                 const unsigned int component = 0) const;
@@ -384,7 +384,7 @@ public:
    * If this is a boundary face (at_boundary() returns true), then
    * $\jump{\nabla^2 u} = \nabla^2 u_{\text{cell0}}$.
    */
-  Tensor<2, dim>
+  Tensor<2, spacedim>
   jump_hessian(const unsigned int interface_dof_index,
                const unsigned int q_point,
                const unsigned int component = 0) const;
@@ -398,7 +398,7 @@ public:
    * If this is a boundary face (at_boundary() returns true), then
    * $\jump{\nabla^3 u} = \nabla^3 u_{\text{cell0}}$.
    */
-  Tensor<3, dim>
+  Tensor<3, spacedim>
   jump_3rd_derivative(const unsigned int interface_dof_index,
                       const unsigned int q_point,
                       const unsigned int component = 0) const;
@@ -423,35 +423,35 @@ private:
   /**
    * The FEFaceValues object for the current cell.
    */
-  FEFaceValues<dim> internal_fe_face_values;
+  FEFaceValues<dim, spacedim> internal_fe_face_values;
 
   /**
    * The FEFaceValues object for the current cell if the cell is refined.
    */
-  FESubfaceValues<dim> internal_fe_subface_values;
+  FESubfaceValues<dim, spacedim> internal_fe_subface_values;
 
   /**
    * The FEFaceValues object for the neighboring cell.
    */
-  FEFaceValues<dim> internal_fe_face_values_neighbor;
+  FEFaceValues<dim, spacedim> internal_fe_face_values_neighbor;
 
   /**
    * The FEFaceValues object for the neighboring cell if the cell is refined.
    */
-  FESubfaceValues<dim> internal_fe_subface_values_neighbor;
+  FESubfaceValues<dim, spacedim> internal_fe_subface_values_neighbor;
 
   /**
    * Pointer to internal_fe_face_values or internal_fe_subface_values,
    * respectively as determined in reinit().
    */
-  FEFaceValuesBase<dim> *fe_face_values;
+  FEFaceValuesBase<dim, spacedim> *fe_face_values;
 
   /**
    * Pointer to internal_fe_face_values_neighbor,
    * internal_fe_subface_values_neighbor, or nullptr, respectively
    * as determined in reinit().
    */
-  FEFaceValuesBase<dim> *fe_face_values_neighbor;
+  FEFaceValuesBase<dim, spacedim> *fe_face_values_neighbor;
 };
 
 
@@ -818,7 +818,7 @@ FEInterfaceValues<dim, spacedim>::average(
 
 
 template <int dim, int spacedim>
-Tensor<1, dim>
+Tensor<1, spacedim>
 FEInterfaceValues<dim, spacedim>::average_gradient(
   const unsigned int interface_dof_index,
   const unsigned int q_point,
@@ -831,7 +831,7 @@ FEInterfaceValues<dim, spacedim>::average_gradient(
                                                       q_point,
                                                       component);
 
-  Tensor<1, dim> value;
+  Tensor<1, spacedim> value;
 
   if (dof_pair[0] != numbers::invalid_unsigned_int)
     value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
@@ -848,7 +848,7 @@ FEInterfaceValues<dim, spacedim>::average_gradient(
 
 
 template <int dim, int spacedim>
-Tensor<2, dim>
+Tensor<2, spacedim>
 FEInterfaceValues<dim, spacedim>::average_hessian(
   const unsigned int interface_dof_index,
   const unsigned int q_point,
@@ -861,7 +861,7 @@ FEInterfaceValues<dim, spacedim>::average_hessian(
                                                          q_point,
                                                          component);
 
-  Tensor<2, dim> value;
+  Tensor<2, spacedim> value;
 
   if (dof_pair[0] != numbers::invalid_unsigned_int)
     value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
@@ -878,7 +878,7 @@ FEInterfaceValues<dim, spacedim>::average_hessian(
 
 
 template <int dim, int spacedim>
-Tensor<1, dim>
+Tensor<1, spacedim>
 FEInterfaceValues<dim, spacedim>::jump_gradient(
   const unsigned int interface_dof_index,
   const unsigned int q_point,
@@ -891,7 +891,7 @@ FEInterfaceValues<dim, spacedim>::jump_gradient(
                                                       q_point,
                                                       component);
 
-  Tensor<1, dim> value;
+  Tensor<1, spacedim> value;
 
   if (dof_pair[0] != numbers::invalid_unsigned_int)
     value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
@@ -908,7 +908,7 @@ FEInterfaceValues<dim, spacedim>::jump_gradient(
 
 
 template <int dim, int spacedim>
-Tensor<2, dim>
+Tensor<2, spacedim>
 FEInterfaceValues<dim, spacedim>::jump_hessian(
   const unsigned int interface_dof_index,
   const unsigned int q_point,
@@ -921,7 +921,7 @@ FEInterfaceValues<dim, spacedim>::jump_hessian(
                                                          q_point,
                                                          component);
 
-  Tensor<2, dim> value;
+  Tensor<2, spacedim> value;
 
   if (dof_pair[0] != numbers::invalid_unsigned_int)
     value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
@@ -937,7 +937,7 @@ FEInterfaceValues<dim, spacedim>::jump_hessian(
 
 
 template <int dim, int spacedim>
-Tensor<3, dim>
+Tensor<3, spacedim>
 FEInterfaceValues<dim, spacedim>::jump_3rd_derivative(
   const unsigned int interface_dof_index,
   const unsigned int q_point,
@@ -950,7 +950,7 @@ FEInterfaceValues<dim, spacedim>::jump_3rd_derivative(
                                                                 q_point,
                                                                 component);
 
-  Tensor<3, dim> value;
+  Tensor<3, spacedim> value;
 
   if (dof_pair[0] != numbers::invalid_unsigned_int)
     value +=
index e83f97cb6f0e3250965a7b370e6d53f44629c5b2..3a18fee1b0db053ba203e0f2f5c383c8d9c26d93 100644 (file)
@@ -22,6 +22,7 @@
 
 #include <deal.II/differentiation/ad.h>
 
+#include <deal.II/fe/fe_interface_values.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
 
@@ -54,9 +55,9 @@ namespace MeshWorker
    * function and with the MeshWorker::mesh_loop function().
    *
    * The ScratchData class has three main goals:
-   * - to create FEValues, FEFaceValues, and FESubfaceValues for the current
-   *   cell and for its neighbor cell *on demand* (only if they are
-   *   necessary for the algorithm provided by the user), and to provide a
+   * - to create FEValues, FEFaceValues, FESubfaceValues, and FEInterfaceValues
+   *   for the current cell and for its neighbor cell *on demand* (only if they
+   *   are necessary for the algorithm provided by the user), and to provide a
    *   uniform interface to access the FEValues objects when assembling cell,
    *   face, or subface contributions
    * - to store arbitrary data types (or references to arbitrary data types),
@@ -70,10 +71,10 @@ namespace MeshWorker
    *   values, gradients, etc., of already computed solution vectors.
    *
    * The methods in the section "Methods to work on current cell"
-   * initialize on demand internal FEValues,
-   * FEFaceValues, and FESubfaceValues objects on the current cell, allowing the
-   * use of this class as a single substitute for three different objects used
-   * to integrate and query finite element values on cells, faces, and subfaces.
+   * initialize on demand internal FEValues, FEFaceValues, FESubfaceValues, and
+   * FEInterfaceValues objects on the current cell, allowing the use of this
+   * class as a single substitute for four different objects used to integrate
+   * and query finite element values on cells, faces, and subfaces.
    *
    * Similarly, the methods in the section "Methods to work on neighbor cell"
    * initialize on demand
@@ -372,6 +373,29 @@ namespace MeshWorker
            const unsigned int face_no,
            const unsigned int subface_no);
 
+    /**
+     * Initialize the internal FEInterfaceValues with the given arguments, and
+     * return a reference to it.
+     *
+     * After calling this function, get_local_dof_indices(),
+     * get_quadrature_points(), get_normal_vectors(), and get_JxW_values() will
+     * be forwarded to the local FEInterfaceValues object. The methods
+     * get_current_fe_values() will return the FEValuesBase associated to the
+     * current cell, while get_neighbor_fe_values() will be associated with the
+     * neighbor cell. The method get_local_dof_indices() will return the
+     * same result of FEInterfaceValues::get_interface_dof_to_dof_indices(),
+     * while the get_neighbor_dof_indices() will return the local dof indices
+     * of the neighbor cell.
+     */
+    const FEInterfaceValues<dim, spacedim> &
+    reinit(const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
+           const unsigned int face_no,
+           const unsigned int sub_face_no,
+           const typename DoFHandler<dim, spacedim>::active_cell_iterator
+             &                cell_neighbor,
+           const unsigned int face_no_neighbor,
+           const unsigned int sub_face_no_neighbor);
+
     /**
      * Get the currently initialized FEValues.
      *
@@ -421,8 +445,8 @@ namespace MeshWorker
      * Initialize the internal neighbor FEValues to use the given @p cell, and
      * return a reference to it.
      *
-     * After calling this function, get_current_fe_values() will return the
-     * same object of this method, as an FEValuesBase reference.
+     * After calling this function, get_current_neighbor_fe_values() will return
+     * the same object of this method, as an FEValuesBase reference.
      */
     const FEValues<dim, spacedim> &
     reinit_neighbor(
@@ -432,8 +456,8 @@ namespace MeshWorker
      * Initialize the internal FEFaceValues to use the given @p face_no on the
      * given @p cell, and return a reference to it.
      *
-     * After calling this function, get_current_fe_values() will return the
-     * same object of this method, as an FEValuesBase reference.
+     * After calling this function, get_current_neighbor_fe_values() will return
+     * the same object of this method, as an FEValuesBase reference.
      */
     const FEFaceValues<dim, spacedim> &
     reinit_neighbor(
@@ -444,8 +468,8 @@ namespace MeshWorker
      * Initialize the internal FESubfaceValues to use the given @p subface_no,
      * on @p face_no, on the given @p cell, and return a reference to it.
      *
-     * After calling this function, get_current_fe_values() will return the
-     * same object of this method, as an FEValuesBase reference.
+     * After calling this function, get_current_neighbor_fe_values() will return
+     * the same object of this method, as an FEValuesBase reference.
      *
      * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function
      * that takes only the @p cell and the @p face_no is called.
@@ -862,6 +886,11 @@ namespace MeshWorker
      */
     std::unique_ptr<FESubfaceValues<dim, spacedim>> neighbor_fe_subface_values;
 
+    /**
+     * Interface values on facets.
+     */
+    std::unique_ptr<FEInterfaceValues<dim, spacedim>> interface_fe_values;
+
     /**
      * Dof indices on the current cell.
      */
@@ -886,13 +915,13 @@ namespace MeshWorker
      * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues
      * object on this cell.
      */
-    SmartPointer<FEValuesBase<dim, spacedim>> current_fe_values;
+    SmartPointer<const FEValuesBase<dim, spacedim>> current_fe_values;
 
     /**
      * A pointer to the last used FEValues/FEFaceValues, or FESubfaceValues
      * object on the neighbor cell.
      */
-    SmartPointer<FEValuesBase<dim, spacedim>> current_neighbor_fe_values;
+    SmartPointer<const FEValuesBase<dim, spacedim>> current_neighbor_fe_values;
   };
 
 #ifndef DOXYGEN
@@ -936,7 +965,7 @@ namespace MeshWorker
     const VectorType & input_vector,
     const Number       dummy)
   {
-    const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell;
+    const unsigned int n_dofs = local_dof_indices.size();
 
     const std::string name =
       get_unique_dofs_name(global_vector_name, n_dofs, dummy);
index f26a4a095c5c8234dbe1e1f1bb2b3871fcb20558..c2229cce32ee60eae73020b29be21d00d45a5b9d 100644 (file)
@@ -136,6 +136,7 @@ namespace MeshWorker
                                                             cell_update_flags);
 
     fe_values->reinit(cell);
+    local_dof_indices.resize(fe_values->dofs_per_cell);
     cell->get_dof_indices(local_dof_indices);
     current_fe_values = fe_values.get();
     return *fe_values;
@@ -154,6 +155,7 @@ namespace MeshWorker
         *mapping, *fe, face_quadrature, face_update_flags);
 
     fe_face_values->reinit(cell, face_no);
+    local_dof_indices.resize(fe->dofs_per_cell);
     cell->get_dof_indices(local_dof_indices);
     current_fe_values = fe_face_values.get();
     return *fe_face_values;
@@ -174,6 +176,7 @@ namespace MeshWorker
           fe_subface_values = std::make_unique<FESubfaceValues<dim, spacedim>>(
             *mapping, *fe, face_quadrature, face_update_flags);
         fe_subface_values->reinit(cell, face_no, subface_no);
+        local_dof_indices.resize(fe->dofs_per_cell);
         cell->get_dof_indices(local_dof_indices);
 
         current_fe_values = fe_subface_values.get();
@@ -185,6 +188,37 @@ namespace MeshWorker
 
 
 
+  template <int dim, int spacedim>
+  const FEInterfaceValues<dim, spacedim> &
+  ScratchData<dim, spacedim>::reinit(
+    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
+    const unsigned int                                              face_no,
+    const unsigned int                                              sub_face_no,
+    const typename DoFHandler<dim, spacedim>::active_cell_iterator
+      &                cell_neighbor,
+    const unsigned int face_no_neighbor,
+    const unsigned int sub_face_no_neighbor)
+  {
+    if (!interface_fe_values)
+      interface_fe_values = std::make_unique<FEInterfaceValues<dim, spacedim>>(
+        *mapping, *fe, face_quadrature, face_update_flags);
+    interface_fe_values->reinit(cell,
+                                face_no,
+                                sub_face_no,
+                                cell_neighbor,
+                                face_no_neighbor,
+                                sub_face_no_neighbor);
+
+    current_fe_values          = &interface_fe_values->get_fe_face_values(0);
+    current_neighbor_fe_values = &interface_fe_values->get_fe_face_values(1);
+
+    cell_neighbor->get_dof_indices(neighbor_dof_indices);
+    local_dof_indices = interface_fe_values->get_interface_dof_indices();
+    return *interface_fe_values;
+  }
+
+
+
   template <int dim, int spacedim>
   const FEValues<dim, spacedim> &
   ScratchData<dim, spacedim>::reinit_neighbor(
diff --git a/tests/meshworker/scratch_data_08.cc b/tests/meshworker/scratch_data_08.cc
new file mode 100644 (file)
index 0000000..6c2c526
--- /dev/null
@@ -0,0 +1,257 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2018 - 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// Solve Laplacian using SIPG + mesh_loop + ScratchData + FEInterfaceData
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/patterns.h>
+#include <deal.II/base/thread_management.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.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/sparse_direct.h>
+
+#include <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  Triangulation<dim, spacedim> tria;
+  FE_DGQ<dim, spacedim>        fe(1);
+  DoFHandler<dim, spacedim>    dh(tria);
+
+  Functions::ConstantFunction<spacedim> rhs_function(1);
+  Functions::ConstantFunction<spacedim> boundary_function(0);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(5);
+  tria.execute_coarsening_and_refinement();
+  dh.distribute_dofs(fe);
+
+
+  SparsityPattern sparsity;
+
+  {
+    DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dh, dsp);
+    sparsity.copy_from(dsp);
+  }
+
+  SparseMatrix<double> matrix;
+  matrix.reinit(sparsity);
+
+  Vector<double> solution(dh.n_dofs());
+  Vector<double> rhs(dh.n_dofs());
+
+  QGauss<dim>     quad(3);
+  QGauss<dim - 1> face_quad(3);
+
+  UpdateFlags cell_flags = update_values | update_gradients |
+                           update_quadrature_points | update_JxW_values;
+  UpdateFlags face_flags = update_values | update_gradients |
+                           update_quadrature_points |
+                           update_face_normal_vectors | update_JxW_values;
+
+  // Stabilization for SIPG
+  double gamma = 100;
+
+  using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+
+  auto cell = dh.begin_active();
+  auto endc = dh.end();
+
+  typedef decltype(cell) Iterator;
+
+  struct CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+  };
+
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<CopyDataFace>            face_data;
+
+    void
+    reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+    }
+  };
+
+  ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+  CopyData    copy;
+
+  auto cell_worker =
+    [&rhs_function](const Iterator &cell, ScratchData &s, CopyData &c) {
+      const auto &fev = s.reinit(cell);
+      const auto &JxW = s.get_JxW_values();
+      const auto &p   = s.get_quadrature_points();
+
+      c.reinit(cell, s.get_local_dof_indices().size());
+
+      for (unsigned int q = 0; q < p.size(); ++q)
+        for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < fev.dofs_per_cell; ++j)
+              {
+                c.cell_matrix(i, j) +=
+                  fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q];
+              }
+            c.cell_rhs(i) +=
+              fev.shape_value(i, q) * rhs_function.value(p[q]) * JxW[q];
+          }
+    };
+
+  auto boundary_worker = [gamma, &boundary_function](const Iterator &    cell,
+                                                     const unsigned int &f,
+                                                     ScratchData &       s,
+                                                     CopyData &          c) {
+    const auto &fev = s.reinit(cell, f);
+    const auto &JxW = s.get_JxW_values();
+    const auto &p   = s.get_quadrature_points();
+    const auto &n   = s.get_normal_vectors();
+
+    for (unsigned int q = 0; q < p.size(); ++q)
+      for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
+        {
+          for (unsigned int j = 0; j < fev.dofs_per_cell; ++j)
+            {
+              c.cell_matrix(i, j) +=
+                (-fev.shape_grad(i, q) * n[q] * fev.shape_value(j, q) +
+                 -fev.shape_grad(j, q) * n[q] * fev.shape_value(i, q) +
+                 gamma / cell->face(f)->diameter() * fev.shape_value(i, q) *
+                   fev.shape_value(j, q)) *
+                JxW[q];
+            }
+          c.cell_rhs(i) +=
+            ((gamma / cell->face(f)->diameter() * fev.shape_value(i, q) -
+              fev.shape_grad(i, q) * n[q]) *
+             boundary_function.value(p[q])) *
+            JxW[q];
+        }
+  };
+
+  auto face_worker = [gamma](const Iterator &    cell,
+                             const unsigned int &f,
+                             const unsigned int &sf,
+                             const Iterator &    ncell,
+                             const unsigned int &nf,
+                             const unsigned int &nsf,
+                             ScratchData &       s,
+                             CopyData &          c) {
+    const auto &fev = s.reinit(cell, f, sf, ncell, nf, nsf);
+    const auto &JxW = s.get_JxW_values();
+
+    const auto &p = s.get_quadrature_points();
+    const auto &n = s.get_normal_vectors();
+
+
+    c.face_data.emplace_back();
+    auto &copy_data_face             = c.face_data.back();
+    auto &face_matrix                = copy_data_face.cell_matrix;
+    copy_data_face.joint_dof_indices = fev.get_interface_dof_indices();
+    const auto n_dofs                = fev.n_current_interface_dofs();
+    face_matrix.reinit(n_dofs, n_dofs);
+
+    const double gh = gamma / cell->face(f)->diameter();
+
+    for (unsigned int q = 0; q < p.size(); ++q)
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        for (unsigned int j = 0; j < n_dofs; ++j)
+          {
+            face_matrix(i, j) +=
+              (-fev.jump_gradient(i, q) * n[q] * fev.average(j, q) -
+               fev.average(i, q) * fev.jump_gradient(j, q) * n[q] +
+               gh * fev.jump(i, q) * fev.jump(j, q)) *
+              JxW[q];
+          }
+  };
+
+  auto copier = [&constraints, &matrix, &rhs](const CopyData &c) {
+    constraints.distribute_local_to_global(
+      c.cell_matrix, c.cell_rhs, c.local_dof_indices, matrix, rhs);
+
+    for (auto &cdf : c.face_data)
+      constraints.distribute_local_to_global(cdf.cell_matrix,
+                                             cdf.joint_dof_indices,
+                                             matrix);
+  };
+
+
+
+  mesh_loop(cell,
+            endc,
+            cell_worker,
+            copier,
+            scratch,
+            copy,
+            assemble_own_cells | assemble_boundary_faces |
+              assemble_own_interior_faces_once,
+            boundary_worker,
+            face_worker);
+
+  SparseDirectUMFPACK inv;
+  inv.initialize(matrix);
+
+  inv.vmult(solution, rhs);
+
+  deallog << "Linfty norm of solution " << solution.linfty_norm() << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  test<2, 2>();
+}
diff --git a/tests/meshworker/scratch_data_08.output b/tests/meshworker/scratch_data_08.output
new file mode 100644 (file)
index 0000000..200ebcc
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Linfty norm of solution 0.0736360

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.