From 5d90c842817dc0f4775e3a4500a8e94f9822f634 Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 8 Jan 2009 00:00:16 +0000
Subject: [PATCH] transformations moved into Mapping

git-svn-id: https://svn.dealii.org/trunk@18127 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/fe/fe_poly.templates.h    |  67 +--
 .../deal.II/include/fe/fe_raviart_thomas.h    |   4 +-
 deal.II/deal.II/include/fe/mapping.h          |  50 ++-
 .../deal.II/include/fe/mapping_cartesian.h    |  27 +-
 deal.II/deal.II/include/fe/mapping_q.h        |  40 --
 deal.II/deal.II/include/fe/mapping_q1.h       |  24 --
 deal.II/deal.II/source/fe/fe_poly_tensor.cc   | 283 +++++--------
 .../deal.II/source/fe/fe_raviart_thomas.cc    |  20 +-
 .../source/fe/fe_raviart_thomas_nodal.cc      |  29 +-
 deal.II/deal.II/source/fe/mapping.cc          |  65 ++-
 .../deal.II/source/fe/mapping_cartesian.cc    | 252 +++++-------
 deal.II/deal.II/source/fe/mapping_q.cc        | 212 ++--------
 deal.II/deal.II/source/fe/mapping_q1.cc       | 381 +++++++-----------
 13 files changed, 546 insertions(+), 908 deletions(-)

diff --git a/deal.II/deal.II/include/fe/fe_poly.templates.h b/deal.II/deal.II/include/fe/fe_poly.templates.h
index b579d00b4e..337cf7fbf1 100644
--- a/deal.II/deal.II/include/fe/fe_poly.templates.h
+++ b/deal.II/deal.II/include/fe/fe_poly.templates.h
@@ -260,14 +260,8 @@ FE_Poly<TensorProductPolynomials<1>,1,2>::fill_fe_values (const Mapping<1,2> &ma
 	  data.shape_values(k,i) = fe_data.shape_values[k][i];
       
       if (flags & update_gradients)
-	{
-	  Assert (data.shape_gradients[k].size() <=
- 		  fe_data.shape_gradients[k].size(), 
- 		  ExcInternalError()); 
- 	  mapping.transform(fe_data.shape_gradients[k],
-			    data.shape_gradients[k], 
-			    mapping_data, mapping_covariant); 
-	}
+	mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], 
+			  mapping_data, mapping_covariant);
     }
   
   //  const typename QProjector<1>::DataSetDescriptor dsd;
@@ -303,14 +297,8 @@ FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_values (const Mapping<2,3> &ma
 	  data.shape_values(k,i) = fe_data.shape_values[k][i];     
 
        if (flags & update_gradients)
-       {
-	 Assert (data.shape_gradients[k].size() <=
-		 fe_data.shape_gradients[k].size(), 
-		 ExcInternalError());
-	 mapping.transform(fe_data.shape_gradients[k],
-			   data.shape_gradients[k],
+	 mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
 			   mapping_data, mapping_covariant);
- 	}
     }
   
 /*   const QProjector<2>::DataSetDescriptor dsd; */
@@ -327,12 +315,13 @@ FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_values (const Mapping<2,3> &ma
 template <>
 inline
 void
-FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (const Mapping<1,2> &/* mapping */,
-				   const Triangulation<1,2>::cell_iterator   &/* cell */,
-				   const Quadrature<1>                       &quadrature,
-				   Mapping<1,2>::InternalDataBase            &/* mapping_data */,
-				   Mapping<1,2>::InternalDataBase            &fedata,
-				   FEValuesData<1,2>                         &data) const
+FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (
+  const Mapping<1,2>& mapping,
+  const Triangulation<1,2>::cell_iterator& /*cell*/,
+  const Quadrature<1>& quadrature,
+  Mapping<1,2>::InternalDataBase& mapping_data,
+  Mapping<1,2>::InternalDataBase& fedata,
+  FEValuesData<1,2>& data) const
 {
 				   // convert data object to internal
 				   // data for this class. fails with
@@ -353,12 +342,8 @@ FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (const Mapping<1,2> &/* mapping
 // TODO: I would think this should work. Guido
       
       if (flags & update_gradients)
-	{
-	  AssertThrow(false, ExcNotImplemented());
-//  	  mapping.transform(fe_data.shape_gradients[k], 0,
-// 			    data.shape_gradients[k],
-// 			    mapping_data, mapping_covariant);
-	}
+	mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
+			  mapping_data, mapping_covariant);
     }
   
 				   //  const typename QProjector<1>::DataSetDescriptor dsd;
@@ -374,10 +359,10 @@ FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (const Mapping<1,2> &/* mapping
 template <>
 inline
 void
-FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values (const Mapping<2,3> &/* mapping */,
+FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values (const Mapping<2,3>& mapping,
 				   const Triangulation<2,3>::cell_iterator &/* cell */,
-				   const Quadrature<2> &quadrature,
-				   Mapping<2,3>::InternalDataBase &/* mapping_data */,
+				   const Quadrature<2>& quadrature,
+				   Mapping<2,3>::InternalDataBase& mapping_data,
 				   Mapping<2,3>::InternalDataBase &fedata,
 				   FEValuesData<2,3> &data) const
 {
@@ -394,15 +379,8 @@ FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values (const Mapping<2,3> &/* mapping
       
 
        if (flags & update_gradients)
- 	{
-	  AssertThrow(false, ExcNotImplemented());
-/*  	  Assert (data.shape_gradients[k].size()  */
-/*  		  fe_data.shape_gradients[k].size(),  */
-/*  		  ExcInternalError());	    */
-/*  	  mapping.transform(fe_data.shape_gradients[k], 0,  */
-/*  				      data.shape_gradients[k],  */
-/*  				      mapping_data, mapping_covariant);  */
- 	}
+	 mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
+  				      mapping_data, mapping_covariant);
     }
   
 /*   const QProjector<2>::DataSetDescriptor dsd; */
@@ -670,14 +648,9 @@ FE_Poly<POLY,dim,spacedim>::fill_fe_subface_values (const Mapping<dim,spacedim>
 	  data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
       
       if (flags & update_gradients)
-	{
-	  Assert (data.shape_gradients[k].size() + offset <=
-		  fe_data.shape_gradients[k].size(),
-		  ExcInternalError());
-	  mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
-			    data.shape_gradients[k],
-			    mapping_data, mapping_covariant);
-	}
+	mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
+			  data.shape_gradients[k],
+			  mapping_data, mapping_covariant);
     }
   
   if (flags & update_hessians)
diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h
index 6895cfc790..fdb283259d 100644
--- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h
+++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -442,6 +442,8 @@ class FE_RaviartThomasNodal
     virtual FiniteElementDomination::Domination
     compare_for_face_domination (const FiniteElement<dim> &fe_other) const;
 
+    virtual UpdateFlags update_once (const UpdateFlags flags) const;
+    virtual UpdateFlags update_each (const UpdateFlags flags) const;
   private:
 				     /**
 				      * Only for internal use. Its
diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h
index 51023ab054..095d4e6204 100644
--- a/deal.II/deal.II/include/fe/mapping.h
+++ b/deal.II/deal.II/include/fe/mapping.h
@@ -64,10 +64,43 @@ enum MappingType
       mapping_covariant = 0x0001,
 /// Contravariant mapping
       mapping_contravariant = 0x0002,
+/// Mapping of the gradient of a covariant vector field
+      mapping_covariant_gradient = 0x0003,
+/// Mapping of the gradient of a contravariant vector field
+      mapping_contravariant_gradient = 0x0004,
 /// The Piola transform usually used for Hdiv elements
+				       /**
+					* Piola transform is the
+					* standard transformation of
+					* vector valued elements in
+					* H<sup>div</sup>. It amounts
+					* to a contravariant
+					* transformation scaled by the
+					* inverse of the volume
+					* element.
+					*
+					* If applied to a rank 2
+					* tensor, the mapping class
+					* will apply the correct
+					* transformation for the
+					* gradient of such a vector
+					* field.
+					*/
       mapping_piola = 0x0100,
-/// The Piola transform for the derivative of an Hdiv element
-      mapping_piola_gradient = 0x0101,
+/// The mapping used for Nedelec elements
+				       /**
+					* curl-conforming elements are
+					* mapped as covariant
+					* vectors. Nevertheless, we
+					* introduce a separate mapping
+					* type, such that we can use
+					* the same flag for the vector
+					* and its gradient (the
+					* transformation of the
+					* gradient differs from the
+					* one used by #mapping_covariant).
+					*/
+      mapping_nedelec = 0x0200,
 /// The mapping used for Raviart-Thomas elements
       mapping_raviart_thomas = mapping_piola,
 /// The mapping used for BDM elements
@@ -332,42 +365,39 @@ class Mapping : public Subscriptor
 				     /**
 				      * @deprecated Use transform() instead.
 				      */
-    virtual
     void
     transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
                          const unsigned int                                    offset,
                          VectorSlice<std::vector<Tensor<1,spacedim> > >        output,
-			 const InternalDataBase &internal) const = 0;
+			 const InternalDataBase &internal) const;
 
                                      /**
 				      * @deprecated Use transform() instead.
                                       */
-    virtual
     void
     transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
                          const unsigned int                 offset,
 			 VectorSlice<std::vector<Tensor<2,spacedim> > >      output,
-			 const InternalDataBase &internal) const = 0;
+			 const InternalDataBase &internal) const;
 
 				     /**
 				      * @deprecated Use transform() instead.
 				      */
-    virtual
     void
     transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
                              const unsigned int                 offset,
 			     VectorSlice<std::vector<Tensor<1,spacedim> > >      output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const = 0;
+			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
 
                                      /**
 				      * @deprecated Use transform() instead.
                                       */
     
-    virtual void
+    void
     transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > intput,
                              const unsigned int                 offset,
 			     const VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const = 0;
+			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
 
 				     /**
 				      * The transformed (generalized)
diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h
index 3056a686c7..e15022238b 100644
--- a/deal.II/deal.II/include/fe/mapping_cartesian.h
+++ b/deal.II/deal.II/include/fe/mapping_cartesian.h
@@ -102,33 +102,8 @@ class MappingCartesian : public Mapping<dim,spacedim>
     transform (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
                VectorSlice<std::vector<Tensor<2,spacedim> > > output,
                const typename Mapping<dim,spacedim>::InternalDataBase &internal,
-	       const MappingType type) const;
-    
-
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			 const typename Mapping<dim, spacedim>::InternalDataBase &internal) const;
-    
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			 const typename Mapping<dim, spacedim>::InternalDataBase &internal) const;
-    
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			     const typename Mapping<dim, spacedim>::InternalDataBase &internal) const;
+	       const MappingType type) const;    
     
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			     const typename Mapping<dim, spacedim>::InternalDataBase &internal) const;
-
     virtual Point<spacedim>
     transform_unit_to_real_cell (
       const typename Triangulation<dim,spacedim>::cell_iterator &cell,
diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h
index 288f12ce32..13e7a4546a 100644
--- a/deal.II/deal.II/include/fe/mapping_q.h
+++ b/deal.II/deal.II/include/fe/mapping_q.h
@@ -128,46 +128,6 @@ class MappingQ : public MappingQ1<dim,spacedim>
                const typename Mapping<dim,spacedim>::InternalDataBase &internal,
 	       const MappingType type) const;
     
-				     /**
-				      * Implementation of the
-				      * interface in Mapping.
-				      */
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<1,spacedim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			 const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
-				     /**
-				      * Implementation of the
-				      * interface in Mapping.
-				      */
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<2,spacedim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			 const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
-				     /**
-				      * Implementation of the
-				      * interface in Mapping.
-				      */
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;    
-
-				     /**
-				      * Implementation of the
-				      * interface in Mapping.
-				      */
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;    
-    
 				     /**
 				      * Return the degree of the
 				      * mapping, i.e. the value which
diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h
index fd9a81d643..0e059bb1dc 100644
--- a/deal.II/deal.II/include/fe/mapping_q1.h
+++ b/deal.II/deal.II/include/fe/mapping_q1.h
@@ -89,30 +89,6 @@ class MappingQ1 : public Mapping<dim,spacedim>
                const typename Mapping<dim,spacedim>::InternalDataBase &internal,
 	       const MappingType type) const;
     
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			 const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
-    virtual void
-    transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-                         const unsigned int                 offset,
-			 VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			 const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
-    virtual void
-    transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-                             const unsigned int                 offset,
-			     VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-			     const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-    
                                      /**
                                       * Return a pointer to a copy of the
                                       * present object. The caller of this
diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc
index 3ca14a3a4e..1ce943a7c7 100644
--- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc
+++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc
@@ -407,47 +407,27 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 	      case mapping_covariant:
 	      case mapping_contravariant:
 	      {
-						 // Use auxiliary vector for
-						 // transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == mapping_covariant)
-		  mapping.transform(fe_data.shape_values[i],
-				    shape_values, mapping_data, mapping_covariant);
-		else
-		  mapping.transform(fe_data.shape_values[i],
-				    shape_values, mapping_data, mapping_contravariant);
-		    
-						 // then copy over to target:
+		mapping.transform(fe_data.shape_values[i],
+				  shape_values, mapping_data, mapping_type);
+		
 		for (unsigned int k=0; k<n_q_points; ++k)
-		  {
-//TODO: [GK] Remove this after change to Piola complete
-		    
-						     // Recompute determinant
-
-		    const double
-		      J = (mapping_type == mapping_contravariant ?
-			   data.JxW_values[k] / quadrature.weight(k)
-			   :
-			   1.0);
-
-		    for (unsigned int d=0; d<dim; ++d)
-		      data.shape_values(first+d,k)
-			= sign_change[i] * (shape_values[k][d] / J);
-		  }
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_values(first+d,k) = shape_values[k][d];
+		
 		break;
 	      }
 	      case mapping_piola:
-		    if (true)
-		      {
-			std::vector<Tensor<1,dim> > shape_values (n_q_points);
-			mapping.transform(fe_data.shape_values[i], shape_values,
-					  mapping_data, mapping_piola);
-			for (unsigned int k=0; k<n_q_points; ++k)
-			  for (unsigned int d=0; d<dim; ++d)
-			    data.shape_values(first+d,k)
-			      = sign_change[i] * shape_values[k][d];
-			break;
-		      }
+	      {
+		std::vector<Tensor<1,dim> > shape_values (n_q_points);
+		mapping.transform(fe_data.shape_values[i], shape_values,
+				  mapping_data, mapping_piola);
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_values(first+d,k)
+		      = sign_change[i] * shape_values[k][d];
+		break;
+	      }
 	      
 	      default:
 		Assert(false, ExcNotImplemented());
@@ -470,47 +450,40 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
 		break;
 	      }
-	      
 	      case mapping_covariant:
 	      {
 		mapping.transform(fe_data.shape_grads[i], shape_grads1,
-				  mapping_data, mapping_covariant);
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
-						 // do second transformation
-		mapping.transform(shape_grads2, shape_grads1,
-				  mapping_data, mapping_covariant);
-						 // transpose back
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
+				  mapping_data, mapping_covariant_gradient);
+		
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
+		
+		break;
+	      }
+	      case mapping_contravariant:
+	      {
+		mapping.transform(fe_data.shape_grads[i], shape_grads1,
+				  mapping_data, mapping_contravariant_gradient);
 		
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
 		    data.shape_gradients[first+d][k] = shape_grads2[k][d];
-
+		
 		break;
 	      }
 	      case mapping_piola:
-	      case mapping_contravariant:
 	      {
 		mapping.transform(fe_data.shape_grads[i], shape_grads1,
-				  mapping_data, mapping_covariant);
+				  mapping_data, mapping_piola);
 		
-		mapping.transform(shape_grads1, shape_grads2,
-				  mapping_data, mapping_contravariant);
-
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
-		    {
-						       // Recompute determinant
-		      const double J = data.JxW_values[k] / quadrature.weight(k);
-		      data.shape_gradients[first+d][k] = sign_change[i] * 
-			shape_grads2[k][d] / J;
-		    }
+		    data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d];
 		
 		break;
 	      }
-	      
+		    
 	      default:
 		Assert(false, ExcNotImplemented());
 	    }
@@ -593,52 +566,26 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 						 // Use auxiliary vector
 						 // for transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == mapping_covariant)
-		  mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-				    shape_values, mapping_data, mapping_covariant);
-		else
-		  mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-				    shape_values, mapping_data, mapping_contravariant);
+		mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+				  shape_values, mapping_data, mapping_type);
 		
-						 // then copy over to target:
 		for (unsigned int k=0; k<n_q_points; ++k)
-		  {
-						     // recompute
-						     // determinant. note that
-						     // we have to take the
-						     // *cell* transformation
-						     // (cell_JxW_values)
-						     // because the usual
-						     // JxW_values field
-						     // contains the
-						     // transformation of the
-						     // face in
-						     // fill_fe_face_values
-		    const double
-		      J = (mapping_type == mapping_contravariant ?
-			   data.cell_JxW_values[k] / quadrature.weight(k)
-			   :
-			   1.0);
-		    
-		    for (unsigned int d=0; d<dim; ++d)
-		      data.shape_values(first+d,k)
-			= sign_change[i] * (shape_values[k][d] / J);
-		  }
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_values(first+d,k) = shape_values[k][d];
 		
 		break;
 	      }
 	      case mapping_piola:
-		    if (true)
-		      {
-			std::vector<Tensor<1,dim> > shape_values (n_q_points);
-			mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-					  shape_values, mapping_data, mapping_piola);
-			for (unsigned int k=0; k<n_q_points; ++k)
-			  for (unsigned int d=0; d<dim; ++d)
-			    data.shape_values(first+d,k)
-			      = sign_change[i] * shape_values[k][d];
-			break;
-		      }
+	      {
+		std::vector<Tensor<1,dim> > shape_values (n_q_points);
+		mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+				  shape_values, mapping_data, mapping_piola);
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_values(first+d,k)
+		      = sign_change[i] * shape_values[k][d];
+		break;
+	      }
 	      
 	      default:
 		Assert(false, ExcNotImplemented());
@@ -665,40 +612,33 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 	      case mapping_covariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
-				  shape_grads1, mapping_data, mapping_covariant);
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
-						 // do second transformation
-		mapping.transform(shape_grads2, shape_grads1,
-				  mapping_data, mapping_covariant);
-						 // transpose back
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
+				  shape_grads1, mapping_data, mapping_covariant_gradient);
 		
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
 		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
 		break;
 	      }
-	      case mapping_piola:
 	      case mapping_contravariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
-				  shape_grads1, mapping_data, mapping_covariant);
-
-		mapping.transform(shape_grads1, shape_grads2,
-				  mapping_data, mapping_contravariant);
+				  shape_grads1, mapping_data, mapping_contravariant_gradient);
 
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
-		    {
-						       // recompute
-						       // determinant in the
-						       // same way as above
-		      const double J = data.cell_JxW_values[k] / quadrature.weight(k);
-		      data.shape_gradients[first+d][k]
-			= sign_change[i] * shape_grads2[k][d] / J;
-		    }
+		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
+		
+		break;
+	      }
+	      case mapping_piola:
+	      {
+		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
+				  shape_grads1, mapping_data, mapping_piola);
+		
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d];
+		
 		break;
 	      }
 	      
@@ -783,52 +723,26 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 						 // Use auxiliary vector for
 						 // transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == mapping_covariant)
-		  mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-				    shape_values, mapping_data, mapping_covariant);
-		else
-		  mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-				    shape_values, mapping_data, mapping_contravariant);
-		    
-						 // then copy over to target:
+		mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+				  shape_values, mapping_data, mapping_type);
+		
 		for (unsigned int k=0; k<n_q_points; ++k)
-		  {
-						     // recompute
-						     // determinant. note that
-						     // we have to take the
-						     // *cell* transformation
-						     // (cell_JxW_values)
-						     // because the usual
-						     // JxW_values field
-						     // contains the
-						     // transformation of the
-						     // face in
-						     // fill_fe_face_values
-		    const double
-		      J = (mapping_type == mapping_contravariant ?
-			   data.cell_JxW_values[k] / quadrature.weight(k)
-			   :
-			   1.0);
-		    
 		    for (unsigned int d=0; d<dim; ++d)
-		      data.shape_values(first+d,k)
-			= sign_change[i] * (shape_values[k][d] / J);
-		  }
-
+		      data.shape_values(first+d,k) = shape_values[k][d];
+		
 		break;
 	      }
 	      case mapping_piola:
-		    if (true)
-		      {
-			std::vector<Tensor<1,dim> > shape_values (n_q_points);
-			mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
-					  shape_values, mapping_data, mapping_piola);
-			for (unsigned int k=0; k<n_q_points; ++k)
-			  for (unsigned int d=0; d<dim; ++d)
-			    data.shape_values(first+d,k)
-			      = sign_change[i] * shape_values[k][d];
-			break;
-		      }
+	      {
+		std::vector<Tensor<1,dim> > shape_values (n_q_points);
+		mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+				  shape_values, mapping_data, mapping_piola);
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_values(first+d,k)
+		      = sign_change[i] * shape_values[k][d];
+		break;
+	      }
 	      default:
 		Assert(false, ExcNotImplemented());
 	    }
@@ -854,46 +768,41 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 	      case mapping_covariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
-				  shape_grads1, mapping_data, mapping_covariant);
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
-						 // do second transformation
-		mapping.transform(shape_grads2, shape_grads1,
-				  mapping_data, mapping_covariant);
-						 // transpose back
-		for (unsigned int q=0; q<n_q_points; ++q)
-		  shape_grads2[q] = transpose(shape_grads1[q]);
-		
+				  shape_grads1, mapping_data, mapping_covariant_gradient);
+
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
 		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
 		break;
 	      }
 
-	      case mapping_piola:
 	      case mapping_contravariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
-				  shape_grads1, mapping_data, mapping_covariant);
-
-		mapping.transform(shape_grads1, shape_grads2,
-				  mapping_data, mapping_contravariant);
-
+				  shape_grads1, mapping_data, mapping_contravariant_gradient);
+		
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
-		    {
-						       // recompute
-						       // determinant in the
-						       // same way as above
-		      const double J = data.cell_JxW_values[k] / quadrature.weight(k);
-		      data.shape_gradients[first+d][k]
-			= sign_change[i] * shape_grads2[k][d] / J;
-		    }
+		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
+		
 		break;
 	      }
 
+	      case mapping_piola:
+	      {
+		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
+				  shape_grads1,
+				  mapping_data, mapping_piola);
+		
+		for (unsigned int k=0; k<n_q_points; ++k)
+		  for (unsigned int d=0; d<dim; ++d)
+		    data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d];
+		
+		break;
+	      }
+	      
 	      default:
-		Assert(false, ExcNotImplemented());
+		    Assert(false, ExcNotImplemented());
 	    }
 	}
     }
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
index 900ebb1c15..cce73d0fcc 100644
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
@@ -50,7 +50,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
   
-  this->mapping_type = mapping_piola;
+  this->mapping_type = mapping_raviart_thomas;
 				   // First, initialize the
 				   // generalized support points and
 				   // quadrature weights, since they
@@ -463,20 +463,14 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
   UpdateFlags out = update_default;
 
   if (flags & update_values)
-    out |= update_values             | update_covariant_transformation
-                                     | update_piola 
-                                     | update_cell_JxW_values
-                                     | update_JxW_values;
+    out |= update_values | update_piola;
+  
   if (flags & update_gradients)
-    out |= update_gradients          | update_covariant_transformation 
-                                     | update_piola
-                                     | update_cell_JxW_values
-                                     | update_JxW_values;
+    out |= update_gradients | update_piola | update_covariant_transformation;
+  
   if (flags & update_hessians)
-    out |= update_hessians
-	   | update_contravariant_transformation
-	   | update_covariant_transformation;
-
+    out |= update_hessians | update_piola | update_covariant_transformation;
+  
   return out;
 }
 
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
index 8af15830ad..9c0087f7cd 100644
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
@@ -42,7 +42,7 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
   
-  this->mapping_type = mapping_none;
+  this->mapping_type = mapping_raviart_thomas;
 				   // These must be done first, since
 				   // they change the evaluation of
 				   // basis functions
@@ -770,6 +770,33 @@ FE_RaviartThomasNodal<dim>::initialize_support_points (const unsigned int deg)
 }
 
 
+template <int dim>
+UpdateFlags
+FE_RaviartThomasNodal<dim>::update_once (const UpdateFlags) const
+{
+  return update_default;
+}
+
+
+template <int dim>
+UpdateFlags
+FE_RaviartThomasNodal<dim>::update_each (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+  
+  if (flags & update_values)
+    out |= update_values | update_piola;
+  
+  if (flags & update_gradients)
+    out |= update_gradients | update_piola | update_covariant_transformation;
+  
+  if (flags & update_hessians)
+    out |= update_hessians | update_piola | update_covariant_transformation;
+  
+  return out;
+}
+
+
 template class FE_RaviartThomasNodal<deal_II_dimension>;
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/fe/mapping.cc b/deal.II/deal.II/source/fe/mapping.cc
index 5a2ae42186..e0d8290ede 100644
--- a/deal.II/deal.II/source/fe/mapping.cc
+++ b/deal.II/deal.II/source/fe/mapping.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2005, 2006, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -22,6 +22,69 @@ Mapping<dim, spacedim>::~Mapping ()
 {}
 
 
+// This function is deprecated and has been replaced by transform above
+template<int dim, int spacedim>
+void
+Mapping<dim,spacedim>::transform_covariant (
+  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
+  const unsigned int                 offset,
+  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
+{
+  Assert (offset == 0, ExcInternalError());
+
+  transform(input, output, mapping_data, mapping_covariant);
+}
+
+
+
+// This function is deprecated and has been replaced by transform above
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::transform_covariant (
+    const VectorSlice<const std::vector<Tensor<2,dim> > > input,
+    const unsigned int                 offset,
+    VectorSlice<std::vector<Tensor<2,spacedim> > > output,
+    const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
+{
+  Assert (offset == 0, ExcInternalError());
+  
+  transform(input, output, mapping_data, mapping_covariant);  
+}
+
+
+
+// This function is deprecated and has been replaced by transform above
+template<int dim, int spacedim>
+void
+Mapping<dim,spacedim>::transform_contravariant (
+  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
+  const unsigned int                 offset,
+  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
+{
+  Assert (offset == 0, ExcInternalError());
+
+  transform(input, output, mapping_data, mapping_contravariant);
+}
+
+
+
+// This function is deprecated and has been replaced by transform above
+template<int dim, int spacedim>
+void
+Mapping<dim,spacedim>::transform_contravariant (
+  const VectorSlice<const std::vector<Tensor<2,dim> > > input,
+  const unsigned int                 offset,
+  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
+{
+  Assert (offset == 0, ExcInternalError());
+
+  transform(input, output, mapping_data, mapping_contravariant);
+}
+
+
 /*------------------------------ InternalDataBase ------------------------------*/
 
 
diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc
index c9f4dbb9c3..20fad4a21d 100644
--- a/deal.II/deal.II/source/fe/mapping_cartesian.cc
+++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc
@@ -534,40 +534,38 @@ MappingCartesian<dim,spacedim>::transform (
   switch (mapping_type)
     {
       case mapping_covariant:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_covariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  for (unsigned int d=0;d<dim;++d)
-		    output[i][d] = input[i][d]/data.length[d];
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d=0;d<dim;++d)
+	    output[i][d] = input[i][d]/data.length[d];
+	return;
+      }
+      
       case mapping_contravariant:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_contravariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  for (unsigned int d=0;d<dim;++d)
-		    output[i][d] = input[i][d]*data.length[d];
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d=0;d<dim;++d)
+	    output[i][d] = input[i][d]*data.length[d];
+	return;
+      }
       case mapping_piola:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_contravariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		Assert (data.update_flags & update_volume_elements,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  for (unsigned int d=0;d<dim;++d)
-		    output[i][d] = input[i][d] * data.length[d] / data.volume_element;
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_volume_elements,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d=0;d<dim;++d)
+	    output[i][d] = input[i][d] * data.length[d] / data.volume_element;
+	return;
+      }
       default:
 	    Assert(false, ExcNotImplemented());
     }
