manifold.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2015 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/grid/tria.h>

DEAL_II_NAMESPACE_OPEN

template <int dim, int space_dim> class Triangulation;


namespace Manifolds
{
  template <typename OBJECT>
  Quadrature<OBJECT::AccessorType::space_dimension>
  get_default_quadrature(const OBJECT &obj, bool with_laplace = false);
}


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


  virtual ~Manifold ();

  virtual
  Point<spacedim>
  get_new_point(const Quadrature<spacedim> &quad) 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;
};


template <int dim, int spacedim=dim>
class FlatManifold: public Manifold<dim, spacedim>
{
public:
  FlatManifold (const Point<spacedim> periodicity=Point<spacedim>(),
                const double tolerance=1e-10);

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


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

  DeclException4(ExcPeriodicBox, int, Point<spacedim>, Point<spacedim>, double,
                 << "The component number " << arg1 << " of the point [ " << arg2
                 << " ]  is not in the interval [ " << -arg4
                 << ", " << arg3[arg4] << "), bailing out.");

  const double tolerance;
};


template <int dim, int spacedim=dim, int chartdim=dim>
class ChartManifold: public Manifold<dim,spacedim>
{
public:
  ChartManifold(const Point<chartdim> periodicity=Point<chartdim>());

  virtual ~ChartManifold ();


  virtual Point<spacedim>
  get_new_point(const Quadrature<spacedim> &quad) 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;

private:
  const FlatManifold<dim,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 OBJECT>
  Quadrature<OBJECT::AccessorType::space_dimension>
  get_default_quadrature(const OBJECT &obj, bool with_laplace)
  {
    const int spacedim = OBJECT::AccessorType::space_dimension;
    const int dim = OBJECT::AccessorType::structure_dimension;

    std::vector<Point<spacedim> > sp;
    std::vector<double> wp;


    // 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:
        sp.resize(2);
        wp.resize(2);
        sp[0] = obj->vertex(0);
        wp[0] = .5;
        sp[1] = obj->vertex(1);
        wp[1] = .5;
        break;
      case 2:
        sp.resize(8);
        wp.resize(8);

        for (unsigned int i=0; i<4; ++i)
          {
            sp[i] = obj->vertex(i);
            sp[4+i] = ( obj->line(i)->has_children() ?
                        obj->line(i)->child(0)->vertex(1) :
                        obj->line(i)->get_manifold().get_new_point_on_line(obj->line(i)) );
          }

        if (with_laplace)
          {
            std::fill(wp.begin(), wp.begin()+4, 1.0/16.0);
            std::fill(wp.begin()+4, wp.end(), 3.0/16.0);
          }
        else
          std::fill(wp.begin(), wp.end(), 1.0/8.0);
        break;
      case 3:
      {
        TriaIterator<TriaAccessor<3, 3, 3> > hex
          = static_cast<TriaIterator<TriaAccessor<3, 3, 3> > >(obj);
        const unsigned int np =
          GeometryInfo<dim>::vertices_per_cell+
          GeometryInfo<dim>::lines_per_cell+
          GeometryInfo<dim>::faces_per_cell;
        sp.resize(np);
        wp.resize(np);
        std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&sp);

        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);
            wp[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)));
            wp[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)));
            wp[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(wp.begin(), wp.end(), 1.0/np);
      }
      break;
      default:
        Assert(false, ExcInternalError());
        break;
      }
    return Quadrature<spacedim>(sp,wp);
  }
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif