From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 12 Oct 2000 13:37:56 +0000 (+0000)
Subject: New transform_real_to_unit_cell function in FiniteElement
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=695898633290e2652af3e19046db53be99c51dfd;p=dealii-svn.git

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 b3d0fc17c0..ebd8319e3a 100644
--- 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 9ee12608c7..74cbf06841 100644
--- 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 ee47fc6e63..62795946a2 100644
--- 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 883a6d7dfe..54979a122b 100644
--- 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 3633c5c784..e0e42046e6 100644
--- 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,