manifold.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2023 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.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_tria_manifold_h
#define dealii_tria_manifold_h
 
 
/*----------------------------   manifold.h     ---------------------------*/
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/subscriptor.h>
 
#include <deal.II/grid/tria.h>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declaration
#ifndef DOXYGEN
template <int, typename>
class Table;
#endif
 
namespace Manifolds
{
  template <typename MeshIteratorType>
  inline constexpr std::size_t
  n_default_points_per_cell()
  {
    // Note that in C++11 a constexpr function can only have a return
    // statement, so we cannot alias the structure dimension
    return GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::
             vertices_per_cell +
           GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::
             lines_per_cell +
           GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::
             quads_per_cell +
           GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::
             hexes_per_cell -
           1; // don't count the cell itself, just the bounding objects
  }
 
  template <typename MeshIteratorType>
  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
                       n_default_points_per_cell<MeshIteratorType>()>,
            std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
  get_default_points_and_weights(const MeshIteratorType &iterator,
                                 const bool with_interpolation = false);
} // namespace Manifolds
 
 
 
template <int dim, int spacedim = dim>
class Manifold : public Subscriptor
{
public:
  // explicitly check for sensible template arguments
  static_assert(dim <= spacedim,
                "The dimension <dim> of a Manifold must be less than or "
                "equal to the space dimension <spacedim> in which it lives.");
 
 
  using FaceVertexNormals =
    std::array<Tensor<1, spacedim>, GeometryInfo<dim>::vertices_per_face>;
 
 
  virtual ~Manifold() override = default;
 
  virtual std::unique_ptr<Manifold<dim, spacedim>>
  clone() const = 0;
 
 
  virtual Point<spacedim>
  get_intermediate_point(const Point<spacedim> &p1,
                         const Point<spacedim> &p2,
                         const double           w) const;
 
  virtual Point<spacedim>
  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
                const ArrayView<const double> &         weights) const;
 
 
  virtual void
  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                 const Table<2, double> &                weights,
                 ArrayView<Point<spacedim>>              new_points) const;
 
  virtual Point<spacedim>
  project_to_manifold(
    const ArrayView<const Point<spacedim>> &surrounding_points,
    const Point<spacedim> &                 candidate) const;
 
  virtual Point<spacedim>
  get_new_point_on_line(
    const typename Triangulation<dim, spacedim>::line_iterator &line) const;
 
  virtual Point<spacedim>
  get_new_point_on_quad(
    const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
 
  virtual Point<spacedim>
  get_new_point_on_hex(
    const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
 
 
  Point<spacedim>
  get_new_point_on_face(
    const typename Triangulation<dim, spacedim>::face_iterator &face) const;
 
 
  Point<spacedim>
  get_new_point_on_cell(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
 
  virtual Tensor<1, spacedim>
  get_tangent_vector(const Point<spacedim> &x1,
                     const Point<spacedim> &x2) const;
 
  virtual Tensor<1, spacedim>
  normal_vector(
    const typename Triangulation<dim, spacedim>::face_iterator &face,
    const Point<spacedim> &                                     p) const;
 
  virtual void
  get_normals_at_vertices(
    const typename Triangulation<dim, spacedim>::face_iterator &face,
    FaceVertexNormals &face_vertex_normals) const;
 
};
 
 
template <int dim, int spacedim = dim>
class FlatManifold : public Manifold<dim, spacedim>
{
public:
  FlatManifold(const Tensor<1, spacedim> &periodicity = Tensor<1, spacedim>(),
               const double               tolerance   = 1e-10);
 
  virtual std::unique_ptr<Manifold<dim, spacedim>>
  clone() const override;
 
  virtual Point<spacedim>
  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
                const ArrayView<const double> &         weights) const override;
 
 
  virtual void
  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                 const Table<2, double> &                weights,
                 ArrayView<Point<spacedim>> new_points) const override;
 
  virtual Point<spacedim>
  project_to_manifold(const ArrayView<const Point<spacedim>> &points,
                      const Point<spacedim> &candidate) const override;
 
  virtual Tensor<1, spacedim>
  get_tangent_vector(const Point<spacedim> &x1,
                     const Point<spacedim> &x2) const override;
 
  virtual Tensor<1, spacedim>
  normal_vector(
    const typename Triangulation<dim, spacedim>::face_iterator &face,
    const Point<spacedim> &p) const override;
 
  virtual void
  get_normals_at_vertices(
    const typename Triangulation<dim, spacedim>::face_iterator &face,
    typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
    const override;
 
  const Tensor<1, spacedim> &
  get_periodicity() const;
 
private:
  const Tensor<1, spacedim> periodicity;
 
  DeclException3(ExcPeriodicBox,
                 int,
                 Point<spacedim>,
                 double,
                 << "The component number " << arg1 << " of the point [ "
                 << arg2 << " ] is not in the interval [ 0, " << arg3
                 << "), bailing out.");
 
  const double tolerance;
};
 
 
template <int dim, int spacedim = dim, int chartdim = dim>
class ChartManifold : public Manifold<dim, spacedim>
{
public:
  // explicitly check for sensible template arguments
  static_assert(dim <= spacedim,
                "The dimension <dim> of a ChartManifold must be less than or "
                "equal to the space dimension <spacedim> in which it lives.");
 
  ChartManifold(const Tensor<1, chartdim> &periodicity = Tensor<1, chartdim>());
 
  virtual ~ChartManifold() override = default;
 
  virtual Point<spacedim>
  get_intermediate_point(const Point<spacedim> &p1,
                         const Point<spacedim> &p2,
                         const double           w) const override;
 
  virtual Point<spacedim>
  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
                const ArrayView<const double> &         weights) const override;
 
  virtual void
  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                 const Table<2, double> &                weights,
                 ArrayView<Point<spacedim>> new_points) const override;
  virtual Point<chartdim>
  pull_back(const Point<spacedim> &space_point) const = 0;
 
  virtual Point<spacedim>
  push_forward(const Point<chartdim> &chart_point) const = 0;
 
  virtual DerivativeForm<1, chartdim, spacedim>
  push_forward_gradient(const Point<chartdim> &chart_point) const;
 
  virtual Tensor<1, spacedim>
  get_tangent_vector(const Point<spacedim> &x1,
                     const Point<spacedim> &x2) const override;
 
  const Tensor<1, chartdim> &
  get_periodicity() const;
 
private:
  const FlatManifold<chartdim, chartdim> sub_manifold;
};
 
 
/* -------------- declaration of explicit specializations ------------- */
 
#ifndef DOXYGEN
 
template <>
Point<1>
Manifold<1, 1>::get_new_point_on_face(
  const Triangulation<1, 1>::face_iterator &) const;
 
template <>
Point<2>
Manifold<1, 2>::get_new_point_on_face(
  const Triangulation<1, 2>::face_iterator &) const;
 
 
template <>
Point<3>
Manifold<1, 3>::get_new_point_on_face(
  const Triangulation<1, 3>::face_iterator &) const;
 
 
template <>
Point<1>
Manifold<1, 1>::get_new_point_on_quad(
  const Triangulation<1, 1>::quad_iterator &) const;
 
template <>
Point<2>
Manifold<1, 2>::get_new_point_on_quad(
  const Triangulation<1, 2>::quad_iterator &) const;
 
 
template <>
Point<3>
Manifold<1, 3>::get_new_point_on_quad(
  const Triangulation<1, 3>::quad_iterator &) const;
 
 
template <>
Point<3>
Manifold<3, 3>::get_new_point_on_hex(
  const Triangulation<3, 3>::hex_iterator &) const;
 
/*---Templated functions---*/
 
