]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extended the mapping classes with functionality to provide the
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 22 Aug 2006 17:39:31 +0000 (17:39 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 22 Aug 2006 17:39:31 +0000 (17:39 +0000)
determinant of the Jacobian of the cell mapping. In particular this
functionality is required for the fe_poly_tensor class, where a Piola
transformation is used.
Furthermore the procedure to determine required changes of the sign of
the Face DoFs of RT and ABF elements was extended to be able to cope
with refined grids.
The last patch implements the third (difficult) case for computing the
interface contraints in the hp case.

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

17 files changed:
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/fe_values.h
deal.II/deal.II/include/fe/mapping.h
deal.II/deal.II/include/fe/mapping_cartesian.h
deal.II/deal.II/include/fe/mapping_q.h
deal.II/deal.II/include/fe/mapping_q1.h
deal.II/deal.II/include/grid/tria_accessor.h
deal.II/deal.II/source/dofs/dof_tools.cc
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_system.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/fe/mapping_cartesian.cc
deal.II/deal.II/source/fe/mapping_q.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index efd05294a1b75f96782b47c3afc517f95857d33f..577cd571e886cc5c2c561fd952840bff5c135fe3 100644 (file)
@@ -261,6 +261,22 @@ class FE_PolyTensor : public FiniteElement<dim>
          contravariant
     };
 
+
+
+                                    /**
+                                     * On noncartesian grids, the
+                                     * sign of the DoFs associated with
+                                     * the faces of the elements has
+                                     * to be changed in some cases.
+                                     * This procedure implements an
+                                     * algorithm, which determines
+                                     * the DoFs, which need this
+                                     * sign change for a given cell.
+                                     */
+    void get_face_sign_change (const typename Triangulation<dim>::cell_iterator &cell,
+                              std::vector<double> &face_sign) const;
+
+
                                     /**
                                      * The mapping type to be used to
                                      * map shape functions from the
index 5bec6b407fe2a71b7a2969737e3e3d54e980f648..a83e7634ea5e7f1c2f15e44a4edc3f171bcf1e99 100644 (file)
@@ -199,7 +199,19 @@ enum UpdateFlags
                                        * transformation defined by
                                        * the Mapping.
                                        */
-      update_transformation_gradients     = 0x1000
+      update_transformation_gradients     = 0x1000,
+
+                                      /**
+                                       * Compute the JxW values
+                                       * on faces for the cell mapping
+                                       * and not the face mapping!
+                                       * This functionality is required
+                                       * for the Piola mapping which
+                                       * is used in conjunction with
+                                       * H_div subspaces like RT and ABF.
+                                       */
+      update_cell_JxW_values              = 0x2000
+
 };
 
 
index 0b6e2d5a7facc5870d89ce7962ccfe0b3ec56ff3..5941ecd692c1eeda082210ff9bfb4027c00a0588 100644 (file)
@@ -183,6 +183,31 @@ class FEValuesData
                                      */
     std::vector<double>       JxW_values;
 
