]> https://gitweb.dealii.org/ - dealii.git/commitdiff
TensorProductMatrixSymmetricSumCollection: allow to disable compression 14331/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 2 Oct 2022 14:38:59 +0000 (16:38 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 10 Dec 2022 06:32:46 +0000 (07:32 +0100)
include/deal.II/lac/tensor_product_matrix.h
include/deal.II/numerics/tensor_product_matrix_creator.h
tests/lac/tensor_product_matrix.h [new file with mode: 0644]
tests/lac/tensor_product_matrix_08.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_08.with_lapack=true.output [new file with mode: 0644]

index a59f68dae011bc3f07ce0b52f50553a0c8d38897..27f27169abad93b7eab0b3ecbda89634fedac379 100644 (file)
@@ -301,8 +301,6 @@ namespace internal
 
         std::bitset<width> mask;
 
-        using ScalarNumber = typename Number::value_type;
-
         for (unsigned int v = 0; v < width; ++v)
           {
             ScalarNumber a = 0.0;
@@ -370,6 +368,28 @@ class TensorProductMatrixSymmetricSumCollection
   using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
 
 public:
+  /**
+   * Struct to configure TensorProductMatrixSymmetricSumCollection.
+   */
+  struct AdditionalData
+  {
+    /**
+     * Constructor.
+     */
+    AdditionalData(const bool compress_matrices = true);
+
+    /**
+     * Try to compress internal matrices. Default: true.
+     */
+    bool compress_matrices;
+  };
+
+  /**
+   * Consturctor.
+   */
+  TensorProductMatrixSymmetricSumCollection(
+    const AdditionalData &additional_data = AdditionalData());
+
   /**
    * Allocate memory. The parameter @p specifies the maximum value
    * of the index used in invert() and apply_inverse().
@@ -382,10 +402,9 @@ public:
    * stiffness matrices @p Ks to the stored data, looking out for
    * compression possibilities.
    */
+  template <typename T>
   void
-  insert(const unsigned int                       index,
-         const std::array<Table<2, Number>, dim> &Ms,
-         const std::array<Table<2, Number>, dim> &Ks);
+  insert(const unsigned int index, const T &Ms, const T &Ks);
 
   /**
    * Finalize setup. This function computes, e.g., the
@@ -421,6 +440,17 @@ public:
   storage_size() const;
 
 private:
+  /**
+   * Try to compress matrices.
+   */
+  const bool compress_matrices;
+
+  /**
+   * Container used to collect 1D matrices if no compression is
+   * requested. The memory is freed during finalize().
+   */
+  std::vector<MatrixPairType> mass_and_derivative_matrices;
+
   /**
    * Container used during setup to determine the unique 1D
    * matrices. The memory is freed during finalize().
@@ -1041,35 +1071,69 @@ TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
 
 
 
+template <int dim, typename Number, int n_rows_1d>
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
+  AdditionalData::AdditionalData(const bool compress_matrices)
+  : compress_matrices(compress_matrices)
+{}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
+  TensorProductMatrixSymmetricSumCollection(
+    const AdditionalData &additional_data)
+  : compress_matrices(additional_data.compress_matrices)
+{}
+
+
+
 template <int dim, typename Number, int n_rows_1d>
 void
 TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::reserve(
   const unsigned int size)
 {
-  indices.assign(size * dim, numbers::invalid_unsigned_int);
+  if (compress_matrices == false)
+    mass_and_derivative_matrices.resize(size * dim);
+  else
+    indices.assign(size * dim, numbers::invalid_unsigned_int);
 }
 
 
 
 template <int dim, typename Number, int n_rows_1d>
+template <typename T>
 void
 TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::insert(
-  const unsigned int                       index,
-  const std::array<Table<2, Number>, dim> &Ms,
-  const std::array<Table<2, Number>, dim> &Ks)
+  const unsigned int index,
+  const T &          Ms_in,
+  const T &          Ks_in)
 {
+  const auto Ms =
+    internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
+  const auto Ks =
+    internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
+
   for (unsigned int d = 0; d < dim; ++d)
     {
       const MatrixPairType matrix(Ms[d], Ks[d]);
 
-      const auto ptr = cache.find(matrix);
-
-      if (ptr != cache.end())
-        indices[index * dim + d] = ptr->second;
+      if (compress_matrices == false)
+        {
+          mass_and_derivative_matrices[index * dim + d] = matrix;
+        }
       else
         {
-          indices[index * dim + d] = cache.size();
-          cache[matrix]            = cache.size();
+          const auto ptr = cache.find(matrix);
+
+          if (ptr != cache.end())
+            indices[index * dim + d] = ptr->second;
+          else
+            {
+              const auto size          = cache.size();
+              indices[index * dim + d] = size;
+              cache[matrix]            = size;
+            }
         }
     }
 }
@@ -1080,9 +1144,6 @@ template <int dim, typename Number, int n_rows_1d>
 void
 TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
 {
-  this->vector_ptr.resize(cache.size() + 1);
-  this->matrix_ptr.resize(cache.size() + 1);
-
   const auto store = [&](const unsigned int    index,
                          const MatrixPairType &M_and_K) {
     std::array<Table<2, Number>, 1> mass_matrix;
@@ -1114,8 +1175,47 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
       }
   };
 
-  if (cache.size() == indices.size())
+  if (compress_matrices == false)
+    {
+      // case 1) no compression requested
+
+      AssertDimension(cache.size(), 0);
+      AssertDimension(indices.size(), 0);
+
+      this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
+      this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
+
+      for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
+        {
+          const auto &M = mass_and_derivative_matrices[i].first;
+
+          this->vector_ptr[i + 1] = M.n_rows();
+          this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
+        }
+
+      for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
+        {
+          this->vector_ptr[i + 1] += this->vector_ptr[i];
+          this->matrix_ptr[i + 1] += this->matrix_ptr[i];
+        }
+
+      this->mass_matrices.resize_fast(matrix_ptr.back());
+      this->derivative_matrices.resize_fast(matrix_ptr.back());
+      this->eigenvectors.resize_fast(matrix_ptr.back());
+      this->eigenvalues.resize_fast(vector_ptr.back());
+
+      for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
+        store(i, mass_and_derivative_matrices[i]);
+
+      mass_and_derivative_matrices.clear();
+    }
+  else if (cache.size() == indices.size())
     {
+      // case 2) compression requested but none possible
+
+      this->vector_ptr.resize(cache.size() + 1);
+      this->matrix_ptr.resize(cache.size() + 1);
+
       std::map<unsigned int, MatrixPairType> inverted_cache;
 
       for (const auto &i : cache)
@@ -1135,18 +1235,24 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
           this->matrix_ptr[i + 1] += this->matrix_ptr[i];
         }
 
-      this->mass_matrices.resize(matrix_ptr.back());
-      this->derivative_matrices.resize(matrix_ptr.back());
-      this->eigenvectors.resize(matrix_ptr.back());
-      this->eigenvalues.resize(vector_ptr.back());
+      this->mass_matrices.resize_fast(matrix_ptr.back());
+      this->derivative_matrices.resize_fast(matrix_ptr.back());
+      this->eigenvectors.resize_fast(matrix_ptr.back());
+      this->eigenvalues.resize_fast(vector_ptr.back());
 
       for (unsigned int i = 0; i < indices.size(); ++i)
         store(i, inverted_cache[indices[i]]);
 
       indices.clear();
+      cache.clear();
     }
   else
     {
+      // case 2) compression requested but none possible
+
+      this->vector_ptr.resize(cache.size() + 1);
+      this->matrix_ptr.resize(cache.size() + 1);
+
       for (const auto &i : cache)
         {
           const auto &M = i.first.first;
@@ -1161,16 +1267,16 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
           this->matrix_ptr[i + 1] += this->matrix_ptr[i];
         }
 
-      this->mass_matrices.resize(matrix_ptr.back());
-      this->derivative_matrices.resize(matrix_ptr.back());
-      this->eigenvectors.resize(matrix_ptr.back());
-      this->eigenvalues.resize(vector_ptr.back());
+      this->mass_matrices.resize_fast(matrix_ptr.back());
+      this->derivative_matrices.resize_fast(matrix_ptr.back());
+      this->eigenvectors.resize_fast(matrix_ptr.back());
+      this->eigenvalues.resize_fast(vector_ptr.back());
 
       for (const auto &i : cache)
         store(i.second, i.first);
-    }
 
-  cache.clear();
+      cache.clear();
+    }
 }
 
 
index 1242e64b80027c93ec07474d15837452f31afd2c..184851b04ff93a03b76abf30b673a7a24ea3bba9 100644 (file)
@@ -99,6 +99,19 @@ namespace TensorProductMatrixCreator
 {
   namespace internal
   {
+    template <typename Number>
+    void
+    clear_row_and_column(const unsigned int  n_dofs_1D_with_overlap,
+                         const unsigned int  n,
+                         FullMatrix<Number> &matrix)
+    {
+      for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
+        {
+          matrix[i][n] = 0.0;
+          matrix[n][i] = 0.0;
+        }
+    }
+
     template <typename Number>
     std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>
     create_reference_mass_and_stiffness_matrices(
@@ -150,7 +163,8 @@ namespace TensorProductMatrixCreator
                  fe_values.JxW(q_index));
             }
 
-      return {mass_matrix_reference, derivative_matrix_reference, false};
+      return std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>{
+        mass_matrix_reference, derivative_matrix_reference, false};
     }
   } // namespace internal
 
@@ -191,14 +205,6 @@ namespace TensorProductMatrixCreator
     std::array<FullMatrix<Number>, dim> Ms;
     std::array<FullMatrix<Number>, dim> Ks;
 
-    const auto clear_row_and_column = [&](const unsigned int n, auto &matrix) {
-      for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
-        {
-          matrix[i][n] = 0.0;
-          matrix[n][i] = 0.0;
-        }
-    };
-
     for (unsigned int d = 0; d < dim; ++d)
       {
         Ms[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
@@ -235,8 +241,12 @@ namespace TensorProductMatrixCreator
               {
                 // left DBC
                 const unsigned i0 = n_overlap - 1;
-                clear_row_and_column(i0, Ms[d]);
-                clear_row_and_column(i0, Ks[d]);
+                internal::clear_row_and_column(n_dofs_1D_with_overlap,
+                                               i0,
+                                               Ms[d]);
+                internal::clear_row_and_column(n_dofs_1D_with_overlap,
+                                               i0,
+                                               Ks[d]);
               }
             else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann)
               {
@@ -268,8 +278,12 @@ namespace TensorProductMatrixCreator
               {
                 // right DBC
                 const unsigned i0 = n_overlap + n_dofs_1D - 2;
-                clear_row_and_column(i0, Ms[d]);
-                clear_row_and_column(i0, Ks[d]);
+                internal::clear_row_and_column(n_dofs_1D_with_overlap,
+                                               i0,
+                                               Ms[d]);
+                internal::clear_row_and_column(n_dofs_1D_with_overlap,
+                                               i0,
+                                               Ks[d]);
               }
             else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann)
               {
diff --git a/tests/lac/tensor_product_matrix.h b/tests/lac/tensor_product_matrix.h
new file mode 100644 (file)
index 0000000..9e311eb
--- /dev/null
@@ -0,0 +1,133 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#pragma once
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
+namespace dealii
+{
+  namespace GridTools
+  {
+    // Compute harminic extent of all non-artificial cells.
+    template <int dim>
+    std::vector<std::array<double, dim>>
+    compute_harmonic_cell_extent(const Mapping<dim> &       mapping,
+                                 const Triangulation<dim> & triangulation,
+                                 const Quadrature<dim - 1> &quadrature)
+    {
+      std::vector<std::array<double, dim>> result(
+        triangulation.n_active_cells());
+
+      FE_Nothing<dim>   fe_nothing;
+      FEFaceValues<dim> fe_face_values_0(mapping,
+                                         fe_nothing,
+                                         quadrature,
+                                         update_quadrature_points);
+      FEFaceValues<dim> fe_face_values_1(mapping,
+                                         fe_nothing,
+                                         quadrature,
+                                         update_quadrature_points);
+
+      for (const auto &cell : triangulation.active_cell_iterators())
+        if (cell->is_artificial() == false)
+          {
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                fe_face_values_0.reinit(cell, 2 * d + 0);
+                fe_face_values_1.reinit(cell, 2 * d + 1);
+
+                double extent = 0.0;
+
+                for (unsigned int q = 0; q < quadrature.size(); ++q)
+                  extent += fe_face_values_0.quadrature_point(q).distance(
+                              fe_face_values_1.quadrature_point(q)) *
+                            quadrature.weight(q);
+
+                result[cell->active_cell_index()][d] = extent;
+              }
+          }
+
+      return result;
+    }
+
+    // Compute harmonic extent of each locally owned cell including of each
+    // of its neighbors. If there is no neigbor, its extent is zero.
+    template <int dim>
+    std::vector<dealii::ndarray<double, dim, 3>>
+    compute_harmonic_patch_extent(const Mapping<dim> &       mapping,
+                                  const Triangulation<dim> & triangulation,
+                                  const Quadrature<dim - 1> &quadrature)
+    {
+      // 1) compute extent of each non-artificial cell
+      const auto harmonic_cell_extents =
+        GridTools::compute_harmonic_cell_extent(mapping,
+                                                triangulation,
+                                                quadrature);
+
+      // 2) accumulate for each face the normal extent for the
+      // neigboring cell(s); here we also consider periodicies
+      std::vector<double> face_extent(triangulation.n_faces(), 0.0);
+
+      for (const auto &cell : triangulation.active_cell_iterators())
+        if (cell->is_artificial() == false)
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              const auto extent =
+                harmonic_cell_extents[cell->active_cell_index()][d];
+
+              const auto add_extent_to_faces = [&](const unsigned int face_no) {
+                face_extent[cell->face(face_no)->index()] += extent;
+
+                if (cell->has_periodic_neighbor(face_no) &&
+                    (cell->periodic_neighbor(face_no)->is_artificial() ==
+                     false))
+                  face_extent[cell->periodic_neighbor(face_no)
+                                ->face(cell->periodic_neighbor_face_no(face_no))
+                                ->index()] += extent;
+              };
+
+              add_extent_to_faces(2 * d + 0); // face 0
+              add_extent_to_faces(2 * d + 1); // face 1
+            }
+
+      // 3) collect cell extent including those of the neighboring
+      // cells, which corresponds to the difference of extent of the
+      // current cell and the face extent
+      std::vector<dealii::ndarray<double, dim, 3>> result(
+        triangulation.n_active_cells());
+
+      for (const auto &cell : triangulation.active_cell_iterators())
+        if (cell->is_locally_owned())
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              const auto cell_extent =
+                harmonic_cell_extents[cell->active_cell_index()][d];
+
+              const auto index = cell->active_cell_index();
+
+              result[index][d][0] =
+                face_extent[cell->face(2 * d + 0)->index()] - cell_extent;
+              result[index][d][1] = cell_extent;
+              result[index][d][2] =
+                face_extent[cell->face(2 * d + 1)->index()] - cell_extent;
+            }
+
+      return result;
+    }
+
+  } // namespace GridTools
+} // namespace dealii
diff --git a/tests/lac/tensor_product_matrix_08.cc b/tests/lac/tensor_product_matrix_08.cc
new file mode 100644 (file)
index 0000000..508869c
--- /dev/null
@@ -0,0 +1,217 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test compression of TensorProductMatrixSymmetricSumCollection.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q_cache.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/tensor_product_matrix_creator.h>
+
+#include "../tests.h"
+
+#include "./tensor_product_matrix.h"
+
+template <int dim, typename Number>
+void
+do_test_mesh(const Mapping<dim> &mapping, const Triangulation<dim> &tria)
+{
+  using FDM = TensorProductMatrixSymmetricSumCollection<dim, Number>;
+
+  const unsigned int fe_degree = 3;
+  const unsigned int n_overlap = 1;
+
+  FE_Q<dim>       fe(fe_degree);
+  FE_Q<1>         fe_1D(fe_degree);
+  QGauss<1>       quadrature_1D(fe_degree + 1);
+  QGauss<dim - 1> quadrature_face(fe_degree + 1);
+
+  const auto harmonic_patch_extent =
+    GridTools::compute_harmonic_patch_extent(mapping, tria, quadrature_face);
+
+  FDM collection_0(typename FDM::AdditionalData(true));
+  FDM collection_1(typename FDM::AdditionalData(false));
+
+  collection_0.reserve(tria.n_active_cells());
+  collection_1.reserve(tria.n_active_cells());
+
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      const auto &patch_extent =
+        harmonic_patch_extent[cell->active_cell_index()];
+
+      const auto M_and_K = TensorProductMatrixCreator::
+        create_laplace_tensor_product_matrix<dim, Number>(
+          cell,
+          std::set<unsigned int>{0},
+          std::set<unsigned int>{},
+          fe_1D,
+          quadrature_1D,
+          patch_extent,
+          n_overlap);
+
+      collection_0.insert(cell->active_cell_index(),
+                          M_and_K.first,
+                          M_and_K.second);
+      collection_1.insert(cell->active_cell_index(),
+                          M_and_K.first,
+                          M_and_K.second);
+    }
+
+  collection_0.finalize();
+  collection_1.finalize();
+
+  deallog << "Storage sizes: " << collection_0.storage_size() << " "
+          << collection_1.storage_size() << std::endl;
+
+  Vector<Number>        src(fe.n_dofs_per_cell());
+  Vector<Number>        dst(fe.n_dofs_per_cell());
+  AlignedVector<Number> tmp;
+  FullMatrix<Number>    matrix_0(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
+  FullMatrix<Number>    matrix_1(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
+
+  for (unsigned int cell = 0; cell < tria.n_active_cells(); ++cell)
+    {
+      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+        {
+          for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+            src[j] = i == j;
+
+          collection_0.apply_inverse(cell,
+                                     make_array_view(dst),
+                                     make_array_view(src),
+                                     tmp);
+          for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+            matrix_0[j][i] = dst[j];
+
+          collection_1.apply_inverse(cell,
+                                     make_array_view(dst),
+                                     make_array_view(src),
+                                     tmp);
+          for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+            matrix_1[j][i] = dst[j];
+        }
+
+
+      FloatingPointComparator<Number> comp(1e-5, false);
+
+      Assert((comp.compare(matrix_0, matrix_1) ==
+              FloatingPointComparator<Number>::ComparisonResult::equal),
+             ExcInternalError());
+    }
+
+  deallog << "OK!" << std::endl;
+  deallog << std::endl;
+}
+
+
+template <typename Number>
+void
+do_test()
+{
+  const int dim = 2;
+
+  MappingQ1<dim> mapping;
+
+  const std::vector<std::tuple<unsigned int, bool, bool>> settings = {
+    std::tuple<unsigned int, bool, bool>{0, false, false}, // hyper-cube
+    std::tuple<unsigned int, bool, bool>{0, true, false}, // hyper-cube with PBC
+    std::tuple<unsigned int, bool, bool>{1, false, false}, // hyper-rectangle
+    std::tuple<unsigned int, bool, bool>{1,
+                                         true,
+                                         false}, // hyper-rectangle with PBC
+    std::tuple<unsigned int, bool, bool>{0, false, true} // deformed hyper-cube
+  };
+
+  for (const auto setting : settings)
+    {
+      Triangulation<dim> tria;
+
+      const auto do_pbc = std::get<1>(setting);
+
+      switch (std::get<0>(setting))
+        {
+          case 0:
+            GridGenerator::hyper_cube(tria, 0.0, 1.0, do_pbc);
+            break;
+          case 1:
+            GridGenerator::subdivided_hyper_rectangle(
+              tria, {1, 2}, {0, 0}, {1.0, 1.0}, do_pbc);
+            break;
+          default:
+            AssertThrow(false, ExcNotImplemented());
+        }
+
+      if (do_pbc)
+        {
+          std::vector<GridTools::PeriodicFacePair<
+            typename Triangulation<dim>::cell_iterator>>
+            periodic_faces;
+          for (unsigned int d = 0; d < dim; ++d)
+            GridTools::collect_periodic_faces(
+              tria, 2 * d, 2 * d + 1, d, periodic_faces);
+          tria.add_periodicity(periodic_faces);
+        }
+
+      tria.refine_global(3);
+
+      MappingQCache<dim> mapping(3);
+
+      mapping.initialize(
+        MappingQ1<dim>(),
+        tria,
+        [&](const typename Triangulation<dim>::cell_iterator &,
+            const Point<dim> &point) -> Point<dim> {
+          Point<dim> result;
+
+          if (std::get<2>(setting) == false) // Cartesian mesh
+            return result;
+
+          for (unsigned int d = 0; d < dim; ++d)
+            result[d] = std::pow(point[d], d + 0.9) *
+                        std::pow(point[(d + 1) % dim], d + 1.1);
+
+          return result;
+        },
+        true);
+
+
+      do_test_mesh<dim, Number>(mapping, tria);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  do_test<double>();
+  do_test<float>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_08.with_lapack=true.output b/tests/lac/tensor_product_matrix_08.with_lapack=true.output
new file mode 100644 (file)
index 0000000..d7f63ad
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Storage sizes: 3 128
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 1 128
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 6 256
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 2 256
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 128 128
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 3 128
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 1 128
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 6 256
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 2 256
+DEAL::OK!
+DEAL::
+DEAL::Storage sizes: 128 128
+DEAL::OK!
+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.