namespace Manifolds
{
  template <typename MeshIteratorType>
  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
                       n_default_points_per_cell<MeshIteratorType>()>,
            std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
  get_default_points_and_weights(const MeshIteratorType &iterator,
                                 const bool              with_interpolation)
  {
    const int dim      = MeshIteratorType::AccessorType::structure_dimension;
    const int spacedim = MeshIteratorType::AccessorType::space_dimension;
    constexpr std::size_t points_per_cell =
      n_default_points_per_cell<MeshIteratorType>();
 
    std::pair<std::array<Point<spacedim>, points_per_cell>,
              std::array<double, points_per_cell>>
      points_weights;
 
 
    // note that the exact weights are chosen such as to minimize the
    // distortion of the four new quads from the optimal shape; their
    // derivation and values is copied over from the
    // interpolation function in the mapping
    switch (dim)
      {
        case 1:
          Assert(points_weights.first.size() == 2, ExcInternalError());
          Assert(points_weights.second.size() == 2, ExcInternalError());
          points_weights.first[0]  = iterator->vertex(0);
          points_weights.second[0] = .5;
          points_weights.first[1]  = iterator->vertex(1);
          points_weights.second[1] = .5;
          break;
        case 2:
          Assert(points_weights.first.size() == 8, ExcInternalError());
          Assert(points_weights.second.size() == 8, ExcInternalError());
 
          for (unsigned int i = 0; i < 4; ++i)
            {
              points_weights.first[i] = iterator->vertex(i);
              points_weights.first[4 + i] =
                (iterator->line(i)->has_children() ?
                   iterator->line(i)->child(0)->vertex(1) :
                   iterator->line(i)->get_manifold().get_new_point_on_line(
                     iterator->line(i)));
            }
 
          if (with_interpolation)
            {
              std::fill(points_weights.second.begin(),
                        points_weights.second.begin() + 4,
                        -0.25);
              std::fill(points_weights.second.begin() + 4,
                        points_weights.second.end(),
                        0.5);
            }
          else
            std::fill(points_weights.second.begin(),
                      points_weights.second.end(),
                      1.0 / 8.0);
          break;
        case 3:
          {
            TriaIterator<TriaAccessor<3, 3, 3>> hex =
              static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
            const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
                                    GeometryInfo<dim>::lines_per_cell +
                                    GeometryInfo<dim>::faces_per_cell;
            Assert(points_weights.first.size() == np, ExcInternalError());
            Assert(points_weights.second.size() == np, ExcInternalError());
            auto *sp3 = reinterpret_cast<
              std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
                *>(&points_weights.first);
 
            unsigned int j = 0;
 
            // note that the exact weights are chosen such as to minimize the
            // distortion of the eight new hexes from the optimal shape through
            // transfinite interpolation from the faces and vertices, see
            // TransfiniteInterpolationManifold for a deeper explanation of the
            // mechanisms
            if (with_interpolation)
              {
                for (unsigned int i = 0;
                     i < GeometryInfo<dim>::vertices_per_cell;
                     ++i, ++j)
                  {
                    (*sp3)[j]                = hex->vertex(i);
                    points_weights.second[j] = 1.0 / 8.0;
                  }
                for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
                     ++i, ++j)
                  {
                    (*sp3)[j] =
                      (hex->line(i)->has_children() ?
                         hex->line(i)->child(0)->vertex(1) :
                         hex->line(i)->get_manifold().get_new_point_on_line(
                           hex->line(i)));
                    points_weights.second[j] = -1.0 / 4.0;
                  }
                for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
                     ++i, ++j)
                  {
                    (*sp3)[j] =
                      (hex->quad(i)->has_children() ?
                         hex->quad(i)->isotropic_child(0)->vertex(3) :
                         hex->quad(i)->get_manifold().get_new_point_on_quad(
                           hex->quad(i)));
                    points_weights.second[j] = 1.0 / 2.0;
                  }
              }
            else
              // Overwrite the weights with 1/np if we don't want to use
              // interpolation.
              std::fill(points_weights.second.begin(),
                        points_weights.second.end(),
                        1.0 / np);
          }
          break;
        default:
          Assert(false, ExcInternalError());
          break;
      }
    return points_weights;
  }
} // namespace Manifolds
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif