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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/subscriptor.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/point.h>
#include <deal.II/base/derivative_form.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);
}

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.");


  typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];


  virtual ~Manifold () = 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;
  }
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif