From ac3c38e839e76d07da0a2621643c037eee3d0bb5 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@gmail.com>
Date: Wed, 6 Aug 2014 16:01:41 +0200
Subject: [PATCH] Added FunctionManifoldChart

---
 include/deal.II/grid/manifold_lib.h | 99 +++++++++++++++++++++++++++++
 source/grid/manifold_lib.cc         | 86 +++++++++++++++++++++++++
 source/grid/manifold_lib.inst.in    |  3 +
 3 files changed, 188 insertions(+)

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index 8f54e8464f..462e137969 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -20,6 +20,8 @@
 
 #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>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -180,6 +182,103 @@ private:
   double tolerance;
 };
 
+
+/**
+ * Manifold description derived from ManifoldChart, based on explicit
+ * Function<spacedim> and Function<chartdim> objects describing the
+ * push_forward() and pull_back() functions. 
+ *
+ * You can use this Manifold object to describe any arbitray shape
+ * domain, as long as you can express it in terms of an invertible
+ * map, for which you provide both the forward expression, and the
+ * inverse expression.
+ *
+ * In debug mode, a check is performed to verify that the
+ * tranformations are actually one the inverse of the other.
+ * 
+ * @ingroup manifold
+ *
+ * @author Luca Heltai, 2014
+ */
+template <int dim, int spacedim=dim, int chartdim=dim>
+class FunctionManifoldChart : public ManifoldChart<dim, spacedim, chartdim> 
+{
+public:
+  /**
+   * Explicit functions constructor. Takes a push_forward function of
+   * spacedim components, and a pull_back function of chartdim
+   * components. See the documentation of the base class ManifoldChart
+   * for the meaning of the optional #periodicity argument.
+   */
+  FunctionManifoldChart(const Function<chartdim> &push_forward_function,
+			const Function<spacedim> &pull_back_function,
+			const Point<chartdim> periodicity=Point<chartdim>());
+
+  /**
+   * Expressions constructor. Takes the expressions of the
+   * push_forward function of spacedim components, and of the
+   * pull_back function of chartdim components. See the documentation
+   * of the base class ManifoldChart for the meaning of the optional
+   * #periodicity argument.
+   *
+   * The strings should be the readable by the default constructor of
+   * the FunctionParser classes. You can specify custom variable
+   * expressions with the last two optional arguments. If you don't,
+   * the default names are used, i.e., "x,y,z"
+   */
+  FunctionManifoldChart(const std::string push_forward_expression,
+			const std::string pull_back_expression,
+			const Point<chartdim> periodicity=Point<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());
+
+  /**
+   * If needed, we delete the pointers we own.
+   */
+  ~FunctionManifoldChart();
+  
+  /**
+   * Given a point in the chartdim coordinate system, uses the
+   * push_forward_function to compute the push_forward of points in
+   * #chardim space dimensions to #spacedim space dimensions.
+   */
+  virtual Point<spacedim>
+  push_forward(const Point<chartdim> &chart_point) const;
+
+  /**
+   * Given a point in the spacedim coordinate system, uses the
+   * pull_back_function to compute the pull_back of points in
+   * #spacedim space dimensions to #chartdim space dimensions.
+   */
+  virtual Point<chartdim>
+  pull_back(const Point<spacedim> &space_point) const;
+
+private:
+  /**
+   * Constants for the FunctionParser classes.
+   */
+  const typename FunctionParser<spacedim>::ConstMap const_map;
+
+  /**
+   * Pointer to the push_forward function.
+   */
+  SmartPointer<const Function<chartdim>,
+	       FunctionManifoldChart<dim,spacedim,chartdim> > push_forward_function;
+  
+  /**
+   * Pointer to the pull_back function.
+   */
+  SmartPointer<const Function<spacedim>, 
+	       FunctionManifoldChart<dim,spacedim,chartdim> > pull_back_function;
+
+  /**
+   * Check ownership of the smart pointers.
+   */
+  const bool owns_pointers;
+};
+  
+  
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index d2781df4d4..80a7e5003d 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -19,6 +19,7 @@
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/base/tensor.h>
+#include <deal.II/lac/vector.h>
 #include <cmath>
 
 DEAL_II_NAMESPACE_OPEN
@@ -202,6 +203,91 @@ get_new_point (const Quadrature<spacedim> &quad) const
             point_on_axis);
 }
 