+                                    /**
+                                     * Store an array of weights
+                                     * times the Jacobi determinant
+                                     * at the quadrature points. This
+                                     * function is reset each time
+                                     * reinit() is called. The
+                                     * Jacobi determinant is actually
+                                     * the reciprocal value of the
+                                     * Jacobi matrices stored in this
+                                     * class, see the general
+                                     * documentation of this class
+                                     * for more information.
+                                     *
+                                     * In contrast to JxW_values,
+                                     * this vector uses the transformation
+                                     * of the cell even in case the
+                                     * corresponding FEValues object
+                                     * refers to a face. Currently the
+                                     * only purpose is to provide
+                                     * exactly this value to the
+                                     * transformation which is used
+                                     * for the PolyTensor elements.
+                                     */
+    std::vector<double>       cell_JxW_values;
+
                                     /**
                                      * Array of quadrature points. This array
                                      * is set up upon calling reinit() and
index e4bfa3aaca1062381e0e822e28a755a8b76c6062..bb0660c0e3c734d3e49ecb74c6087e084083f47d 100644 (file)
@@ -472,7 +472,8 @@ class Mapping : public Subscriptor
                         std::vector<Point<dim> >        &quadrature_points,
                         std::vector<double>                      &JxW_values,
                         std::vector<Tensor<1,dim> >     &boundary_form,
-                        std::vector<Point<dim> >        &normal_vectors) const = 0;
+                        std::vector<Point<dim> >        &normal_vectors,
+                        std::vector<double>             &cell_JxW_values) const = 0;
 
                                     /**
                                      * See above.
@@ -486,7 +487,8 @@ class Mapping : public Subscriptor
                            std::vector<Point<dim> >        &quadrature_points,
                            std::vector<double>                      &JxW_values,
                            std::vector<Tensor<1,dim> >     &boundary_form,
-                           std::vector<Point<dim> >        &normal_vectors) const = 0;
+                           std::vector<Point<dim> >        &normal_vectors,
+                           std::vector<double>             &cell_JxW_values) const = 0;
 
                                     /**
                                      * Give class @p FEValues access
index 804ab76a9e5ae6613ed522f46e19774711db7945..967b58c862dcb60d1c0dc9ad1de585e1edaf6c84 100644 (file)
@@ -68,7 +68,8 @@ class MappingCartesian : public Mapping<dim>
                         std::vector<Point<dim> >        &quadrature_points,
                         std::vector<double>             &JxW_values,
                         std::vector<Tensor<1,dim> >        &boundary_form,
-                        std::vector<Point<dim> >        &normal_vectors) const ;
+                        std::vector<Point<dim> >        &normal_vectors,
+                        std::vector<double>             &cell_JxW_values) const ;
     virtual void
     fill_fe_subface_values (const typename Triangulation<dim>::cell_iterator &cell,
                            const unsigned int face_no,
@@ -78,7 +79,8 @@ class MappingCartesian : public Mapping<dim>
                            std::vector<Point<dim> >        &quadrature_points,
                            std::vector<double>             &JxW_values,
                            std::vector<Tensor<1,dim> >        &boundary_form,
-                           std::vector<Point<dim> >        &normal_vectors) const ;
+                           std::vector<Point<dim> >        &normal_vectors,
+                           std::vector<double>             &cell_JxW_values) const ;
 
                                     /**
                                      * Implementation of the
@@ -232,7 +234,8 @@ template <> void MappingCartesian<1>::fill_fe_face_values (
   std::vector<Point<1> >&,
   std::vector<double>&,
   std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const;
+  std::vector<Point<1> >&,
+  std::vector<double>&) const;
 
 template <> void MappingCartesian<1>::fill_fe_subface_values (
   const Triangulation<1>::cell_iterator &,
@@ -243,7 +246,8 @@ template <> void MappingCartesian<1>::fill_fe_subface_values (
   std::vector<Point<1> >&,
   std::vector<double>&,
   std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const;
+  std::vector<Point<1> >&,
+  std::vector<double>&) const;
 
 #endif // DOXYGEN  
 
index f6f4d01ba8526b1d15b193211bfbf4b682811863..a98a841d533ec93bb1c3d4015a06718a9bccf13f 100644 (file)
@@ -220,7 +220,8 @@ class MappingQ : public MappingQ1<dim>
                         typename std::vector<Point<dim> >        &quadrature_points,
                         std::vector<double>             &JxW_values,
                         typename std::vector<Tensor<1,dim> >        &exterior_form,
-                        typename std::vector<Point<dim> >        &normal_vectors) const ;
+                        typename std::vector<Point<dim> >        &normal_vectors,
+                        std::vector<double>             &cell_JxW_values) const ;
 
                                     /**
                                      * Implementation of the interface in
@@ -235,7 +236,8 @@ class MappingQ : public MappingQ1<dim>
                            typename std::vector<Point<dim> >        &quadrature_points,
                            std::vector<double>             &JxW_values,
                            typename std::vector<Tensor<1,dim> >        &exterior_form,
-                           typename std::vector<Point<dim> >        &normal_vectors) const ;
+                           typename std::vector<Point<dim> >        &normal_vectors,
+                           std::vector<double>             &cell_JxW_values) const ;
 
                                     /**
                                      * For <tt>dim=2,3</tt>. Append the
index 2a9fa19e080191154eb67e1683f90328a2b63df2..415e236fa4a94e664e487ac98901a2b4698f21d7 100644 (file)
@@ -397,7 +397,8 @@ class MappingQ1 : public Mapping<dim>
                         typename std::vector<Point<dim> >        &quadrature_points,
                         std::vector<double>             &JxW_values,
                         typename std::vector<Tensor<1,dim> >        &boundary_form,
-                        typename std::vector<Point<dim> >        &normal_vectors) const ;
+                        typename std::vector<Point<dim> >        &normal_vectors,
+                        std::vector<double>             &cell_JxW_values) const ;
 
                                     /**
                                      * Implementation of the interface in
@@ -412,7 +413,8 @@ class MappingQ1 : public Mapping<dim>
                            typename std::vector<Point<dim> >        &quadrature_points,
                            std::vector<double>             &JxW_values,
                            typename std::vector<Tensor<1,dim> >        &boundary_form,
-                           typename std::vector<Point<dim> >        &normal_vectors) const ;
+                           typename std::vector<Point<dim> >        &normal_vectors,
+                           std::vector<double>             &cell_JxW_values) const ;
     
                                     /**
                                      * Compute shape values and/or
@@ -483,7 +485,8 @@ class MappingQ1 : public Mapping<dim>
                            std::vector<Point<dim> >    &quadrature_points,
                            std::vector<double>         &JxW_values,
                            std::vector<Tensor<1,dim> > &boundary_form,
-                           std::vector<Point<dim> > &normal_vectors) const;
+                           std::vector<Point<dim> > &normal_vectors,
+                           std::vector<double>         &cell_JxW_values) const;
     
                                     /**
                                      * Compute shape values and/or
@@ -712,7 +715,8 @@ template <> void MappingQ1<1>::compute_fill_face (
   std::vector<Point<1> > &,
   std::vector<double> &,
   std::vector<Tensor<1,1> > &,
-  std::vector<Point<1> > &) const;
+  std::vector<Point<1> > &,
+  std::vector<double> &) const;
 
 template <> void MappingQ1<1>::fill_fe_face_values (
   const Triangulation<1>::cell_iterator &,
@@ -722,7 +726,8 @@ template <> void MappingQ1<1>::fill_fe_face_values (
   std::vector<Point<1> >&,
   std::vector<double>&,
   std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const;
+  std::vector<Point<1> >&,
+  std::vector<double>&) const;
 
 template <> void MappingQ1<1>::fill_fe_subface_values (
   const Triangulation<1>::cell_iterator &,
@@ -733,7 +738,8 @@ template <> void MappingQ1<1>::fill_fe_subface_values (
   std::vector<Point<1> >&,
   std::vector<double>&,
   std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const;
+  std::vector<Point<1> >&,
+  std::vector<double>&) const;
 
 #endif // DOXYGEN
 
index 06a8834568cbe2bec4fed85b84d5375c019d7f69..de217b11d8d9ecdc39c1aee125d09bc60a130467 100644 (file)
@@ -2538,6 +2538,11 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                      *  @p ith neighbor.  If the
                                      *  neighbor does not exist, an
                                      *  invalid iterator is returned.
+                                     *  
+                                     *  <bf>Note</bf> (cf. TriaLevel<0>):
+                                     *  The neighbor of a cell has at most the
+                                     *  same level as this cell, i.e. it may
+                                     *  or may not be refined.
                                      */
     TriaIterator<dim,CellAccessor<dim> >
     neighbor (const unsigned int i) const;
index 1e3b832cf74a491897502ac9555d63b64fdae254..d7a231fb65a2112ec267ee02a5157b25d001ad62 100644 (file)
@@ -2287,13 +2287,29 @@ namespace internal
                                                 // that the other fine elements will
                                                 // be constrained to the coarse element.
                {
-                 Assert (false, ExcNotImplemented ());           
+                 //              Assert (false, ExcNotImplemented ());
 
                  typename DH::active_cell_iterator neighbor_child
                    = cell->neighbor_child_on_subface (face, min_degree_subface);
                  const unsigned int n_dofs_on_children = neighbor_child->get_fe().dofs_per_face;
                  dofs_on_children.resize (n_dofs_on_children);
 
+
+                                                      // Get DoFs on child cell with
+                                                      // lowest polynomial degree.
+                                                       // All other DoFs will be constrained
+                                                       // to the DoFs of this face.               
+                 const unsigned int subface_fe_index
+                   = neighbor_child->active_fe_index();
+                     
+                                                      // Same procedure as for the
+                                                      // mother cell. Extract the face
+                                                      // DoFs from the cell DoFs.
+                 cell->face(face)->child(min_degree_subface)
+                   ->get_dof_indices (dofs_on_children,
+                                      subface_fe_index);
+
+
                                                   // The idea is to introduce
                                                   // a "virtual" intermediate coarse
                                                   // level face with the lowest
@@ -2361,12 +2377,88 @@ namespace internal
                  
                                                   // Now create the element
                                                   // constraint for this subface.
-                 FullMatrix<double> face_constraints_m (n_dofs_on_children,
-                                                        n_dofs_on_children);
+                 FullMatrix<double> fc_sface_ipol (n_dofs_on_children,
+                                                   n_dofs_on_children);
+                 FullMatrix<double> fc_ipol_sface (n_dofs_on_children,
+                                                   n_dofs_on_children);
                  neighbor_child->get_fe().get_subface_interpolation_matrix (neighbor_child->get_fe (),
                                                                             min_degree_subface,
-                                                                            face_constraints_m);
-//TODO: Continue ...
+                                                                            fc_sface_ipol); 
+                 // Invert it, to get a mapping from the DoFs of the
+                 // "Master-subface" to the intermediate layer.
+                 fc_ipol_sface.invert (fc_sface_ipol);
+
+                 // Create constraint matrix for the mother face.
+                 const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face;
+                 dofs_on_mother.resize (n_dofs_on_mother);
+                 cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ());
+
+                 FullMatrix<double> fc_mother_ipol (n_dofs_on_children,
+                                                    n_dofs_on_mother);
+                 FullMatrix<double> fc_mother_sface (n_dofs_on_children,
+                                                     n_dofs_on_mother);
+                 neighbor_child->get_fe ().get_face_interpolation_matrix (cell->get_fe(),
+                                                                          fc_mother_ipol);
+                 fc_ipol_sface.mmult (fc_mother_sface, fc_mother_ipol);
+
+                 // Add constraints to global constraint
+                 // matrix.
+                 filter_constraints (dofs_on_children,
+                                     dofs_on_mother,
+                                     fc_mother_sface,
+                                     constraints);
+
+                                                  // Now create constraint matrices for
+                                                  // the subfaces and assemble them
+                 for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+                   {
+                     // As the "Master-subface" does not need constraints, skip it.
+                     if (c != min_degree_subface)
+                       {
+                         typename DH::active_cell_iterator neighbor_child_slave
+                           = cell->neighbor_child_on_subface (face, c);
+                         const unsigned int n_dofs_on_mother = neighbor_child_slave->get_fe().dofs_per_face;
+                         dofs_on_mother.resize (n_dofs_on_mother);
+
+                                                      // Find face number on the finer
+                                                      // neighboring cell, which is
+                                                      // shared the face with the
+                                                      // face of the coarser cell.
+                         const unsigned int neighbor2=
+                           cell->neighbor_of_neighbor(face);
+                         Assert (neighbor_child_slave->face(neighbor2) == cell->face(face)->child(c),
+                                 ExcInternalError());
+                     
+                                                      // Same procedure as for the
+                                                      // mother cell. Extract the face
+                                                      // DoFs from the cell DoFs.
+                         const unsigned int subface_fe_index
+                           = neighbor_child_slave->active_fe_index();
+                     
+                         cell->face(face)->child(c)
+                           ->get_dof_indices (dofs_on_mother,
+                                              subface_fe_index);                     
+                                             
+                                                      // Now create the element
+                                                      // constraint for this subface.
+                         FullMatrix<double> fc_child_sface_ipol (n_dofs_on_children,
+                                                                 n_dofs_on_mother);
+                         FullMatrix<double> fc_child_sface_sface (n_dofs_on_children,
+                                                                  n_dofs_on_mother);
+                         neighbor_child->get_fe ().get_subface_interpolation_matrix 
+                           (neighbor_child_slave->get_fe(),
+                            c, fc_child_sface_ipol);
+
+                         fc_ipol_sface.mmult (fc_child_sface_sface, fc_child_sface_ipol);
+
+                                                      // Add constraints to global constraint
+                                                      // matrix.
+                         filter_constraints (dofs_on_children,
+                                             dofs_on_mother,
+                                             fc_child_sface_sface,
+                                             constraints);
+                       }
+                   }
                }
            }
          else
