]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
use new Piola transform for vector valued functions
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Jan 2009 21:20:12 +0000 (21:20 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Jan 2009 21:20:12 +0000 (21:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@18102 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_poly_tensor.h
deal.II/deal.II/include/fe/fe_update_flags.h
deal.II/deal.II/include/fe/mapping.h
deal.II/deal.II/include/fe/mapping_cartesian.h
deal.II/deal.II/source/fe/fe_abf.cc
deal.II/deal.II/source/fe/fe_poly_tensor.cc
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
deal.II/deal.II/source/fe/mapping_cartesian.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index 6e44ccbd1f768aaaf880d56c73eb58ca4a088707..e4dfabb33abfc53bd6006f4465634a4d20ab661d 100644 (file)
@@ -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;
index d854cf729994b0b135a5be5fde5ed8a9c509d7d0..b81bd855ff50c8a684fe25ec803388f2f7fe6488 100644 (file)
@@ -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
 };
 
 
index 95779156ee6b0a8c1e0034fab2c70c0b79b71811..51023ab054494bb86dca46c563369710d1e71e90 100644 (file)
@@ -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
index 6967368e8d8a2f5961192994c95e14c0a31dde98..3056a686c7b488ff644a57c5ebd898ada34c9eb4 100644 (file)
@@ -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
index ce017586d3fbae1aad6a36aea969f9cf82e49d94..55d99e23e33d976591bdb58744c6e705cfa5780f 100644 (file)
@@ -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
index 5b53f9613a13ec897e19074475b686afca20870e..3ca14a3a4ed96520e3f6d444f2590315c0e4283d 100644 (file)
@@ -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;
 
index 6b0b24da4eecd3edc4e576adfb865235fd66a4a7..900ebb1c157f852a6ebbd4f5810117e3515f2be8 100644 (file)
@@ -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)
index bf20b13275c6efe63e2c461ddfdb4a9b927835e9..8af15830ad7b1f5ee9664101a4e83bee460c221e 100644 (file)
@@ -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
index 9613be68526c215c7a2efd729358ce149638f8a4..c9f4dbb9c37a9132ef3622ebbaf33549ebcab4dc 100644 (file)
@@ -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());
     }
index 5c0b52828cc1e54145e7252243504f8b80fe2a36..d229bda80231b86c4383073e62839feb89e36bc0 100644 (file)
@@ -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());
     }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.