tensor_product_manifold.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 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__tensor_product_manifold_h
#define dealii__tensor_product_manifold_h

#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/point.h>
#include <deal.II/grid/manifold.h>

DEAL_II_NAMESPACE_OPEN



template <int dim,
          int dim_A, int spacedim_A, int chartdim_A,
          int dim_B, int spacedim_B, int chartdim_B>
class TensorProductManifold :
  public ChartManifold<dim,spacedim_A+spacedim_B,chartdim_A+chartdim_B>
{
public:
  static const unsigned int chartdim = chartdim_A+chartdim_B;
  static const unsigned int spacedim = spacedim_A+spacedim_B;

  TensorProductManifold (
    const ChartManifold<dim_A, spacedim_A, chartdim_A> &manifold_A,
    const ChartManifold<dim_B, spacedim_B, chartdim_B> &manifold_B);

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

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

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

private:
  SmartPointer<const ChartManifold<dim_A, spacedim_A, chartdim_A>,
               TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B> > manifold_A;

  SmartPointer<const ChartManifold<dim_B, spacedim_B, chartdim_B>,
               TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B> > manifold_B;
};



/*------------------Template Implementations------------------------*/



namespace internal
{
  namespace TensorProductManifold
  {
    template <int dim1, int dim2>
    Tensor<1,dim1+dim2> concat(const Tensor<1,dim1> &p1, const Tensor<1,dim2> &p2)
    {
      Tensor<1,dim1+dim2> r;
      for (unsigned int d=0; d<dim1; ++d)
        r[d] = p1[d];
      for (unsigned int d=0; d<dim2; ++d)
        r[dim1+d] = p2[d];
      return r;
    }

    template <int dim1, int dim2>
    Point<dim1+dim2> concat(const Point<dim1> &p1, const Point<dim2> &p2)
    {
      Point<dim1+dim2> r;
      for (unsigned int d=0; d<dim1; ++d)
        r[d] = p1[d];
      for (unsigned int d=0; d<dim2; ++d)
        r[dim1+d] = p2[d];
      return r;
    }

    template <int dim1, int dim2>
    void split_point(const Point<dim1+dim2> &source, Point<dim1> &p1, Point<dim2> &p2)
    {
      for (unsigned int d=0; d<dim1; ++d)
        p1[d] = source[d];
      for (unsigned int d=0; d<dim2; ++d)
        p2[d] = source[dim1+d];
    }

  }
}

template <int dim,
          int dim_A, int spacedim_A, int chartdim_A,
          int dim_B, int spacedim_B, int chartdim_B>
TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>
::TensorProductManifold(
  const ChartManifold<dim_A, spacedim_A, chartdim_A> &manifold_A,
  const ChartManifold<dim_B, spacedim_B, chartdim_B> &manifold_B)
  : ChartManifold<dim,spacedim_A+spacedim_B,chartdim_A+chartdim_B> (
    internal::TensorProductManifold::concat(
      manifold_A.get_periodicity(),
      manifold_B.get_periodicity())),
  manifold_A (&manifold_A),
  manifold_B (&manifold_B)
{}

template <int dim,
          int dim_A, int spacedim_A, int chartdim_A,
          int dim_B, int spacedim_B, int chartdim_B>
Point<TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::chartdim>
TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>
::pull_back(const Point<TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::spacedim> &space_point) const
{
  Point<spacedim_A> space_point_A;
  Point<spacedim_B> space_point_B;
  internal::TensorProductManifold::split_point(space_point, space_point_A, space_point_B);

  Point<chartdim_A> result_A = manifold_A->pull_back(space_point_A);
  Point<chartdim_B> result_B = manifold_B->pull_back(space_point_B);

  return internal::TensorProductManifold::concat(result_A, result_B);
}

template <int dim,
          int dim_A, int spacedim_A, int chartdim_A,
          int dim_B, int spacedim_B, int chartdim_B>
Point<TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::spacedim>
TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>
::push_forward(const Point<TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::chartdim> &chart_point) const
{
  Point<chartdim_A> chart_point_A;
  Point<chartdim_B> chart_point_B;
  internal::TensorProductManifold::split_point(chart_point, chart_point_A, chart_point_B);

  Point<spacedim_A> result_A = manifold_A->push_forward(chart_point_A);
  Point<spacedim_B> result_B = manifold_B->push_forward(chart_point_B);

  return internal::TensorProductManifold::concat(result_A, result_B);
}

template <int dim,
          int dim_A, int spacedim_A, int chartdim_A,
          int dim_B, int spacedim_B, int chartdim_B>
DerivativeForm<1,
               TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::chartdim,
               TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::spacedim>

               TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>
               ::push_forward_gradient(const Point<TensorProductManifold<dim, dim_A, spacedim_A, chartdim_A, dim_B, spacedim_B, chartdim_B>::chartdim> &chart_point) const
{
  Point<chartdim_A> chart_point_A;
  Point<chartdim_B> chart_point_B;
  internal::TensorProductManifold::split_point(chart_point, chart_point_A, chart_point_B);

  DerivativeForm<1,chartdim_A,spacedim_A> result_A
    = manifold_A->push_forward_gradient(chart_point_A);
  DerivativeForm<1,chartdim_B,spacedim_B> result_B
    = manifold_B->push_forward_gradient(chart_point_B);


  DerivativeForm<1,chartdim,spacedim> result;
  for (unsigned int i = 0; i<chartdim_A; ++i)
    for (unsigned int j = 0; j<spacedim_A; ++j)
      result[j][i] = result_A[j][i];
  for (unsigned int i = 0; i<chartdim_B; ++i)
    for (unsigned int j = 0; j<spacedim_B; ++j)
      result[j+spacedim_A][i+chartdim_A] = result_B[j][i];

  return result;
}




DEAL_II_NAMESPACE_CLOSE

#endif