index c16dc5a6fa1c4be85d1d3ff109b5c5d46fb9750c..ea3a6929177809c0422baac97f59f7ee90b6b9a8 100644 (file)
@@ -62,9 +62,6 @@ FE_ABF<dim>::FE_ABF (const unsigned int deg)
   FullMatrix<double> M(n_dofs, n_dofs);
   FETools::compute_node_matrix(M, *this);
 
-  //TODO: Remove debugging output
-  //  M.print (std::cout);
-
   this->inverse_node_matrix.reinit(n_dofs, n_dofs);
   this->inverse_node_matrix.invert(M);
                                   // From now on, the shape functions
@@ -532,10 +529,12 @@ FE_ABF<dim>::update_each (const UpdateFlags flags) const
   if (flags & update_values)
     out |= update_values             | update_covariant_transformation
                                      | update_contravariant_transformation 
+                                     | update_cell_JxW_values
                                      | update_JxW_values;
   if (flags & update_gradients)
     out |= update_gradients          | update_covariant_transformation 
                                      | update_contravariant_transformation
+                                     | update_cell_JxW_values
                                      | update_JxW_values;
   //TODO: Set update flags appropriately and figure out, how the second
   // derivatives for the RT elements can be computed correctly.
index fa8812ced23ce13a8cb363df7db889c3374e23f1..cf97b95d1e0e4601542d64c06f64a6edbe1dfb56 100644 (file)
@@ -154,6 +154,76 @@ FE_PolyTensor<POLY,dim>::shape_grad_grad_component (
 
 
 
+//---------------------------------------------------------------------------
+// Utility method, which is used to determine the change of sign for
+// the DoFs on the faces of the given cell.
+//---------------------------------------------------------------------------
+
+template <class POLY, int dim>
+void
+FE_PolyTensor<POLY,dim>::get_face_sign_change (
+     const typename Triangulation<dim>::cell_iterator &cell,
+     std::vector<double> &face_sign) const
+{
+  Assert (face_sign.size () == this->dofs_per_cell,
+         ExcDimensionMismatch (face_sign.size (), this->dofs_per_cell));
+
+  // Default is no sign change. I.e. multiply by one.
+  std::fill (face_sign.begin (), face_sign.end (), 1.0);
+
+#if deal_II_dimension > 1
+  const unsigned int dofs_per_face = this->dofs_per_face;
+
+  if (dim == 2)
+    {
+      for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2; 
+          f < GeometryInfo<dim>::faces_per_cell; ++f)
+       {
+         typename Triangulation<dim>::face_iterator face = cell->face (f);
+         if (!face->at_boundary ())
+           {
+             const unsigned int neighbor_level = cell->neighbor (f)->level ();
+             const unsigned int cell_level = cell->level ();
+
+             unsigned int nn = (unsigned int) -1;
+
+             // Same level, the easy case.
+             if (neighbor_level == cell_level)
+               nn = cell->neighbor_of_neighbor (f);
+             else
+               {
+                 // Neighbor is more refined (due to the internal restrictions,
+                 // this case should never be encountered!
+                 if (neighbor_level > cell_level)
+                   {
+                     Assert (false, 
+                             ExcMessage ("neighbor_level is larger than cell level!"));
+                   }
+                 // Neighbor is less refined
+                 else
+                   nn = cell->neighbor_of_coarser_neighbor (f).first;
+               }
+             
+             Assert (nn != (unsigned int) -1,
+                     ExcInternalError ());                   
+             
+             if (nn < GeometryInfo<dim>::faces_per_cell / 2)
+               {
+                 for (unsigned int j = 0; j < dofs_per_face; ++j)
+                   face_sign[f * dofs_per_face + j] = -1.0;
+               }
+           }
+       }
+    }
+  else
+    {
+      // TODO: Think about 3D!.
+    }
+#endif
+
+}  
+
+
 //---------------------------------------------------------------------------
 // Data field initialization
 //---------------------------------------------------------------------------
