manifold.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 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/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>
  Quadrature<MeshIteratorType::AccessorType::space_dimension>
  get_default_quadrature(const MeshIteratorType &iterator,
                         const bool              with_laplace = false) DEAL_II_DEPRECATED;

  template <typename MeshIteratorType>
  std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
      std::vector<double> >
      get_default_points_and_weights(const MeshIteratorType &iterator,
                                     const bool              with_laplace = false);
}


template <int dim, int spacedim=dim>
class Manifold : public Subscriptor
{
public:

  // explicitly check for sensible template arguments
#ifdef DEAL_II_WITH_CXX11
  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.");
#endif


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


  virtual ~Manifold ();


  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 Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;

  virtual
  Point<spacedim>
  get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
                 const std::vector<double>           &weights) const;

  virtual
  void
  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
                  const Table<2,double>               &weights,
                  std::vector<Point<spacedim> >       &new_points) const;

  virtual
  Point<spacedim> project_to_manifold (const std::vector<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
  Point<spacedim>
  get_new_point(const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;

  virtual
  Point<spacedim>
  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                const std::vector<double>           &weights) const;

  virtual
  void
  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
                  const Table<2,double>               &weights,
                  std::vector<Point<spacedim> >       &new_points) const;

  virtual
  Point<spacedim>
  project_to_manifold (const std::vector<Point<spacedim> > &points,
                       const Point<spacedim> &candidate) const;

  virtual
  Tensor<1,spacedim>
  get_tangent_vector (const Point<spacedim> &x1,
                      const Point<spacedim> &x2) const;

  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
#ifdef DEAL_II_WITH_CXX11
  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.");
#endif

  ChartManifold(const Tensor<1,chartdim> &periodicity = Tensor<1,chartdim>());

  virtual ~ChartManifold ();


  virtual
  Point<spacedim>
  get_new_point(const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;

  virtual
  Point<spacedim>
  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                const std::vector<double>           &weights) const;

  virtual
  void
  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
                  const Table<2,double>               &weights,
                  std::vector<Point<spacedim> >       &new_points) const;

  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;

  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_laplace)
  {
    const std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
          std::vector<double> > points_and_weights = get_default_points_and_weights(iterator,
                                                     with_laplace);
    return Quadrature<MeshIteratorType::AccessorType::space_dimension>(points_and_weights.first,
           points_and_weights.second);
  }

  template <typename MeshIteratorType>
  std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
      std::vector<double> >
      get_default_points_and_weights(const MeshIteratorType &iterator,
                                     const bool              with_laplace)
  {
    const int spacedim = MeshIteratorType::AccessorType::space_dimension;
    const int dim = MeshIteratorType::AccessorType::structure_dimension;

    std::pair<std::vector<Point<spacedim> >,
        std::vector<double> > 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
    // @p{MappingQ::set_laplace_on_vector} function
    switch (dim)
      {
      case 1:
        points_weights.first.resize(2);
        points_weights.second.resize(2);
        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:
        points_weights.first.resize(8);
        points_weights.second.resize(8);

        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_laplace)
          {
            std::fill(points_weights.second.begin(), points_weights.second.begin()+4, 1.0/16.0);
            std::fill(points_weights.second.begin()+4, points_weights.second.end(), 3.0/16.0);
          }
        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;
        points_weights.first.resize(np);
        points_weights.second.resize(np);
        std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&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; their
        // derivation and values is copied over from the
        // @p{MappingQ::set_laplace_on_vector} function
        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
          {
            (*sp3)[j] = hex->vertex(i);
            points_weights.second[j] = 1.0/128.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] = 7.0/192.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/12.0;
          }
        // Overwrite the weights with 1/np if we don't want to use
        // laplace vectors.
        if (with_laplace == false)
          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