+
+// ============================================================
+// FunctionManifoldChart
+// ============================================================
+template <int dim, int spacedim, int chartdim>
+FunctionManifoldChart<dim,spacedim,chartdim>::FunctionManifoldChart
+(const Function<chartdim> &push_forward_function,
+ const Function<spacedim> &pull_back_function,
+ const Point<chartdim> periodicity):
+  ManifoldChart<dim,spacedim,chartdim>(periodicity),
+  push_forward_function(&push_forward_function),
+  pull_back_function(&pull_back_function),
+  owns_pointers(false)
+{
+  AssertDimension(push_forward_function.n_components, spacedim);
+  AssertDimension(pull_back_function.n_components, chartdim);
+}
+
+template <int dim, int spacedim, int chartdim>
+FunctionManifoldChart<dim,spacedim,chartdim>::FunctionManifoldChart
+(const std::string push_forward_expression,
+ const std::string pull_back_expression, 
+ const Point<chartdim> periodicity, 
+ const typename FunctionParser<spacedim>::ConstMap const_map,
+ const std::string chart_vars, 
+ const std::string space_vars) :
+  ManifoldChart<dim,spacedim,chartdim>(periodicity),
+  const_map(const_map),
+  owns_pointers(true)
+{
+  FunctionParser<chartdim> * pf = new FunctionParser<chartdim>(spacedim);
+  FunctionParser<spacedim> * pb = new FunctionParser<spacedim>(chartdim);
+  pf->initialize(chart_vars, push_forward_expression, const_map);
+  pb->initialize(space_vars, pull_back_expression, const_map);
+  push_forward_function = pf;
+  pull_back_function = pb;
+}
+
+template <int dim, int spacedim, int chartdim>
+FunctionManifoldChart<dim,spacedim,chartdim>::~FunctionManifoldChart() {
+  if(owns_pointers == true) {
+    const Function<chartdim> * pf = push_forward_function;
+    push_forward_function = 0;
+    delete pf;
+    
+    const Function<spacedim> * pb = pull_back_function;
+    pull_back_function = 0;
+    delete pb;
+  }
+}
+								   
+template <int dim, int spacedim, int chartdim>
+Point<spacedim>
+FunctionManifoldChart<dim,spacedim,chartdim>::push_forward(const Point<chartdim> &chart_point) const
+{
+  Vector<double> pf(spacedim);
+  Point<spacedim> result;
+  push_forward_function->vector_value(chart_point, pf);
+  for (unsigned int i=0; i<spacedim; ++i)
+    result[i] = pf[i];
+
+#ifdef DEBUG
+  Vector<double> pb(chartdim);
+  pull_back_function->vector_value(result, pb);
+  for (unsigned int i=0; i<chartdim; ++i)
+    Assert(abs(pb[i]-chart_point[i]) < 1e-10*chart_point.norm(),
+           ExcMessage("The push forward is not the inverse of the pull back! Bailing out."));
+#endif
+
+  return result;
+}
+
+
+template <int dim, int spacedim, int chartdim>
+Point<chartdim>
+FunctionManifoldChart<dim,spacedim,chartdim>::pull_back(const Point<spacedim> &space_point) const
+{
+  Vector<double> pb(chartdim);
+  Point<chartdim> result;
+  pull_back_function->vector_value(space_point, pb);
+  for (unsigned int i=0; i<chartdim; ++i)
+    result[i] = pb[i];
+  return result;
+}
+
 // explicit instantiations
 #include "manifold_lib.inst"
 
diff --git a/source/grid/manifold_lib.inst.in b/source/grid/manifold_lib.inst.in
index 0d3616eb91..0945e6d66a 100644
--- a/source/grid/manifold_lib.inst.in
+++ b/source/grid/manifold_lib.inst.in
@@ -22,6 +22,9 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 #endif
 #if deal_II_dimension <= deal_II_space_dimension
     template class CylindricalManifold<deal_II_dimension, deal_II_space_dimension>;
+    template class FunctionManifoldChart<deal_II_dimension, deal_II_space_dimension, 1>;			  
+    template class FunctionManifoldChart<deal_II_dimension, deal_II_space_dimension, 2>;
+    template class FunctionManifoldChart<deal_II_dimension, deal_II_space_dimension, 3>;
 #endif
   }
 
-- 
2.39.5