std::bitset<width> mask;
- using ScalarNumber = typename Number::value_type;
-
for (unsigned int v = 0; v < width; ++v)
{
ScalarNumber a = 0.0;
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().
* 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
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().
+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;
+ }
}
}
}
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;
}
};
- 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)
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;
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();
+ }
}
{
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(
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
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);
{
// 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)
{
{
// 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)
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff