From ec2319d145e33a813081b75126eb44be611a4a81 Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 6 Jan 2009 21:20:12 +0000
Subject: [PATCH] use new Piola transform for vector valued functions

git-svn-id: https://svn.dealii.org/trunk@18102 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/include/fe/fe_poly_tensor.h   |  53 +------
 deal.II/deal.II/include/fe/fe_update_flags.h  |  27 ++--
 deal.II/deal.II/include/fe/mapping.h          |  21 ++-
 .../deal.II/include/fe/mapping_cartesian.h    |   5 +
 deal.II/deal.II/source/fe/fe_abf.cc           |   8 +-
 deal.II/deal.II/source/fe/fe_poly_tensor.cc   | 131 ++++++++++--------
 .../deal.II/source/fe/fe_raviart_thomas.cc    |   8 +-
 .../source/fe/fe_raviart_thomas_nodal.cc      |   4 +-
 .../deal.II/source/fe/mapping_cartesian.cc    |  43 ++++--
 deal.II/deal.II/source/fe/mapping_q1.cc       |  23 +++
 10 files changed, 186 insertions(+), 137 deletions(-)

diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h
index 6e44ccbd1f..e4dfabb33a 100644
--- a/deal.II/deal.II/include/fe/fe_poly_tensor.h
+++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 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
@@ -106,7 +106,7 @@ DEAL_II_NAMESPACE_OPEN
  * in #mapping_type. Therefore, each constructor should contain a line
  * like:
  * @verbatim
- * this->mapping_type = this->independent_on_cartesian;
+ * this->mapping_type = this->mapping_none;
  * @endverbatim
  *
  * @see PolynomialsBDM, PolynomialsRaviartThomas
@@ -219,57 +219,10 @@ class FE_PolyTensor : public FiniteElement<dim,spacedim>
     virtual UpdateFlags update_each (const UpdateFlags flags) const;
     
   protected:
-				     /**
-				      * Different options for
-				      * transforming the basis
-				      * functions from the reference
-				      * cell to the actual mesh cell.
-				      *
-				      * Most vector valued elements
-				      * either transform shape
-				      * functions to keep node values
-				      * on edges meaningful. Still, in
-				      * special cases, it may be
-				      * possible to avoid the mapping.
-				      */
-    enum MappingType {
-					   /**
-					    * No mapping has been
-					    * selected, throw an error
-					    * if needed.
-					    */
-	  no_mapping,
-					   /**
-					    * Shape functions do not
-					    * depend on actual mesh
-					    * cell
-					    */
-	  independent,
-					   /**
-					    * Shape functions do not
-					    * depend on actual mesh
-					    * cell. The mapping class
-					    * must be
-					    * MappingCartesian.
-					    */
-	  independent_on_cartesian,
-					   /**
-					    * Shape functions are
-					    * transformed covariant.
-					    */ 
-	  covariant,
-					   /**
-					    * Shape functions are
-					    * transformed
-					    * contravariant.
-					    */
-	  contravariant
-    };
-
 				     /**
 				      * The mapping type to be used to
 				      * map shape functions from the
-				      * reference cell to the msh
+				      * reference cell to the mesh
 				      * cell.
 				      */
     MappingType mapping_type;
diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h
index d854cf7299..b81bd855ff 100644
--- a/deal.II/deal.II/include/fe/fe_update_flags.h
+++ b/deal.II/deal.II/include/fe/fe_update_flags.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -221,7 +221,12 @@ enum UpdateFlags
 					* the Mapping.
 					*/
       update_transformation_gradients     = 0x2000,
-
+				       //! Determinant of the Jacobian
+				       /**
+					* Compute the volume element
+					* in each quadrature point.
+					*/
+      update_volume_elements              = 0x4000,
 				       /**
 					* Compute the JxW values
 					* on faces for the cell mapping
@@ -231,25 +236,25 @@ enum UpdateFlags
 					* is used in conjunction with
 					* H_div subspaces like RT and ABF.
 					*/
-      update_cell_JxW_values              = 0x4000,
+      update_cell_JxW_values              = 0x8000,
 				       /**
 					* Update the location of the
 					* mapped generalized support
 					* points of the element.
 					*/
-      update_support_points               = 0x8000,
+      update_support_points               = 0x10000,
 				       /**
 					* Update the Jacobian of the
 					* mapping in generalized
 					* support points.
 					*/
-      update_support_jacobians            = 0x10000,
+      update_support_jacobians            = 0x20000,
 				       /**
 					* Update the inverse Jacobian
 					* of the mapping in
 					* generalized support points.
 					*/
-      update_support_inverse_jacobians    = 0x20000,
+      update_support_inverse_jacobians    = 0x40000,
 				       /**
 					* @deprecated Update
 					* quadrature points
@@ -266,8 +271,14 @@ enum UpdateFlags
 					* vectors. Use
 					* update_face_normal_vectors
 					*/
-      update_normal_vectors               = update_face_normal_vectors
-
+      update_normal_vectors               = update_face_normal_vectors,
+				       //! Values needed for Piola transform
+				       /**
+					* Combination of the flags
+					* needed for Piola transform
+					* of Hdiv elements.
+					*/
+      update_piola = update_volume_elements | update_contravariant_transformation
 };
 
 
diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h
index 95779156ee..51023ab054 100644
--- a/deal.II/deal.II/include/fe/mapping.h
+++ b/deal.II/deal.II/include/fe/mapping.h
@@ -32,8 +32,6 @@ template <int dim, int spacedim> class FEValues;
 template <int dim, int spacedim> class FEFaceValues;
 template <int dim, int spacedim> class FESubfaceValues;
 
-//TODO: Offset in transform functions should be replaced by initializing VectorSlice correctly
-
                                      /**
                                       * The transformation type used
                                       * for the Mapping::transform() functions.
@@ -55,18 +53,21 @@ template <int dim, int spacedim> class FESubfaceValues;
 				      * for vectors. If such a
 				      * MappingType is applied to a
 				      * rank 2 tensor, it is implied
-				      * that the resulting Tensor
-				      * corresponds to the derivative
-				      * of the vector.
+				      * that the mapping is applied to
+				      * each column.
                                       */
 enum MappingType
 {
+/// No mapping
+      mapping_none = 0x0000,
 /// Covariant mapping
       mapping_covariant = 0x0001,
 /// Contravariant mapping
       mapping_contravariant = 0x0002,
 /// The Piola transform usually used for Hdiv elements
-      mapping_piola = 0x0003,
+      mapping_piola = 0x0100,
+/// The Piola transform for the derivative of an Hdiv element
+      mapping_piola_gradient = 0x0101,
 /// The mapping used for Raviart-Thomas elements
       mapping_raviart_thomas = mapping_piola,
 /// The mapping used for BDM elements
@@ -261,6 +262,14 @@ class Mapping : public Subscriptor
 					  * object.
 					  */
 	virtual unsigned int memory_consumption () const;
+
+					 /**
+					  * The determinant of the
+					  * Jacobian in each
+					  * quadrature point. Filled
+					  * if #update_volume_elements.
+					  */
+	std::vector<double> volume_elements;
 	
 					 /**
 					  * The positions of the
diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h
index 6967368e8d..3056a686c7 100644
--- a/deal.II/deal.II/include/fe/mapping_cartesian.h
+++ b/deal.II/deal.II/include/fe/mapping_cartesian.h
@@ -187,6 +187,11 @@ class MappingCartesian : public Mapping<dim,spacedim>
 					  */
 	Tensor<1,dim> length;
 
+					 /**
+					  * The volume element
+					  */
+	double volume_element;
+	
 					 /**
 					  * Vector of all quadrature
 					  * points. Especially, all
diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc
index ce017586d3..55d99e23e3 100644
--- a/deal.II/deal.II/source/fe/fe_abf.cc
+++ b/deal.II/deal.II/source/fe/fe_abf.cc
@@ -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
@@ -56,7 +56,7 @@ FE_ABF<dim>::FE_ABF (const unsigned int deg)
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
   
-  this->mapping_type = this->contravariant;
+  this->mapping_type = mapping_piola;
 				   // First, initialize the
 				   // generalized support points and
 				   // quadrature weights, since they
@@ -539,12 +539,12 @@ FE_ABF<dim>::update_each (const UpdateFlags flags) const
 
   if (flags & update_values)
     out |= update_values             | update_covariant_transformation
-                                     | update_contravariant_transformation 
+                                     | update_piola
                                      | update_cell_JxW_values
                                      | update_JxW_values;
   if (flags & update_gradients)
     out |= update_gradients          | update_covariant_transformation 
-                                     | update_contravariant_transformation
+                                     | update_piola
                                      | update_cell_JxW_values
                                      | update_JxW_values;
   //TODO: Set update flags appropriately and figure out, how the second
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 5b53f9613a..3ca14a3a4e 100644
--- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc
+++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc
@@ -257,9 +257,10 @@ FE_PolyTensor<POLY,dim,spacedim>::shape_grad_grad_component (const unsigned int
 
 template <class POLY, int dim, int spacedim>
 typename Mapping<dim,spacedim>::InternalDataBase *
-FE_PolyTensor<POLY,dim,spacedim>::get_data (const UpdateFlags      update_flags,
-			     const Mapping<dim,spacedim>    &mapping,
-			     const Quadrature<dim> &quadrature) const
+FE_PolyTensor<POLY,dim,spacedim>::get_data (
+  const UpdateFlags      update_flags,
+  const Mapping<dim,spacedim>    &mapping,
+  const Quadrature<dim> &quadrature) const
 {
 				   // generate a new data object and
 				   // initialize some fields
@@ -371,11 +372,6 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 				   // possible
   InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
 
-//   Assert(mapping_type == independent
-// 	 || ( mapping_type == independent_on_cartesian
-// 	      && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
-// 	 ExcNotImplemented());
-
   const unsigned int n_q_points = quadrature.size();
   const UpdateFlags flags(fe_data.current_update_flags());
   
@@ -400,8 +396,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 	{
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
@@ -409,13 +404,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 		break;
 	      }
 	      
-	      case covariant:
-	      case contravariant:
+	      case mapping_covariant:
+	      case mapping_contravariant:
 	      {
 						 // Use auxiliary vector for
 						 // transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == covariant)
+		if (mapping_type == mapping_covariant)
 		  mapping.transform(fe_data.shape_values[i],
 				    shape_values, mapping_data, mapping_covariant);
 		else
@@ -425,9 +420,12 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 						 // then copy over to target:
 		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 == contravariant ?
+		      J = (mapping_type == mapping_contravariant ?
 			   data.JxW_values[k] / quadrature.weight(k)
 			   :
 			   1.0);
@@ -438,6 +436,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 		  }
 		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;
+		      }
 	      
 	      default:
 		Assert(false, ExcNotImplemented());
@@ -451,8 +461,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		mapping.transform(fe_data.shape_grads[i], shape_grads1,
 				  mapping_data, mapping_covariant);
@@ -462,7 +471,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 		break;
 	      }
 	      
-	      case covariant:
+	      case mapping_covariant:
 	      {
 		mapping.transform(fe_data.shape_grads[i], shape_grads1,
 				  mapping_data, mapping_covariant);
@@ -481,8 +490,8 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 
 		break;
 	      }
-	      
-	      case contravariant:
+	      case mapping_piola:
+	      case mapping_contravariant:
 	      {
 		mapping.transform(fe_data.shape_grads[i], shape_grads1,
 				  mapping_data, mapping_covariant);
@@ -494,11 +503,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
 		  for (unsigned int d=0; d<dim; ++d)
 		    {
 						       // Recompute determinant
-		      const double
-			J = (mapping_type == contravariant ?
-			     data.JxW_values[k] / quadrature.weight(k)
-			     :
-			     1.0);
+		      const double J = data.JxW_values[k] / quadrature.weight(k);
 		      data.shape_gradients[first+d][k] = sign_change[i] * 
 			shape_grads2[k][d] / J;
 		    }
@@ -574,8 +579,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 	{
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
@@ -583,13 +587,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 		break;
 	      }
 	      
-	      case covariant:
-	      case contravariant:
+	      case mapping_covariant:
+	      case mapping_contravariant:
 	      {
 						 // Use auxiliary vector
 						 // for transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == covariant)
+		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
@@ -611,7 +615,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 						     // face in
 						     // fill_fe_face_values
 		    const double
-		      J = (mapping_type == contravariant ?
+		      J = (mapping_type == mapping_contravariant ?
 			   data.cell_JxW_values[k] / quadrature.weight(k)
 			   :
 			   1.0);
@@ -623,6 +627,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 		
 		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;
+		      }
 	      
 	      default:
 		Assert(false, ExcNotImplemented());
@@ -636,8 +652,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
 				  shape_grads1, mapping_data, mapping_covariant);
@@ -647,7 +662,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 		break;
 	      }
 	      
-	      case covariant:
+	      case mapping_covariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
 				  shape_grads1, mapping_data, mapping_covariant);
@@ -665,8 +680,8 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
 		break;
 	      }
-	      
-	      case contravariant:
+	      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);
@@ -680,8 +695,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
 						       // recompute
 						       // determinant in the
 						       // same way as above
-		      const double
-			J = data.cell_JxW_values[k] / quadrature.weight(k);
+		      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;
 		    }
@@ -755,8 +769,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 	{
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		for (unsigned int k=0; k<n_q_points; ++k)
 		  for (unsigned int d=0; d<dim; ++d)
@@ -764,13 +777,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 		break;
 	      }
 	      
-	      case covariant:
-	      case contravariant:
+	      case mapping_covariant:
+	      case mapping_contravariant:
 	      {
 						 // Use auxiliary vector for
 						 // transformation
 		std::vector<Tensor<1,dim> > shape_values (n_q_points);
-		if (mapping_type == covariant)
+		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
@@ -792,7 +805,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 						     // face in
 						     // fill_fe_face_values
 		    const double
-		      J = (mapping_type == contravariant ?
+		      J = (mapping_type == mapping_contravariant ?
 			   data.cell_JxW_values[k] / quadrature.weight(k)
 			   :
 			   1.0);
@@ -804,7 +817,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 
 		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;
+		      }
 	      default:
 		Assert(false, ExcNotImplemented());
 	    }
@@ -817,8 +841,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 
 	  switch (mapping_type)
 	    {
-	      case independent:
-	      case independent_on_cartesian:
+	      case mapping_none:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
 				  shape_grads1, mapping_data, mapping_covariant);
@@ -828,7 +851,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 		break;
 	      }
 	      
-	      case covariant:
+	      case mapping_covariant:
 	      {
 		mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
 				  shape_grads1, mapping_data, mapping_covariant);
@@ -846,8 +869,9 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 		    data.shape_gradients[first+d][k] = shape_grads1[k][d];
 		break;
 	      }
-	      
-	      case contravariant:
+
+	      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);
@@ -861,8 +885,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 						       // recompute
 						       // determinant in the
 						       // same way as above
-		      const double
-			J = data.cell_JxW_values[k] / quadrature.weight(k);
+		      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;
 		    }
@@ -913,8 +936,7 @@ template <class POLY, int dim, int spacedim>
 UpdateFlags
 FE_PolyTensor<POLY,dim,spacedim>::update_once (const UpdateFlags flags) const
 {
-  Assert (mapping_type != no_mapping, ExcNotInitialized());
-  const bool values_once = (mapping_type == independent);
+  const bool values_once = (mapping_type == mapping_none);
   
   UpdateFlags out = update_default;
 
@@ -929,8 +951,7 @@ template <class POLY, int dim, int spacedim>
 UpdateFlags
 FE_PolyTensor<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const
 {
-  Assert (mapping_type != no_mapping, ExcNotInitialized());
-  const bool values_once = (mapping_type == independent);
+  const bool values_once = (mapping_type == mapping_none);
   
   UpdateFlags out = update_default;
 
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 6b0b24da4e..900ebb1c15 100644
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
@@ -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
@@ -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 = this->contravariant;
+  this->mapping_type = mapping_piola;
 				   // First, initialize the
 				   // generalized support points and
 				   // quadrature weights, since they
@@ -464,12 +464,12 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
 
   if (flags & update_values)
     out |= update_values             | update_covariant_transformation
-                                     | update_contravariant_transformation 
+                                     | update_piola 
                                      | update_cell_JxW_values
                                      | update_JxW_values;
   if (flags & update_gradients)
     out |= update_gradients          | update_covariant_transformation 
-                                     | update_contravariant_transformation
+                                     | update_piola
                                      | update_cell_JxW_values
                                      | update_JxW_values;
   if (flags & update_hessians)
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 bf20b13275..8af15830ad 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
@@ -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
@@ -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 = this->independent_on_cartesian;
+  this->mapping_type = mapping_none;
 				   // These must be done first, since
 				   // they change the evaluation of
 				   // basis functions
diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc
index 9613be6852..c9f4dbb9c3 100644
--- a/deal.II/deal.II/source/fe/mapping_cartesian.cc
+++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc
@@ -86,7 +86,7 @@ MappingCartesian<dim, spacedim>::get_data (const UpdateFlags      update_flags,
   data->update_once = update_once(update_flags);
   data->update_each = update_each(update_flags);
   data->update_flags = data->update_once | data->update_each;
-
+  
   return data;
 }
 
@@ -103,7 +103,7 @@ MappingCartesian<dim, spacedim>::get_face_data (const UpdateFlags update_flags,
   data->update_once = update_once(update_flags);
   data->update_each = update_each(update_flags);
   data->update_flags = data->update_once | data->update_each;
-
+  
   return data;
 }
 
@@ -120,7 +120,7 @@ MappingCartesian<dim, spacedim>::get_subface_data (const UpdateFlags update_flag
   data->update_once = update_once(update_flags);
   data->update_each = update_each(update_flags);
   data->update_flags = data->update_once | data->update_each;
-
+  
   return data;
 }
 
@@ -325,13 +325,15 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator& cell,
 				   // equal and are the product of the
 				   // local lengths in each coordinate
 				   // direction
-  if (data.current_update_flags() & update_JxW_values)
+  if (data.current_update_flags() & (update_JxW_values | update_volume_elements))
     {
       double J = data.length[0];
       for (unsigned int d=1;d<dim;++d)
 	J *= data.length[d];
-      for (unsigned int i=0; i<JxW_values.size();++i)
-	JxW_values[i] = J * q.weight(i);
+      data.volume_element = J;
+      if (data.current_update_flags() & update_JxW_values)
+	for (unsigned int i=0; i<JxW_values.size();++i)
+	  JxW_values[i] = J * q.weight(i);
     }
 				   // "compute" Jacobian at the quadrature
 				   // points, which are all the same
@@ -402,7 +404,13 @@ MappingCartesian<dim, spacedim>::fill_fe_face_values (const typename Triangulati
     for (unsigned int i=0; i<boundary_forms.size();++i)
       boundary_forms[i] = J * normal_vectors[i];
 
-  //TODO: Implement cell_JxW_values
+  if (data.current_update_flags() & update_volume_elements)
+    {
+      J = data.length[0];
+      for (unsigned int d=1;d<dim;++d)
+	J *= data.length[d];
+      data.volume_element = J;
+    }
 }
 
 
@@ -464,7 +472,13 @@ MappingCartesian<dim, spacedim>::fill_fe_subface_values (const typename Triangul
     for (unsigned int i=0; i<boundary_forms.size();++i)
       boundary_forms[i] = J * normal_vectors[i];
 
-  //TODO: Implement cell_JxW_values
+  if (data.current_update_flags() & update_volume_elements)
+    {
+      J = data.length[0];
+      for (unsigned int d=1;d<dim;++d)
+	J *= data.length[d];
+      data.volume_element = J;
+    }
 }
 
 
@@ -541,6 +555,19 @@ MappingCartesian<dim,spacedim>::transform (
 		    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;
+	      }
       default:
 	    Assert(false, ExcNotImplemented());
     }
diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc
index 5c0b52828c..d229bda802 100644
--- a/deal.II/deal.II/source/fe/mapping_q1.cc
+++ b/deal.II/deal.II/source/fe/mapping_q1.cc
@@ -393,6 +393,7 @@ MappingQ1<dim,spacedim>::update_each (const UpdateFlags in) const
 				      | update_cell_JxW_values
 				      | update_boundary_forms
 				      | update_normal_vectors
+				      | update_volume_elements
 				      | update_jacobians
 				      | update_jacobian_grads
 				      | update_inverse_jacobians));
@@ -478,6 +479,9 @@ MappingQ1<dim,spacedim>::compute_data (const UpdateFlags      update_flags,
   if (flags & update_contravariant_transformation)
     data.contravariant.resize(n_original_q_points);
 
+  if (flags & update_volume_elements)
+    data.volume_elements.resize(n_original_q_points);
+
   if (flags & update_jacobian_grads)
     data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
   
@@ -675,6 +679,9 @@ MappingQ1<dim,spacedim>::compute_fill (const typename Triangulation<dim,spacedim
 	      data.contravariant[point][i][j] += data_derv[j] * supp_pts[i];
 	}
 
+  if (update_flags & update_volume_elements)
+    for (unsigned int point=0; point<n_q_points; ++point)
+      data.volume_elements[point] = determinant(data.contravariant[point]);
 
 				    // invert contravariant for
 			            // covariant transformation
@@ -1204,6 +1211,22 @@ MappingQ1<dim,spacedim>::transform (
 		  }
 		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;
+	      }
       default:
 	    Assert(false, ExcNotImplemented());
     }
-- 
2.39.5