manifold_lib.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_manifold_lib_h
#define dealii_manifold_lib_h


#include <deal.II/base/config.h>
#include <deal.II/grid/manifold.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>

#include <boost/container/small_vector.hpp>

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim = dim>
class PolarManifold : public ChartManifold<dim, spacedim, spacedim>
{
public:
  PolarManifold(const Point<spacedim> center = Point<spacedim>());

  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;

  virtual Point<spacedim>
  pull_back(const Point<spacedim> &space_point) const override;

  virtual Point<spacedim>
  push_forward(const Point<spacedim> &chart_point) const override;

  virtual
  DerivativeForm<1,spacedim,spacedim>
  push_forward_gradient(const Point<spacedim> &chart_point) const override;

  const Point<spacedim> center;
private:

  static Tensor<1,spacedim> get_periodicity();
};



template <int dim, int spacedim = dim>
class SphericalManifold : public Manifold<dim, spacedim>
{
public:
  SphericalManifold(const Point<spacedim> center = Point<spacedim>());

  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;

  virtual
  Point<spacedim>
  get_intermediate_point(const Point<spacedim> &p1,
                         const Point<spacedim> &p2,
                         const double w) 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;

  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>
  get_new_point (const ArrayView<const Point<spacedim>> &vertices,
                 const ArrayView<const double>          &weights) const override;

  const Point<spacedim> center;

private:
  std::pair<double, Tensor<1,spacedim> >
  guess_new_point(const ArrayView<const Tensor<1,spacedim>> &directions,
                  const ArrayView<const double> &distances,
                  const ArrayView<const double> &weights) const;

  Point<spacedim>
  get_new_point (const ArrayView<const Tensor<1,spacedim>> &directions,
                 const ArrayView<const double> &distances,
                 const ArrayView<const double> &weights,
                 const Point<spacedim> &candidate_point) const;

  virtual
  void
  get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
                  const ArrayView<const double>          &weights,
                  ArrayView<Point<spacedim>>              new_points) const;

  const PolarManifold<spacedim> polar_manifold;
};

template <int dim, int spacedim = dim>
class CylindricalManifold : public ChartManifold<dim,spacedim,3>
{
public:
  CylindricalManifold (const unsigned int axis = 0,
                       const double tolerance = 1e-10);

  CylindricalManifold (const Tensor<1, spacedim> &direction,
                       const Point<spacedim> &point_on_axis,
                       const double tolerance = 1e-10);

  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;

  virtual Point<3>
  pull_back(const Point<spacedim> &space_point) const override;

  virtual Point<spacedim>
  push_forward(const Point<3> &chart_point) const override;

  virtual DerivativeForm<1, 3, spacedim>
  push_forward_gradient(const Point<3> &chart_point) const override;

  virtual Point<spacedim>
  get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
                 const ArrayView<const double>          &weights) const override;

protected:
  const Tensor<1,spacedim> normal_direction;

  const Tensor<1,spacedim> direction;

  const Point<spacedim> point_on_axis;

private:
  double tolerance;

};



template <int dim, int spacedim=dim, int chartdim=dim>
class FunctionManifold : public ChartManifold<dim, spacedim, chartdim>
{
public:
  FunctionManifold(const Function<chartdim> &push_forward_function,
                   const Function<spacedim> &pull_back_function,
                   const Tensor<1,chartdim> &periodicity=Tensor<1,chartdim>(),
                   const double tolerance=1e-10);

  FunctionManifold(const std::string push_forward_expression,
                   const std::string pull_back_expression,
                   const Tensor<1,chartdim> &periodicity=Tensor<1,chartdim>(),
                   const typename FunctionParser<spacedim>::ConstMap = typename FunctionParser<spacedim>::ConstMap(),
                   const std::string chart_vars=FunctionParser<chartdim>::default_variable_names(),
                   const std::string space_vars=FunctionParser<spacedim>::default_variable_names(),
                   const double tolerance=1e-10,
                   const double h=1e-8);

  virtual ~FunctionManifold() override;

  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;

  virtual Point<spacedim>
  push_forward(const Point<chartdim> &chart_point) const override;

  virtual
  DerivativeForm<1,chartdim,spacedim>
  push_forward_gradient(const Point<chartdim> &chart_point) const override;

  virtual Point<chartdim>
  pull_back(const Point<spacedim> &space_point) const override;

private:
  const typename FunctionParser<spacedim>::ConstMap const_map;

  SmartPointer<const Function<chartdim>,
               FunctionManifold<dim,spacedim,chartdim> > push_forward_function;

  SmartPointer<const Function<spacedim>,
               FunctionManifold<dim,spacedim,chartdim> > pull_back_function;

  const double tolerance;

  const bool owns_pointers;

  const std::string push_forward_expression;

  const std::string pull_back_expression;

  const std::string chart_vars;

  const std::string space_vars;

  const double finite_difference_step;
};



template <int dim>
class TorusManifold : public ChartManifold<dim,3,3>
{
public:
  static const int chartdim = 3;
  static const int spacedim = 3;

  TorusManifold (const double R, const double r);

  virtual std::unique_ptr<Manifold<dim, 3> > clone() const override;

  virtual Point<3>
  pull_back(const Point<3> &p) const override;

  virtual Point<3>
  push_forward(const Point<3> &chart_point) const override;

  virtual
  DerivativeForm<1,3,3>
  push_forward_gradient(const Point<3> &chart_point) const override;

private:
  double r, R;
};



template <int dim, int spacedim=dim>
class TransfiniteInterpolationManifold : public Manifold<dim,spacedim>
{
public:
  TransfiniteInterpolationManifold();

  virtual ~TransfiniteInterpolationManifold() override;

  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;

  void initialize (const Triangulation<dim,spacedim> &triangulation);

  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;

private:
  std::array<unsigned int, 20>
  get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &surrounding_points) const;

  typename Triangulation<dim,spacedim>::cell_iterator
  compute_chart_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                       ArrayView<Point<dim>>                   chart_points) const;

  Point<dim>
  pull_back(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
            const Point<spacedim> &p,
            const Point<dim>      &initial_guess) const;

  Point<spacedim>
  push_forward(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
               const Point<dim> &chart_point) const;

  DerivativeForm<1,dim,spacedim>
  push_forward_gradient(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                        const Point<dim>      &chart_point,
                        const Point<spacedim> &pushed_forward_chart_point) const;

  const Triangulation<dim,spacedim> *triangulation;

  int level_coarse;

  std::vector<bool> coarse_cell_is_flat;

  FlatManifold<dim> chart_manifold;

  boost::signals2::connection clear_signal;
};

DEAL_II_NAMESPACE_CLOSE

#endif