@@ -578,136 +576,82 @@ void
 MappingCartesian<dim,spacedim>::transform (
   const VectorSlice<const std::vector<Tensor<2,dim> > > input,
   VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
+	  ExcInternalError());
+  const InternalData &data = static_cast<const InternalData&>(mapping_data);
+  
   switch (mapping_type)
     {
       case mapping_covariant:
-	    transform_covariant(input, 0, output, internal);
-	    return;
-       case mapping_contravariant:
- 	    transform_contravariant(input, 0, output, internal);
- 	    return;
-      default:
-	    Assert(false, ExcNotImplemented());
-    }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingCartesian<dim, spacedim>::transform_covariant (
-  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
-
-  Assert (data.update_flags & update_covariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-  
-				   // simply scale by inverse Jacobian
-				   // (which is diagonal here)
-  for (unsigned int i=0; i<output.size(); ++i)
-    for (unsigned int d=0;d<dim;++d)
-      output[i][d] = input[i+offset][d]/data.length[d];
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingCartesian<dim, spacedim>::transform_covariant (
-  const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
-
-  Assert (data.update_flags & update_covariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-  
-  Assert (output.size() + offset <= input.size(), ExcInternalError());
-  
-				   // simply scale by inverse Jacobian
-				   // (which is diagonal here)
-  for (unsigned int i=0; i<output.size(); ++i)
-    for (unsigned int d=0; d<dim; ++d)
-      for (unsigned int p=0; p<dim; ++p)
-      output[i][d][p] = input[i+offset][d][p]/data.length[p];
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingCartesian<dim, spacedim>::transform_contravariant (
-  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-				   // convert data object to internal
-				   // data for this class. fails with
-				   // an exception if that is not
-				   // possible
-  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
-
-  Assert (data.update_flags & update_contravariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
-  Assert (output.size() + offset <= input.size(), ExcInternalError());
-  
-				   // simply scale by inverse Jacobian
-				   // (which is diagonal here)
-  for (unsigned int i=0; i<output.size(); ++i)
-    for (unsigned int d=0;d<dim;++d)
-      output[i][d] = input[i+offset][d] * data.length[d];
-}
-
-
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d1=0;d1<dim;++d1)
+	    for (unsigned int d2=0;d2<dim;++d2)
+	      output[i][d1][d2] = input[i][d1][d2] / data.length[d2];
+	return;
+      }
+      
+      case mapping_contravariant:
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d1=0;d1<dim;++d1)
+	    for (unsigned int d2=0;d2<dim;++d2)
+	      output[i][d1][d2] = input[i][d1][d2] * data.length[d2];
+	return;
+      }
 
-template<int dim, int spacedim>
-void
-MappingCartesian<dim, spacedim>::transform_contravariant (
-  const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-				   // convert data object to internal
-				   // data for this class. fails with
-				   // an exception if that is not
-				   // possible
-  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
+      case mapping_covariant_gradient:
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d1=0;d1<dim;++d1)
+	    for (unsigned int d2=0;d2<dim;++d2)
+	      output[i][d1][d2] = input[i][d1][d2] / data.length[d2] / data.length[d1];
+	return;
+      }
 
-  Assert (data.update_flags & update_contravariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+      case mapping_contravariant_gradient:
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d1=0;d1<dim;++d1)
+	    for (unsigned int d2=0;d2<dim;++d2)
+	      output[i][d1][d2] = input[i][d1][d2] * data.length[d2] / data.length[d1];
+	return;
+      }
 
-  Assert (output.size() + offset <= input.size(), ExcInternalError());
-  
-				   // simply scale by inverse Jacobian
-				   // (which is diagonal here)
-  for (unsigned int i=0; i<output.size(); ++i)
-    for (unsigned int d=0; d<dim; ++d)
-      for (unsigned int p=0; p<dim; ++p)
-      output[i][d][p] = input[i+offset][d][p] * data.length[p];
+      case mapping_piola:
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_volume_elements,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  for (unsigned int d1=0;d1<dim;++d1)
+	    for (unsigned int d2=0;d2<dim;++d2)
+	      output[i][d1][d2] = input[i][d1][d2] * data.length[d2]
+				  / data.length[d1] / data.volume_element;
+	return;
+      }
+	    
+      default:
+	    Assert(false, ExcNotImplemented());
+    }
 }
 
 
diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc
index 53e84d17b0..22940c6486 100644
--- a/deal.II/deal.II/source/fe/mapping_q.cc
+++ b/deal.II/deal.II/source/fe/mapping_q.cc
@@ -1207,198 +1207,74 @@ void
 MappingQ<dim,spacedim>::transform (
   const VectorSlice<const std::vector<Tensor<1,dim> > > input,
   VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
-  switch (mapping_type)
-    {
-      case mapping_covariant:
-	    transform_covariant(input, 0, output, internal);
-	    return;
-       case mapping_contravariant:
- 	    transform_contravariant(input, 0, output, internal);
- 	    return;
-      default:
-	    Assert(false, ExcNotImplemented());
-    }
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::transform (
-  const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &internal,
-  const MappingType mapping_type) const
-{
-  switch (mapping_type)
-    {
-      case mapping_covariant:
-	    transform_covariant(input, 0, output, internal);
-	    return;
-       case mapping_contravariant:
- 	    transform_contravariant(input, 0, output, internal);
- 	    return;
-      default:
-	    Assert(false, ExcNotImplemented());
-    }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::transform_covariant (
-  const VectorSlice<const std::vector<Tensor<1,spacedim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-  const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
-    dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
-  Assert(q1_data!=0, ExcInternalError());
-
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
-
-  if (q1_data->is_mapping_q1_data)
-    tensor = q1_data->covariant.begin();
-  else
-    {
-      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
-      Assert(data!=0, ExcInternalError());
-
-      if (data->use_mapping_q1_on_current_cell)
-	tensor = data->mapping_q1_data.covariant.begin();
-      else
-	tensor = data->covariant.begin();    
-    }
-  
-  Tensor<1, spacedim> auxiliary;
-  
-  for (unsigned int i=0; i<output.size(); ++i)
-    {
-      for (unsigned int d=0;d<dim;++d)
-	auxiliary[d] = input[i][d];
-      contract (output[i], auxiliary, *(tensor++));
-    }
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::transform_covariant (
-  const VectorSlice<const std::vector<Tensor<2,spacedim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-  const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
-    dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
-  Assert(q1_data!=0, ExcInternalError());
-
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
-
-  if (q1_data->is_mapping_q1_data)
-    tensor = q1_data->covariant.begin();
-  else
-    {
-      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
-      Assert(data!=0, ExcInternalError());
-
-      if (data->use_mapping_q1_on_current_cell)
-	tensor = data->mapping_q1_data.covariant.begin();
-      else
-	tensor = data->covariant.begin();
-    }
-
-  for (unsigned int i=0; i<output.size(); ++i)
-    contract (output[i], input[i+offset], *(tensor++));
-}
-
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::transform_contravariant (
-  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
   AssertDimension (input.size(), output.size());
-  
+				   // The data object may be jsut a
+				   // MappingQ1::InternalData, so we
+				   // have to test for this first.
   const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
     dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
   Assert(q1_data!=0, ExcInternalError());
-  
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
 
-  if (q1_data->is_mapping_q1_data)
-    tensor = q1_data->contravariant.begin();
-  else
-    {
-      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
-      Assert(data!=0, ExcInternalError());
-
-      if (data->use_mapping_q1_on_current_cell)
-	tensor = data->mapping_q1_data.contravariant.begin();
-      else
-	tensor = data->contravariant.begin();    
-    }
-  
-  Tensor<1, spacedim> auxiliary;
-  
-  for (unsigned int i=0; i<output.size(); ++i)
+				   // If it is a genuine
+				   // MappingQ::InternalData, we have
+				   // to test further
+  if (!q1_data->is_mapping_q1_data)
     {
-      for (unsigned int d=0;d<dim;++d)
-	auxiliary[d] = input[i][d];
-      contract (output[i], *(tensor++), auxiliary);
+      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
+	      ExcInternalError());
+      const InternalData &data = static_cast<const InternalData&>(mapping_data);
+				       // If we only use the
+				       // Q1-portion, we have to
+				       // extract that data object
+      if (data.use_mapping_q1_on_current_cell)
+	q1_data = &data.mapping_q1_data;
     }
+				   // Now, q1_data should have the
+				   // right tensors in it and we call
+				   // the base classes transform
+				   // function
+  MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);  
 }
 
 
-
 template<int dim, int spacedim>
 void
-MappingQ<dim,spacedim>::transform_contravariant (
+MappingQ<dim,spacedim>::transform (
   const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-  const unsigned int                 offset,
   VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
 {
-  Assert (offset == 0, ExcInternalError());
   AssertDimension (input.size(), output.size());
-  
+				   // The data object may be jsut a
+				   // MappingQ1::InternalData, so we
+				   // have to test for this first.
   const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
     dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
   Assert(q1_data!=0, ExcInternalError());
-  
-  typename std::vector<Tensor<2,dim> >::const_iterator tensor;
 
-  if (q1_data->is_mapping_q1_data)
-    tensor = q1_data->contravariant.begin();
-  else
+				   // If it is a genuine
+				   // MappingQ::InternalData, we have
+				   // to test further
+  if (!q1_data->is_mapping_q1_data)
     {
-      const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
-      Assert(data!=0, ExcInternalError());
-
-      if (data->use_mapping_q1_on_current_cell)
-	tensor = data->mapping_q1_data.contravariant.begin();
-      else
-	tensor = data->contravariant.begin();    
+      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
+	      ExcInternalError());
+      const InternalData &data = static_cast<const InternalData&>(mapping_data);
+				       // If we only use the
+				       // Q1-portion, we have to
+				       // extract that data object
+      if (data.use_mapping_q1_on_current_cell)
+	q1_data = &data.mapping_q1_data;
     }
-
-  for (unsigned int i=0; i<output.size(); ++i)
-    contract (output[i], *(tensor++), input[i+offset]);
+				   // Now, q1_data should have the
+				   // right tensors in it and we call
+				   // the base classes transform
+				   // function
+  MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);  
 }
 
 
diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc
index d229bda802..740274e8c6 100644
--- a/deal.II/deal.II/source/fe/mapping_q1.cc
+++ b/deal.II/deal.II/source/fe/mapping_q1.cc
@@ -729,15 +729,16 @@ MappingQ1<dim,spacedim>::compute_mapping_support_points(
 
 template<int dim, int spacedim>
 void
-MappingQ1<dim,spacedim>::fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-				const Quadrature<dim>                     &q,
- 			        typename Mapping<dim,spacedim>::InternalDataBase   &mapping_data,
-				std::vector<Point<spacedim> >             &quadrature_points,
-				std::vector<double>                       &JxW_values,
-				std::vector<Tensor<2,spacedim> >          &jacobians,
-				std::vector<Tensor<3,spacedim> >          &jacobian_grads,
-				std::vector<Tensor<2,spacedim> >          &inverse_jacobians,
-				std::vector<Point<spacedim> >             &cell_normal_vectors) const
+MappingQ1<dim,spacedim>::fill_fe_values (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+  const Quadrature<dim>                     &q,
+  typename Mapping<dim,spacedim>::InternalDataBase   &mapping_data,
+  std::vector<Point<spacedim> >             &quadrature_points,
+  std::vector<double>                       &JxW_values,
+  std::vector<Tensor<2,spacedim> >          &jacobians,
+  std::vector<Tensor<3,spacedim> >          &jacobian_grads,
+  std::vector<Tensor<2,spacedim> >          &inverse_jacobians,
+  std::vector<Point<spacedim> >             &cell_normal_vectors) const
 {
   // ensure that the following cast is really correct:
   Assert (dynamic_cast<InternalData *>(&mapping_data) != 0, 
@@ -1184,49 +1185,49 @@ MappingQ1<dim,spacedim>::transform (
   switch (mapping_type)
     {
       case mapping_covariant:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_covariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  {
-		    for (unsigned int d=0;d<dim;++d)
-		      auxiliary[d] = input[i][d];
-		    contract (output[i], auxiliary, data.covariant[i]);
-		  }
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d=0;d<dim;++d)
+	      auxiliary[d] = input[i][d];
+	    contract (output[i], auxiliary, data.covariant[i]);
+	  }
+	return;
+      }
+      
       case mapping_contravariant:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_contravariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  {
-		    for (unsigned int d=0;d<dim;++d)
-		      auxiliary[d] = input[i][d];
-		    contract (output[i], data.contravariant[i], auxiliary);
-		  }
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d=0;d<dim;++d)
+	      auxiliary[d] = input[i][d];
+	    contract (output[i], data.contravariant[i], auxiliary);
+	  }
+	return;
+      }
+      
       case mapping_piola:
