New transform_real_to_unit_cell function in FiniteElement

git-svn-id: https://svn.dealii.org/trunk@3433 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h
index b3d0fc17c055fdee53070fb3e2ee268eec0c6844..ebd8319e3ad2b649e0332500de07f7967bdc860e 100644 (file)
--- a/deal.II/deal.II/include/fe/fe.h
+++ b/deal.II/deal.II/include/fe/fe.h
@@ -1086,6 +1086,15 @@ class FiniteElement : public FiniteElementBase<dim>
     virtual Point<dim> transform_unit_to_real_cell (const DoFHandler<dim>::cell_iterator cell,
                                                    const Point<dim> &p) const = 0;
     
+                                    /**
+                                     * Transforms the point @p{p} on
+                                     * the real cell to the point
+                                     * @p{p_unit} on the unit cell
+                                     * @p{cell} and returns @p{p_unit}.
+                                     */
+    virtual Point<dim> transform_real_to_unit_cell (const DoFHandler<dim>::cell_iterator cell,
+                                                   const Point<dim> &p) const = 0;
+
                                     /**
                                      * Return the value of the
                                      * @p{i}th shape function of the

diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h
index 9ee12608c762adfecbdf256fdb334dc9151b78b5..74cbf0684186904adc1e06eeddafb77a718fb06b 100644 (file)
--- a/deal.II/deal.II/include/fe/fe_system.h
+++ b/deal.II/deal.II/include/fe/fe_system.h
@@ -267,6 +267,22 @@ class FESystem : public FiniteElement<dim>
     virtual Point<dim> transform_unit_to_real_cell (const DoFHandler<dim>::cell_iterator cell,
                                                    const Point<dim> &p) const;
 
+                                    /**
+                                     * Transforms the point @p{p} on
+                                     * the real cell to the point
+                                     * @p{p_unit} on the unit cell
+                                     * @p{cell} and returns
+                                     * @p{p_unit}. As the
+                                     * transformation mapping of each
+                                     * @p{FiniteElement} of this
+                                     * @p{FESystem} should be the
+                                     * same, this function just calls
+                                     * the @p{transform} function of
+                                     * @p{base_element(0)}.
+                                     */
+    virtual Point<dim> transform_real_to_unit_cell (const DoFHandler<dim>::cell_iterator cell,
+                                                   const Point<dim> &p) const;
+
                                     /**
                                      * Return the value of the @p{i}th shape
                                      * function of the transformation mapping

diff --git a/deal.II/deal.II/include/fe/q1_mapping.h b/deal.II/deal.II/include/fe/q1_mapping.h
index ee47fc6e63448a29f1eec64279ed6f30e6b9458b..62795946a2a3de9519081afc1506c553c459d6ce 100644 (file)
--- a/deal.II/deal.II/include/fe/q1_mapping.h
+++ b/deal.II/deal.II/include/fe/q1_mapping.h
@@ -56,6 +56,15 @@ class FEQ1Mapping : public FiniteElement<dim>
     virtual Point<dim> transform_unit_to_real_cell (const DoFHandler<dim>::cell_iterator cell,
                                                    const Point<dim> &p) const;
     
+                                    /**
+                                     * Transforms the point @p{p} on
+                                     * the real cell to the point
+                                     * @p{p_unit} on the unit cell
+                                     * @p{cell} and returns @p{p_unit}.
+                                     */
+    virtual Point<dim> transform_real_to_unit_cell (const DoFHandler<dim>::cell_iterator cell,
+                                                   const Point<dim> &p) const;
+    
                                     /**
                                      * Return the value of the @p{i}th shape
                                      * function at point @p{p} on the unit cell.

diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc
index 883a6d7dfe3a113b430ef89e6de88b0b8dd809d6..54979a122b6f899dcf2e05b58515744f10245014 100644 (file)
--- a/deal.II/deal.II/source/fe/fe_system.cc
+++ b/deal.II/deal.II/source/fe/fe_system.cc
@@ -933,6 +933,15 @@ Point<dim> FESystem<dim>::transform_unit_to_real_cell (
 };
 
 
+template <int dim>
+Point<dim> FESystem<dim>::transform_real_to_unit_cell (
+  const DoFHandler<dim>::cell_iterator cell,
+  const Point<dim> &p) const
+{
+  return base_elements[0].first->transform_real_to_unit_cell(cell, p);
+};
+
+
 template <int dim>
 double FESystem<dim>::shape_value_transform (const unsigned int i,
                                             const Point<dim>  &p) const

diff --git a/deal.II/deal.II/source/fe/q1_mapping.cc b/deal.II/deal.II/source/fe/q1_mapping.cc
index 3633c5c78443eaab4786d38bc007fc8b6c4e17a9..e0e42046e6a9f2d6a40fd14b512bd29233407d23 100644 (file)
--- a/deal.II/deal.II/source/fe/q1_mapping.cc
+++ b/deal.II/deal.II/source/fe/q1_mapping.cc
@@ -14,6 +14,8 @@
 
 #include <fe/q1_mapping.h>
 #include <base/quadrature.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 
@@ -543,6 +545,59 @@ Point<dim> FEQ1Mapping<dim>::transform_unit_to_real_cell (
 }
 
 
+template <int dim>
+Point<dim> FEQ1Mapping<dim>::transform_real_to_unit_cell (
+  const DoFHandler<dim>::cell_iterator cell,
+  const Point<dim> &p) const
+{
+                                  // Newton iteration to solve
+                                  // f(x)=p(x)-p=0
+                                  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+  
+                                  // start value in center of unit
+                                  // cell (p_unit stands for x)
+  Point<dim> p_unit;
+  for (unsigned int i=0; i<dim; ++i)
+    p_unit(i)=0.5;
+
+                                  // f(x)
+  Point<dim> p_real(transform_unit_to_real_cell(cell, p_unit));
+  Point<dim> f = p_real-p;
+
+  double eps=1e-15*cell->diameter();
+  
+  while (f.square()>eps*eps)
+    {
+                                      // f'(x)
+      Tensor<2,dim> df;
+      for (unsigned int k=0; k<GeometryInfo<dim>::vertices_per_cell; ++k)
+       {
+         Tensor<1,dim> grad_transform(shape_grad_transform (k, p_unit));
+         Point<dim> &vertex=cell->vertex(k);
+
+         
+         for (unsigned int i=0; i<dim; ++i)
+           for (unsigned int j=0; j<dim; ++j)
+             df[i][j]+=vertex[i]*grad_transform[j];
+       }
+      
+                                      // Solve  [f'(x)]d=f(x)
+      Point<dim> d;
+      Tensor<2,dim> df_1;
+
+      df_1 = invert(df);
+      contract (d, df_1, f);
+      
+      p_unit -= d;
+                                  // f(x)
+      p_real=transform_unit_to_real_cell(cell, p_unit);
+      f = p_real-p;
+    }
+  
+  return p_unit;
+}
+
+  
 template <int dim>
 void FEQ1Mapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                       const vector<Point<dim> >            &unit_points,