@@ -289,32 +359,7 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
   // Compute eventual sign changes depending on the neighborhood
   // between two faces.
   std::vector<double> sign_change (n_dofs, 1.0);
-  
-#if deal_II_dimension > 1
-  const unsigned int dofs_per_face = this->dofs_per_face;
-
-  if (dim == 2)
-    {
-      for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2; 
-          f < GeometryInfo<dim>::faces_per_cell; ++f)
-       {
-         typename Triangulation<dim>::face_iterator face = cell->face (f);
-         if (!face->at_boundary ())
-           {
-             unsigned int nn = cell->neighbor_of_neighbor (f);
-             if (nn < GeometryInfo<dim>::faces_per_cell / 2)
-               {
-                 for (unsigned int j = 0; j < dofs_per_face; ++j)
-                   sign_change[f * dofs_per_face + j] = -1.0;
-               }
-           }
-       }
-    }
-  else
-    {
-      // TODO: Think about 3D!.
-    }
-#endif
+  get_face_sign_change (cell, sign_change);
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
@@ -467,31 +512,7 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
   // Compute eventual sign changes depending on the neighborhood
   // between two faces.
   std::vector<double> sign_change (n_dofs, 1.0);
-#if deal_II_dimension > 1
-  const unsigned int dofs_per_face = this->dofs_per_face;
-  
-  if (dim == 2)
-    {
-      for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2; 
-          f < GeometryInfo<dim>::faces_per_cell; ++f)
-       {
-         typename Triangulation<dim>::face_iterator face = cell->face (f);
-         if (!face->at_boundary ())
-           {
-             unsigned int nn = cell->neighbor_of_neighbor (f);
-             if (nn < GeometryInfo<dim>::faces_per_cell / 2)
-               {
-                 for (unsigned int j = 0; j < dofs_per_face; ++j)
-                   sign_change[f * dofs_per_face + j] = -1.0;
-               }
-           }
-       }
-    }
-  else
-    {
-      // TODO: Think about 3D!.
-    }
-#endif
+  get_face_sign_change (cell, sign_change);
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
@@ -527,7 +548,7 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
 //TODO: Note that the Jacobian we compute here is that of the transformation from the reference FACE to the real FACE. However, what we need is the Jacobian of the transformation of the reference to the real CELL, which differs by one power of h. Consequently, the tests fe/rt_{11,13,14,15} presently compute the wrong results
                        double J = 1.0;
                        if (mapping_type == contravariant)
-                         J = data.JxW_values[k] / quadrature.weight(k);
+                         J = data.cell_JxW_values[k] / quadrature.weight(k);
 
                        for (unsigned int d=0; d<dim; ++d)
                          data.shape_values(first+d,k) = sign_change[i] * (shape_values[k][d] / J);
@@ -585,7 +606,8 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
                    {
                      // Recompute determinant
 //TODO: Note that the Jacobian we compute here is that of the transformation from the reference FACE to the real FACE. However, what we need is the Jacobian of the transformation of the reference to the real CELL, which differs by one power of h. Consequently, the tests fe/rt_{11,13,14,15} presently compute the wrong results
-                     double J = data.JxW_values[k] / quadrature.weight(k);
+
+                     double J = data.cell_JxW_values[k] / quadrature.weight(k);
                      data.shape_gradients[first+d][k] = sign_change[i] * 
                        shape_grads2[k][d] / J;
                    }
@@ -639,6 +661,8 @@ FE_PolyTensor<POLY,dim>::fill_fe_subface_values (
 //           && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
 //      ExcNotImplemented());
 //TODO: Size assertions
+
+//TODO: Sign change for the face DoFs!
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
index 788b303a18b23e618234a9aff9e9e5214e38ca17..7f0d221eacfb4fd52e4461e71a2524b7d200582a 100644 (file)
@@ -451,10 +451,12 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
   if (flags & update_values)
     out |= update_values             | update_covariant_transformation
                                      | update_contravariant_transformation 
+                                     | update_cell_JxW_values
                                      | update_JxW_values;
   if (flags & update_gradients)
     out |= update_gradients          | update_covariant_transformation 
                                      | update_contravariant_transformation
+                                     | update_cell_JxW_values
                                      | update_JxW_values;
   //TODO: Set update flags appropriately and figure out, how the second
   // derivatives for the RT elements can be computed correctly.
index cd4f6982d308dd8980c3ada6c395e832042d7e43..0509aad8ec9c5a493879eef34236112f96f1069a 100644 (file)
@@ -832,6 +832,7 @@ compute_fill (const Mapping<dim>                   &mapping,
                                           // the reference element to the global
                                           // coordinate system.
          base_data.JxW_values = data.JxW_values;
+         base_data.cell_JxW_values = data.cell_JxW_values;
 
 
                                            // Make sure that in the
index ad6160d4a92436c705c7f7d023cf476e1cd561db..9f7e38aaf155f659d23b5847813ee86fddf0d4d4 100644 (file)
@@ -312,6 +312,9 @@ FEValuesData<dim>::initialize (const unsigned int        n_quadrature_points,
 
   if (flags & update_normal_vectors)
     this->normal_vectors.resize(n_quadrature_points);
+
+  if (flags & update_cell_JxW_values)
+    this->cell_JxW_values.resize(n_quadrature_points);
 }
 
 
@@ -1011,6 +1014,7 @@ FEValuesBase<dim>::memory_consumption () const
          MemoryConsumption::memory_consumption (this->quadrature_points) +
          MemoryConsumption::memory_consumption (this->normal_vectors) +
          MemoryConsumption::memory_consumption (this->boundary_forms) +
+         MemoryConsumption::memory_consumption (this->cell_JxW_values) +
          sizeof(this->update_flags) +
          MemoryConsumption::memory_consumption (n_quadrature_points) +
          MemoryConsumption::memory_consumption (dofs_per_cell) +
@@ -1520,7 +1524,8 @@ void FEFaceValues<dim>::do_reinit (const unsigned int face_no)
                                          this->quadrature_points,
                                          this->JxW_values,
                                          this->boundary_forms,
-                                         this->normal_vectors);
+                                         this->normal_vectors,
+                                         this->cell_JxW_values);
   
   this->get_fe().fill_fe_face_values(this->get_mapping(),
                                     *this->present_cell, face_no,
@@ -1770,7 +1775,8 @@ void FESubfaceValues<dim>::do_reinit (const unsigned int face_no,
                                             this->quadrature_points,
                                             this->JxW_values,
                                             this->boundary_forms,
-                                            this->normal_vectors);
+                                            this->normal_vectors,
+                                            this->cell_JxW_values);
   
   this->get_fe().fill_fe_subface_values(this->get_mapping(),
                                        *this->present_cell,
index 35e43ed841300533c5161b4248ebd3be65e4a692..58a0d82229aa78d2ec14a81addbbf04b0bc35b1a 100644 (file)
@@ -322,7 +322,8 @@ MappingCartesian<dim>::fill_fe_face_values (const typename Triangulation<dim>::c
                                            std::vector<Point<dim> >     &quadrature_points,
                                            std::vector<double>          &JxW_values,
                                            std::vector<Tensor<1,dim> >  &boundary_forms,
-                                           std::vector<Point<dim> >     &normal_vectors) const
+                                           std::vector<Point<dim> >     &normal_vectors,
+                                           std::vector<double>          &/*cell_JxW_values*/) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -351,6 +352,8 @@ MappingCartesian<dim>::fill_fe_face_values (const typename Triangulation<dim>::c
   if (data.current_update_flags() & update_boundary_forms)
     for (unsigned int i=0; i<boundary_forms.size();++i)
       boundary_forms[i] = J * normal_vectors[i];
+
+  //TODO: Implement cell_JxW_values
 }
 
 