-	    if (true)
-	      {
-		Assert (data.update_flags & update_contravariant_transformation,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		Assert (data.update_flags & update_volume_elements,
-			typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-		
-		for (unsigned int i=0; i<output.size(); ++i)
-		  {
-		    for (unsigned int d=0;d<dim;++d)
-		      auxiliary[d] = input[i][d] / data.volume_elements[i];
-		    contract (output[i], data.contravariant[i], auxiliary);
-		  }
-		return;
-	      }
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_volume_elements,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d=0;d<dim;++d)
+	      auxiliary[d] = input[i][d] / data.volume_elements[i];
+	    contract (output[i], data.contravariant[i], auxiliary);
+	  }
+	return;
+      }
+      
       default:
 	    Assert(false, ExcNotImplemented());
     }
@@ -1238,206 +1239,114 @@ void
 MappingQ1<dim,spacedim>::transform (
   const VectorSlice<const std::vector<Tensor<2, dim> > > input,
   VectorSlice<std::vector<Tensor<2, spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
   const MappingType mapping_type) const
 {
   AssertDimension (input.size(), output.size());
-  
-  switch (mapping_type)
-    {
-      case mapping_covariant:
-	    transform_covariant(input, 0, output, internal);
-	    return;
-       case mapping_contravariant:
- 	    transform_contravariant(input, 0, output, internal);
- 	    return;
-      default:
-	    Assert(false, ExcNotImplemented());
-    }
-}
-
-
-// This function is deprecated and has been replaced by transform above
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::transform_covariant (
-  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-  // ensure that the following cast is really correct:
-  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
-	  ExcInternalError());
-  const InternalData &data = static_cast<const InternalData&>(mapping_data);
-
-  Assert (data.update_flags & update_covariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
-  Tensor<1, spacedim> auxiliary;
-  
-  for (unsigned int i=0; i<output.size(); ++i)
-    {
-      for (unsigned int d=0;d<dim;++d)
-	auxiliary[d] = input[i][d];
-      contract (output[i], auxiliary, data.covariant[i]);
-    }
-}
-
-
-
-// This function is deprecated and has been replaced by transform above
-template <int dim, int spacedim>
-void
-MappingQ1<dim, spacedim>::transform_covariant (
-    const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-    const unsigned int                 offset,
-    VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-    const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  Assert (offset == 0, ExcInternalError());
-  AssertDimension (input.size(), output.size());
-  
-				   // ensure that the following cast
-				   // is really correct:
-  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
-	  ExcInternalError());
-  const InternalData &data = static_cast<const InternalData&>(mapping_data);
-
-  Assert (data.update_flags & update_covariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-  
-  for (unsigned int i=0; i<output.size(); ++i)
-    contract (output[i], input[i+offset], data.covariant[i]);
-}
-
-#if deal_II_dimension == 1
-
-
-// This function is deprecated and has been replaced by transform above
-template <>
-void
-MappingQ1<1, 2>::transform_covariant (
-    const VectorSlice<const std::vector<Tensor<2,1> > > ,
-    const unsigned int,
-    VectorSlice<std::vector<Tensor<2,2> > > ,
-    const Mapping<1,2>::InternalDataBase &) const
-{
-    Assert(false, ExcNotImplemented());
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-
-// This function is deprecated and has been replaced by transform above
-template <>
-void
-MappingQ1<2, 3>::transform_covariant (
-    const VectorSlice<const std::vector<Tensor<2,2> > > ,
-    const unsigned int,
-    VectorSlice<std::vector<Tensor<2,3> > > ,
-    const Mapping<2,3>::InternalDataBase &) const
-{
-    Assert(false, ExcNotImplemented());
-}
-
-#endif
-
-
-// This function is deprecated and has been replaced by transform above
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::transform_contravariant (
-  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  // ensure that the following cast is really correct:
   Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
 	  ExcInternalError());
   const InternalData &data = static_cast<const InternalData&>(mapping_data);
 
-  Assert (data.update_flags & update_contravariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
-  Assert (output.size() + offset <= input.size(), ExcInternalError());
-  
-  Tensor<1, spacedim> auxiliary;
+  Tensor<2, spacedim> aux1;
+  Tensor<2, spacedim> aux2;
   
-  for (unsigned int i=0; i<output.size(); ++i)
+  switch (mapping_type)
     {
-      for (unsigned int d=0;d<dim;++d)
-	auxiliary[d] = input[i][d];      
-      contract (output[i], data.contravariant[i], auxiliary);
+      case mapping_covariant:
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d1=0;d1<dim;++d1)
+	      for (unsigned int d2=0;d2<dim;++d2)
+		aux1[d1][d2] = input[i][d1][d2];
+	    contract (output[i], aux1, data.covariant[i]);
+	  }
+	return;
+      }
+      
+      case mapping_contravariant:
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d1=0;d1<dim;++d1)
+	      for (unsigned int d2=0;d2<dim;++d2)
+		aux1[d1][d2] = input[i][d1][d2];
+	    contract (output[i], data.contravariant[i], aux1);
+	  }
+	return;
+      }
+      
+      case mapping_covariant_gradient:
+      {
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d1=0;d1<dim;++d1)
+	      for (unsigned int d2=0;d2<dim;++d2)
+		aux1[d1][d2] = input[i][d1][d2];
+	    
+	    contract(aux2, aux1, data.covariant[i]);
+	    contract(output[i], data.covariant[i], aux2);
+	  }
+	
+	return;
+      }
+      
+      case mapping_contravariant_gradient:
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d1=0;d1<dim;++d1)
+	      for (unsigned int d2=0;d2<dim;++d2)
+		aux1[d1][d2] = input[i][d1][d2];
+	    
+	    contract(aux2, aux1, data.covariant[i]);
+	    contract(output[i], data.contravariant[i], aux2);
+	  }
+	
+	return;
+      }
+      
+      case mapping_piola:
+      {
+	Assert (data.update_flags & update_covariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_contravariant_transformation,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	Assert (data.update_flags & update_volume_elements,
+		typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+	
+	for (unsigned int i=0; i<output.size(); ++i)
+	  {
+	    for (unsigned int d1=0;d1<dim;++d1)
+	      for (unsigned int d2=0;d2<dim;++d2)
+		aux1[d1][d2] = input[i][d1][d2] / data.volume_elements[i];
+	    
+	    contract (aux2, aux1, data.covariant[i]);
+	    contract (output[i], data.contravariant[i], aux2);
+	  }
+	return;
+      }
+      default:
+	    Assert(false, ExcNotImplemented());
     }
 }
 
 
 
