composition_manifold.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_composition_manifold_h
#define dealii_composition_manifold_h


/*----------------------------   composition_manifold.h     ------------*/

#include <deal.II/base/config.h>

#include <deal.II/base/derivative_form.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/thread_management.h>

#include <deal.II/grid/manifold.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>

DEAL_II_NAMESPACE_OPEN

template <int dim,
          int spacedim         = dim,
          int chartdim         = dim,
          int intermediate_dim = dim,
          int dim1             = dim,
          int dim2             = dim>
class CompositionManifold : public ChartManifold<dim, spacedim, chartdim>
{
public:
  CompositionManifold(const ChartManifold<dim1, intermediate_dim, chartdim> &F,
                      const ChartManifold<dim2, spacedim, intermediate_dim> &G);

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

  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<dim1, intermediate_dim, chartdim>,
    CompositionManifold<dim, spacedim, chartdim, dim1, dim2, intermediate_dim>>
    F;


  SmartPointer<
    const ChartManifold<dim2, spacedim, intermediate_dim>,
    CompositionManifold<dim, spacedim, chartdim, dim1, dim2, intermediate_dim>>
    G;
};


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

template <int dim,
          int spacedim,
          int chartdim,
          int intermediate_dim,
          int dim1,
          int dim2>
CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>::
  CompositionManifold(const ChartManifold<dim1, intermediate_dim, chartdim> &F,
                      const ChartManifold<dim2, spacedim, intermediate_dim> &G)
  : ChartManifold<dim, spacedim, chartdim>(F.get_periodicity())
  , F(&F)
  , G(&G)
{
  // We don't know what to do with a periodicity in the second manifold, so
  // throw an assertion if the second manifold is periodic
  Assert(G.get_periodicity().norm() == 0.0,
         ExcMessage("The second manifold cannot be periodic."));
}


template <int dim,
          int spacedim,
          int chartdim,
          int intermediate_dim,
          int dim1,
          int dim2>
std::unique_ptr<Manifold<dim, spacedim>>
CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>::
  clone() const
{
  return std_cxx14::make_unique<
    CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>>(
    *F, *G);
}


template <int dim,
          int spacedim,
          int chartdim,
          int intermediate_dim,
          int dim1,
          int dim2>
Point<chartdim>
CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>::
  pull_back(const Point<spacedim> &space_point) const
{
  return F->pull_back(G->pull_back(space_point));
}



template <int dim,
          int spacedim,
          int chartdim,
          int intermediate_dim,
          int dim1,
          int dim2>
Point<spacedim>
CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>::
  push_forward(const Point<chartdim> &chart_point) const
{
  return G->push_forward(F->push_forward(chart_point));
}



template <int dim,
          int spacedim,
          int chartdim,
          int intermediate_dim,
          int dim1,
          int dim2>
DerivativeForm<1, chartdim, spacedim>
CompositionManifold<dim, spacedim, chartdim, intermediate_dim, dim1, dim2>::
  push_forward_gradient(const Point<chartdim> &chart_point) const
{
  DerivativeForm<1, chartdim, intermediate_dim> DF =
    F->push_forward_gradient(chart_point);

  DerivativeForm<1, intermediate_dim, spacedim> DG =
    G->push_forward_gradient(F->push_forward(chart_point));

  DerivativeForm<1, chartdim, spacedim> DF_DG;

  for (unsigned int d = 0; d < spacedim; ++d)
    for (unsigned int c = 0; c < chartdim; ++c)
      for (unsigned int s = 0; s < intermediate_dim; ++s)
        DF_DG[d][c] += DG[d][s] * DF[s][c];

  return DF_DG;
}


DEAL_II_NAMESPACE_CLOSE

#endif