manifold.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 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/base/thread_management.h>

#include <deal.II/grid/tria.h>

DEAL_II_NAMESPACE_OPEN

// forward declaration
template <int, typename>
class Table;

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>
  DEAL_II_DEPRECATED Quadrature<MeshIteratorType::AccessorType::space_dimension>
                     get_default_quadrature(const MeshIteratorType &iterator,
                                            const bool              with_interpolation = false);

  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>
  Quadrature<MeshIteratorType::AccessorType::space_dimension>
  get_default_quadrature(const MeshIteratorType &iterator,
                         const bool              with_interpolation)
  {
    const auto points_and_weights =
      get_default_points_and_weights(iterator, with_interpolation);
    static const int spacedim = MeshIteratorType::AccessorType::space_dimension;
    return Quadrature<spacedim>(
      std::vector<Point<spacedim>>(points_and_weights.first.begin(),
                                   points_and_weights.first.end()),
      std::vector<double>(points_and_weights.second.begin(),
                          points_and_weights.second.end()));
  }



  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