-// This function is deprecated and has been replaced by transform above
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::transform_contravariant (
-  const VectorSlice<const std::vector<Tensor<2,dim> > > input,
-  const unsigned int                 offset,
-  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const
-{
-  // ensure that the following cast is really correct:
-  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, 
-	  ExcInternalError());
-  const InternalData &data = static_cast<const InternalData&>(mapping_data);
-
-  Assert (data.update_flags & update_contravariant_transformation,
-	  typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
-  Assert (output.size() + offset <= input.size(), ExcInternalError());
-
-  for (unsigned int i=0; i<output.size(); ++i)
-    contract (output[i], data.contravariant[i], input[i+offset]);
-}
-
-
-#if deal_II_dimension == 1
-
-// This function is deprecated and has been replaced by transform above
-template <>
-void
-MappingQ1<1, 2>::transform_contravariant (
-    const VectorSlice<const std::vector<Tensor<2,1> > > ,
-    const unsigned int,
-    VectorSlice<std::vector<Tensor<2,2> > > ,
-    const Mapping<1,2>::InternalDataBase &) const
-{
-    Assert(false, ExcNotImplemented());
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-// This function is deprecated and has been replaced by transform above
-template <>
-void
-MappingQ1<2, 3>::transform_contravariant (
-    const VectorSlice<const std::vector<Tensor<2,2> > > ,
-    const unsigned int,
-    VectorSlice<std::vector<Tensor<2,3> > > ,
-    const Mapping<2,3>::InternalDataBase &) const
-{
-    Assert(false, ExcNotImplemented());
-}
-
-#endif
-
-
-
 template<int dim, int spacedim>
 Point<spacedim>
 MappingQ1<dim,spacedim>::transform_unit_to_real_cell (
-- 
2.39.5