@@ -365,7 +368,8 @@ MappingCartesian<dim>::fill_fe_subface_values (const typename Triangulation<dim>
                                               std::vector<Point<dim> >     &quadrature_points,
                                               std::vector<double>          &JxW_values,
                                               std::vector<Tensor<1,dim> >  &boundary_forms,
-                                              std::vector<Point<dim> >     &normal_vectors) const
+                                              std::vector<Point<dim> >     &normal_vectors,
+                                              std::vector<double>          &/*cell_JxW_values*/) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -394,6 +398,8 @@ MappingCartesian<dim>::fill_fe_subface_values (const typename Triangulation<dim>
   if (data.current_update_flags() & update_boundary_forms)
     for (unsigned int i=0; i<boundary_forms.size();++i)
       boundary_forms[i] = J * normal_vectors[i];
+
+  //TODO: Implement cell_JxW_values
 }
 
 
@@ -408,7 +414,8 @@ MappingCartesian<1>::fill_fe_face_values (const Triangulation<1>::cell_iterator
                                          std::vector<Point<1> >&,
                                          std::vector<double>&,
                                          std::vector<Tensor<1,1> >&,
-                                         std::vector<Point<1> >&) const
+                                         std::vector<Point<1> >&,
+                                         std::vector<double>          &/*cell_JxW_values*/) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -424,7 +431,8 @@ MappingCartesian<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterat
                                             std::vector<Point<1> >&,
                                             std::vector<double>&,
                                             std::vector<Tensor<1,1> >&,
-                                            std::vector<Point<1> >&) const
+                                            std::vector<Point<1> >&,
+                                            std::vector<double>          &/*cell_JxW_values*/) const
 {
   Assert(false, ExcNotImplemented());
 }
index c715dc13765fa6a272cb22e986e6acc3032d7fe4..794c1277202f4f060a97b165fca68fe020312488 100644 (file)
@@ -342,7 +342,8 @@ MappingQ<dim>::fill_fe_face_values (const typename Triangulation<dim>::cell_iter
                                    std::vector<Point<dim> >     &quadrature_points,
                                    std::vector<double>          &JxW_values,
                                    std::vector<Tensor<1,dim> >  &exterior_forms,
-                                   std::vector<Point<dim> >     &normal_vectors) const
+                                   std::vector<Point<dim> >     &normal_vectors,
+                                   std::vector<double>          &cell_JxW_values) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -385,7 +386,10 @@ MappingQ<dim>::fill_fe_face_values (const typename Triangulation<dim>::cell_iter
                            q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
-                           exterior_forms, normal_vectors);
+                           exterior_forms, normal_vectors,
+                          cell_JxW_values);
+
+  // TODO: Verify implementation of cell_JxW_values
 }
 
 
@@ -399,7 +403,8 @@ MappingQ<dim>::fill_fe_subface_values (const typename Triangulation<dim>::cell_i
                                       std::vector<Point<dim> >     &quadrature_points,
                                       std::vector<double>          &JxW_values,
                                       std::vector<Tensor<1,dim> >  &exterior_forms,
-                                      std::vector<Point<dim> >     &normal_vectors) const
+                                      std::vector<Point<dim> >     &normal_vectors,
+                                      std::vector<double>          &cell_JxW_values) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -443,7 +448,10 @@ MappingQ<dim>::fill_fe_subface_values (const typename Triangulation<dim>::cell_i
                            q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
-                           exterior_forms, normal_vectors);
+                           exterior_forms, normal_vectors,
+                          cell_JxW_values);
+
+  // TODO: Verify implementation of cell_JxW_values ...
 }
 
 
index e4acb4ce35702745077105ca9bb3efa50edaa736..b5bafeccbf3f713d451668e74ec8b0dd774059ef 100644 (file)
@@ -249,6 +249,7 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
                                      | update_covariant_transformation
                                      | update_contravariant_transformation
                                      | update_JxW_values
+                                     | update_cell_JxW_values
                                      | update_boundary_forms
                                      | update_normal_vectors));
 
@@ -271,6 +272,15 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
             | update_normal_vectors))
     out |= update_contravariant_transformation;
 
+                                  // The contravariant transformation
+                                  // is a Piola transformation, which
+                                   // requires the determinant of the
+                                   // Jacobi matrix of the transformation.
+                                   // Therefore these values have to
+                                   // updated for each cell.
+  if (out & update_contravariant_transformation)
+    out |= update_JxW_values | update_cell_JxW_values;
+
   return out;
 }
 
@@ -582,7 +592,8 @@ MappingQ1<dim>::compute_fill_face (const typename Triangulation<dim>::cell_itera
                                   std::vector<Point<dim> >    &quadrature_points,
                                   std::vector<double>         &JxW_values,
                                   std::vector<Tensor<1,dim> > &boundary_forms,
-                                  std::vector<Point<dim> >    &normal_vectors) const
+                                  std::vector<Point<dim> >    &normal_vectors,
+                                  std::vector<double>         &cell_JxW_values) const
 {
   compute_fill (cell, n_q_points, data_set, data, quadrature_points);
 
@@ -658,6 +669,23 @@ MappingQ1<dim>::compute_fill_face (const typename Triangulation<dim>::cell_itera
               normal_vectors[i] = boundary_forms[i] / f;
          }
     }
+
+  // If the Piola transformation is to be used, the
+  // new flag, update_cell_JxW_values has to be set.
+  // It triggers the creation of values for the determinant
+  // of the mapping function on 
+  if (update_flags & update_cell_JxW_values)
+    {      
+      Assert (cell_JxW_values.size() == n_q_points,
+             ExcDimensionMismatch(cell_JxW_values.size(), n_q_points));
+      Assert (data.contravariant.size() == n_q_points,
+             ExcDimensionMismatch(data.contravariant.size(), n_q_points));
+      Assert (weights.size() == n_q_points,
+             ExcDimensionMismatch(weights.size(), n_q_points));
+      for (unsigned int point=0; point<n_q_points; ++point)
+       cell_JxW_values[point]
+         = determinant(data.contravariant[point])*weights[point];
+    }
 }
 
 
@@ -670,7 +698,8 @@ MappingQ1<dim>::fill_fe_face_values (const typename Triangulation<dim>::cell_ite
                                     std::vector<Point<dim> >     &quadrature_points,
                                     std::vector<double>          &JxW_values,
                                     std::vector<Tensor<1,dim> >  &boundary_forms,
-                                    std::vector<Point<dim> >     &normal_vectors) const
+                                    std::vector<Point<dim> >     &normal_vectors,
+                                    std::vector<double>          &cell_JxW_values) const
 {
   InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
   Assert(data_ptr!=0, ExcInternalError());
@@ -688,7 +717,8 @@ MappingQ1<dim>::fill_fe_face_values (const typename Triangulation<dim>::cell_ite
                     quadrature_points,
                     JxW_values,
                     boundary_forms,
-                    normal_vectors);
+                    normal_vectors,
+                    cell_JxW_values);
 }
 
 
@@ -702,7 +732,8 @@ MappingQ1<dim>::fill_fe_subface_values (const typename Triangulation<dim>::cell_
                                        std::vector<Point<dim> >     &quadrature_points,
                                        std::vector<double>          &JxW_values,
                                        std::vector<Tensor<1,dim> >  &boundary_forms,
-                                       std::vector<Point<dim> >     &normal_vectors) const
+                                       std::vector<Point<dim> >     &normal_vectors,
+                                       std::vector<double>          &cell_JxW_values) const
 {
   InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
   Assert(data_ptr!=0, ExcInternalError());
@@ -720,7 +751,8 @@ MappingQ1<dim>::fill_fe_subface_values (const typename Triangulation<dim>::cell_
                     quadrature_points,
                     JxW_values,
                     boundary_forms,
-                    normal_vectors);
+                    normal_vectors,
+                    cell_JxW_values);
 }
 
 
@@ -738,7 +770,8 @@ MappingQ1<1>::compute_fill_face (const Triangulation<1>::cell_iterator &,
                                 std::vector<Point<1> > &,
                                 std::vector<double> &,
                                 std::vector<Tensor<1,1> > &,
-                                std::vector<Point<1> > &) const
+                                std::vector<Point<1> > &,
+                                std::vector<double>          &/*cell_JxW_values*/) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -753,7 +786,8 @@ MappingQ1<1>::fill_fe_face_values (const Triangulation<1>::cell_iterator &,
                                   std::vector<Point<1> >&,
                                   std::vector<double>&,
                                   std::vector<Tensor<1,1> >&,
-                                  std::vector<Point<1> >&) const
+                                  std::vector<Point<1> >&,
+                                  std::vector<double>          &/*cell_JxW_values*/) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -769,7 +803,8 @@ MappingQ1<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterator &,
                                      std::vector<Point<1> >&,
                                      std::vector<double>&,
                                      std::vector<Tensor<1,1> >&,
-                                     std::vector<Point<1> >&) const
+                                     std::vector<Point<1> >&,
+                                     std::vector<double>          &/*cell_JxW_values*/) const
 {
   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.