]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move computation of offsets in QProjections into the class that does the projections...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Oct 2003 19:20:41 +0000 (19:20 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Oct 2003 19:20:41 +0000 (19:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@8108 0785d39b-7218-0410-832d-ea1e28bc413d

20 files changed:
deal.II/base/include/base/quadrature.h
deal.II/base/source/quadrature.cc
deal.II/deal.II/include/fe/mapping_q1.h
deal.II/deal.II/include/grid/geometry_info.h
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/deal.II/source/fe/fe_dgp.cc
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_nedelec.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/fe_q_hierarchical.cc
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_tools.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
deal.II/deal.II/source/numerics/error_estimator.cc
deal.II/doc/news/c-4-0.html
deal.II/examples/step-12/step-12.cc
deal.II/examples/step-14/step-14.cc

index 1a1d59946b7b0767d7ad21deae2315b718154fe7..36a4bf13163679bce0b9ff0fdc29df2cc0a65f41 100644 (file)
@@ -267,7 +267,7 @@ class QIterated : public Quadrature<dim>
 /**
  *  This class is a helper class to facilitate the usage of quadrature formulae
  *  on faces or subfaces of cells. It computes the locations of quadrature
- *  points on the unit cell from a quadrature object for a mannifold of
+ *  points on the unit cell from a quadrature object for a manifold of
  *  one dimension less than that of the cell and the number of the face.
  *  For example, giving the Simpson rule in one dimension and using the
  *  @p{project_to_face} function with face number 1, the returned points will
@@ -280,19 +280,34 @@ class QIterated : public Quadrature<dim>
  *  unit cell), the same applies as above. Note the order in which the
  *  children of a face are numbered, which in two dimensions coincides
  *  with the orientation of the face.
+ *
+ *  The second set of functions generates a quadrature formula by
+ *  projecting a given quadrature rule on @em{all} faces and
+ *  subfaces. This is used in the @ref{FEFaceValues} and
+ *  @ref{FESubfaceValues} classes. Since we now have the quadrature
+ *  points of all faces and subfaces in one array, we need to have a
+ *  way to find the starting index of the points and weights
+ *  corresponding to one face or subface within this array. This is
+ *  done through the @ref{DataSetDescriptor} member class.
  *  
- *  The different functions are grouped into a common class to avoid putting
- *  them into global namespace (and to make documentation easier, since
- *  presently the documentation tool can only handle classes, not global
- *  functions). However, since they have no local data, all functions are
- *  declared @p{static} and can be called without creating an object of this
- *  class.
+ *  The different functions are grouped into a common class to avoid
+ *  putting them into global namespace. However, since they have no
+ *  local data, all functions are declared @p{static} and can be
+ *  called without creating an object of this class.
+ *
+ *  For the 3d case, you should note that the orientation of faces is
+ *  even more intricate than for two dimensions. Quadrature formulae
+ *  are projected upon the faces in their standard orientation, not to
+ *  the inside or outside of the hexahedron. Refer to the
+ *  documentation of the @p{Triangulation} class for a description of
+ *  the orientation of the different faces. To make things more
+ *  complicated, in 3d we allow faces in two orientations (which can
+ *  be identified using @p{cell->face_orientation(face)}), so we have
+ *  to project quadrature formula onto faces and subfaces in two
+ *  orientations. The @ref{DataSetDescriptor} member class is used to
+ *  identify where each dataset starts.
  *
- *  For the 3d case, you should note that the orientation of faces is even
- *  more intricate than for two dimensions. Quadrature formulae are projected
- *  upon the faces in their standard orientation, not to the inside or outside
- *  of the hexahedron! Refer to the documentation of the @p{Triangulation} class
- *  for a description of the orientation of the different faces.
+ *  @author Wolfgang Bangerth, 1998, 1999, 2003
  */
 template <int dim>
 class QProjector
@@ -399,6 +414,147 @@ class QProjector
     project_to_child (const Quadrature<dim>  &quadrature,
                      const unsigned int      child_no);
 
+                                     /**
+                                      * Since the
+                                      * @p{project_to_all_faces} and
+                                      * @p{project_to_all_subfaces}
+                                      * functions chain together the
+                                      * quadrature points and weights
+                                      * of all projections of a face
+                                      * quadrature formula to the
+                                      * faces or subfaces of a cell,
+                                      * we need a way to identify
+                                      * where the starting index of
+                                      * the points and weights for a
+                                      * particular face or subface
+                                      * is. This class provides this:
+                                      * there are static member
+                                      * functions that generate
+                                      * objects of this type, given
+                                      * face or subface indices, and
+                                      * you can then use the generated
+                                      * object in place of an integer
+                                      * that denotes the offset of a
+                                      * given dataset.
+                                      *
+                                      * @author Wolfgang Bangerth, 2003
+                                      */
+    class DataSetDescriptor
+    {
+      public:
+                                         /**
+                                          * Default constructor. This
+                                          * doesn't do much except
+                                          * generating an invalid
+                                          * index, since you didn't
+                                          * give a valid descriptor of
+                                          * the cell, face, or subface
+                                          * you wanted.
+                                          */
+        DataSetDescriptor ();
+
+                                         /**
+                                          * Static function to
+                                          * generate the offset of a
+                                          * cell. Since we only have
+                                          * one cell per quadrature
+                                          * object, this offset is of
+                                          * course zero, but we carry
+                                          * this function around for
+                                          * consistency with the other
+                                          * static functions.
+                                          */
+        static DataSetDescriptor cell ();
+      
+                                         /**
+                                          * Static function to
+                                          * generate an offset object
+                                          * for a given face of a cell
+                                          * with the given face
+                                          * orientation. This function
+                                          * of course is only allowed
+                                          * if @p{dim>=2}, and the
+                                          * face orientation is
+                                          * ignored if the space
+                                          * dimension equals 2.
+                                          *
+                                          * The last argument denotes
+                                          * the number of quadrature
+                                          * points the
+                                          * lower-dimensional face
+                                          * quadrature formula (the
+                                          * one that has been
+                                          * projected onto the faces)
+                                          * has.
+                                          */
+        static
+        DataSetDescriptor
+        face (const unsigned int face_no,
+              const bool         face_orientation,
+              const unsigned int n_quadrature_points);
+
+                                         /**
+                                          * Static function to
+                                          * generate an offset object
+                                          * for a given subface of a
+                                          * cell with the given face
+                                          * orientation. This function
+                                          * of course is only allowed
+                                          * if @p{dim>=2}, and the
+                                          * face orientation is
+                                          * ignored if the space
+                                          * dimension equals 2.
+                                          *
+                                          * The last argument denotes
+                                          * the number of quadrature
+                                          * points the
+                                          * lower-dimensional face
+                                          * quadrature formula (the
+                                          * one that has been
+                                          * projected onto the faces)
+                                          * has.
+                                          */
+        static
+        DataSetDescriptor
+        sub_face (const unsigned int face_no,
+                  const unsigned int subface_no,
+                  const bool         face_orientation,
+                  const unsigned int n_quadrature_points);
+
+                                         /**
+                                          * Conversion operator to an
+                                          * integer denoting the
+                                          * offset of the first
+                                          * element of this dataset in
+                                          * the set of quadrature
+                                          * formulas all projected
+                                          * onto faces and
+                                          * subfaces. This conversion
+                                          * operator allows us to use
+                                          * offset descriptor objects
+                                          * in place of integer
+                                          * offsets.
+                                          */
+        operator unsigned int () const;
+      
+      private:
+                                         /**
+                                          * Store the integer offset
+                                          * for a given cell, face, or
+                                          * subface.
+                                          */
+        const unsigned int dataset_offset;
+
+                                         /**
+                                          * This is the real
+                                          * constructor, but it is
+                                          * private and thus only
+                                          * available to the static
+                                          * member functions above.
+                                          */
+        DataSetDescriptor (const unsigned int dataset_offset);
+    };
+
   private:
                                      /**
                                       * Given a quadrature object in
index 182c84d4ac5dbf43b9d2ddd6c258dda99223e6b3..b591d3c5c53cd26d908facf30689385f7e1c38be 100644 (file)
@@ -878,6 +878,120 @@ QProjector<dim>::project_to_child (const Quadrature<dim>    &quadrature,
 
 
 
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::cell ()
+{
+  return 0;
+}
+
+
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::
+face (const unsigned int face_no,
+      const bool         face_orientation,
+      const unsigned int n_quadrature_points)
+{
+  Assert (dim != 1, ExcInternalError());
+  Assert (face_no < GeometryInfo<dim>::faces_per_cell,
+          ExcInternalError());
+  Assert (n_quadrature_points > 0, ExcInternalError());
+    
+  switch (dim)
+    {
+      case 1:
+      case 2:
+            return face_no * n_quadrature_points;
+
+                                             // in 3d, we have to
+                                             // account for faces
+                                             // that have reverse
+                                             // orientation. thus,
+                                             // we have to store
+                                             // _two_ data sets per
+                                             // face or subface
+      case 3:
+            return ((face_no +
+                     (face_orientation == true ?
+                      0 : GeometryInfo<dim>::faces_per_cell))
+                    * n_quadrature_points);
+
+      default:
+            Assert (false, ExcInternalError());
+    }
+  return static_cast<unsigned int>(-1);
+}
+
+
+
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::
+sub_face (const unsigned int face_no,
+          const unsigned int subface_no,
+          const bool         face_orientation,
+          const unsigned int n_quadrature_points)
+{
+  Assert (dim != 1, ExcInternalError());
+  Assert (face_no < GeometryInfo<dim>::faces_per_cell,
+          ExcInternalError());
+                                   // the trick with +1 prevents
+                                   // that we get a warning in 1d
+  Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
+          ExcInternalError());
+  Assert (n_quadrature_points > 0, ExcInternalError());
+
+  switch (dim)
+    {
+      case 1:
+      case 2:
+            return ((face_no * GeometryInfo<dim>::subfaces_per_face +
+                     subface_no)
+                    * n_quadrature_points);
+
+                                             // for 3d, same as above:
+      case 3:
+            return (((face_no * GeometryInfo<dim>::subfaces_per_face +
+                      subface_no)
+                     + (face_orientation == true ?
+                        0 :
+                        GeometryInfo<dim>::faces_per_cell *
+                        GeometryInfo<dim>::subfaces_per_face)
+                     )
+                    * n_quadrature_points);
+      default:
+            Assert (false, ExcInternalError());
+    }
+  return static_cast<unsigned int>(-1);              
+}
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::operator unsigned int () const
+{
+  return dataset_offset;
+}
+
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::
+DataSetDescriptor (const unsigned int dataset_offset)
+                :
+                dataset_offset (dataset_offset)
+{}
+
+
+template <int dim>
+QProjector<dim>::DataSetDescriptor::
+DataSetDescriptor ()
+                :
+                dataset_offset (static_cast<unsigned int>(-1))
+{}
+
+
+
 // ------------------------------------------------------------ //
 
 
index fc7080e53a850529b3a21488afff5f65b11a90c4..f3dce621a59d2a9b4aee8f54fa24be1c659e3885 100644 (file)
 
 #include <base/config.h>
 #include <base/table.h>
-#include <cmath>
+#include <base/quadrature.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 #include <fe/mapping.h>
 
+#include <cmath>
 
-/*!@addtogroup fe */
-/*@{*/
 
-namespace internal
-{
-  template <int dim>
-  class DataSetDescriptor
-  {
-    public:
-      DataSetDescriptor ();
-      
-      static DataSetDescriptor cell ();
-      
-      static
-      DataSetDescriptor
-      face (const typename DoFHandler<dim>::cell_iterator &cell,
-            const unsigned int face_no,
-            const unsigned int n_quadrature_points);
-
-      static
-      DataSetDescriptor
-      sub_face (const typename DoFHandler<dim>::cell_iterator &cell,
-                const unsigned int face_no,
-                const unsigned int subface_no,
-                const unsigned int n_quadrature_points);
-
-      unsigned int offset () const;
-      
-    private:
-      const unsigned int dataset_offset;
-
-      DataSetDescriptor (const unsigned int dataset_offset);
-  };
-}
 
 
 
+/*!@addtogroup fe */
+/*@{*/
+
+
 /**
  * Mapping of general quadrilateral/hexahedra by d-linear shape
  * functions.
@@ -385,12 +357,17 @@ class MappingQ1 : public Mapping<dim>
     DeclException0 (ExcAccessToUninitializedField);
 
   protected:
+
                                      /**
-                                      * Typedef the right data set
-                                      * descriptor to a local type to
-                                      * make use somewhat simpler.
+                                      * Declare a convenience typedef
+                                      * for the class that describes
+                                      * offsets into quadrature
+                                      * formulas projected onto faces
+                                      * and subfaces.
                                       */
-    typedef internal::DataSetDescriptor<dim> DataSetDescriptor;    
+    typedef
+    typename QProjector<dim>::DataSetDescriptor
+    DataSetDescriptor;
     
                                     /**
                                      * Implementation of the interface in
index ba01b9872509cbd8812ee45d98e4bce0f051e63e..fe10fd3ca87f313c126e4085863b444ec834c25e 100644 (file)
@@ -159,9 +159,34 @@ struct GeometryInfo
                                      * 3 and 2, etc.
                                      *
                                      * For three spatial dimensions,
-                                     * the exact order of the children is
-                                     * laid down in the documentation of
-                                     * the @ref{Triangulation} class.
+                                     * the exact order of the
+                                     * children is laid down in the
+                                     * documentation of the
+                                     * @ref{Triangulation} class.
+                                     * However, it must be noted that
+                                     * this class and function only
+                                     * deals with faces in standard
+                                     * orientation. In 3d, faces can
+                                     * exist in two orientations,
+                                     * though, and if a face is in
+                                     * the wrong orientation, then
+                                     * this function may not give you
+                                     * what you want. You can inquire
+                                     * about the face orientation
+                                     * using the
+                                     * @p{cell->face_orientation}
+                                     * function, and the function to
+                                     * ask for a neighbor's cell
+                                     * behind a given face and
+                                     * subface is
+                                     * @p{cell->neighbor_child_on_subface}.
+                                     * The latter function, in
+                                     * contrast to the present one,
+                                     * also takes into account the
+                                     * actual orientation of the
+                                     * faces of a cell and will
+                                     * return the correct result in
+                                     * all cases.
                                      */
     static unsigned int child_cell_on_face (const unsigned int face,
                                            const unsigned int subface);
@@ -510,9 +535,35 @@ struct GeometryInfo<1>
                                      * 3 and 2, etc.
                                      *
                                      * For three spatial dimensions,
-                                     * the exact order of the children is
-                                     * laid down in the documentation of
-                                     * the @ref{Triangulation} class.
+                                     * the exact order of the
+                                     * children is laid down in the
+                                     * documentation of the
+                                     * @ref{Triangulation}
+                                     * class. However, it must be
+                                     * noted that this class and
+                                     * function only deals with faces
+                                     * in standard orientation. In
+                                     * 3d, faces can exist in two
+                                     * orientations, though, and if a
+                                     * face is in the wrong
+                                     * orientation, then this
+                                     * function may not give you what
+                                     * you want. You can inquire
+                                     * about the face orientation
+                                     * using the
+                                     * @p{cell->face_orientation}
+                                     * function, and the function to
+                                     * ask for a neighbor's cell
+                                     * behind a given face and
+                                     * subface is
+                                     * @p{cell->neighbor_child_on_subface}.
+                                     * The latter function, in
+                                     * contrast to the present one,
+                                     * also takes into account the
+                                     * actual orientation of the
+                                     * faces of a cell and will
+                                     * return the correct result in
+                                     * all cases.
                                      */
     static unsigned int child_cell_on_face (const unsigned int face,
                                            const unsigned int subface);
@@ -769,9 +820,35 @@ struct GeometryInfo<2>
                                      * 3 and 2, etc.
                                      *
                                      * For three spatial dimensions,
-                                     * the exact order of the children is
-                                     * laid down in the documentation of
-                                     * the @ref{Triangulation} class.
+                                     * the exact order of the
+                                     * children is laid down in the
+                                     * documentation of the
+                                     * @ref{Triangulation}
+                                     * class. However, it must be
+                                     * noted that this class and
+                                     * function only deals with faces
+                                     * in standard orientation. In
+                                     * 3d, faces can exist in two
+                                     * orientations, though, and if a
+                                     * face is in the wrong
+                                     * orientation, then this
+                                     * function may not give you what
+                                     * you want. You can inquire
+                                     * about the face orientation
+                                     * using the
+                                     * @p{cell->face_orientation}
+                                     * function, and the function to
+                                     * ask for a neighbor's cell
+                                     * behind a given face and
+                                     * subface is
+                                     * @p{cell->neighbor_child_on_subface}.
+                                     * The latter function, in
+                                     * contrast to the present one,
+                                     * also takes into account the
+                                     * actual orientation of the
+                                     * faces of a cell and will
+                                     * return the correct result in
+                                     * all cases.
                                      */
     static unsigned int child_cell_on_face (const unsigned int face,
                                            const unsigned int subface);
@@ -1030,9 +1107,34 @@ struct GeometryInfo<3>
                                      * 3 and 2, etc.
                                      *
                                      * For three spatial dimensions,
-                                     * the exact order of the children is
-                                     * laid down in the documentation of
-                                     * the @ref{Triangulation} class.
+                                     * the exact order of the
+                                     * children is laid down in the
+                                     * documentation of the
+                                     * @ref{Triangulation}
+                                     * class. However, it must be
+                                     * noted that this class and
+                                     * function only deals with faces
+                                     * in standard orientation. In
+                                     * 3d, faces can exist in two
+                                     * orientations, though, and if a
+                                     * face is in the wrong
+                                     * orientation, then this
+                                     * function may not give you what
+                                     * you want. You can inquire
+                                     * about the face orientation
+                                     * using the
+                                     * @p{cell->face_orientation}
+                                     * function, and the function to
+                                     * ask for a neighbor's cell
+                                     * behind a given face and
+                                     * subface is
+                                     * @p{cell->neighbor_child_on_subface}.
+                                     * The latter function, in contrast
+                                     * to the present one, also takes
+                                     * into account the actual
+                                     * orientation of the faces of a
+                                     * cell and will return the
+                                     * correct result in all cases.
                                      */
     static unsigned int child_cell_on_face (const unsigned int face,
                                            const unsigned int subface);
index de6186001ecae47df32c142d18f99c0063c23b51..aaa7d1aaea224e427c8d741d33a00e8a15428ff0 100644 (file)
@@ -373,7 +373,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
              if (neighbor->level() < cell->level())
                continue;
 
-             unsigned int neighbor_face = cell->neighbor_of_neighbor(face);
+             const unsigned int neighbor_face
+                = cell->neighbor_of_neighbor(face);
 
              if (neighbor->has_children())
                {
@@ -381,35 +382,35 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                       sub_nr != GeometryInfo<dim>::subfaces_per_face;
                       ++sub_nr)
                    {
-                     typename DoFHandler<dim>::cell_iterator sub_neighbor
-                       = neighbor->
-                       child(GeometryInfo<dim>::child_cell_on_face(neighbor_face, sub_nr));
+                     const typename DoFHandler<dim>::cell_iterator
+                        sub_neighbor
+                       = cell->neighbor_child_on_subface (face, sub_nr);
 
                      sub_neighbor->get_dof_indices (dofs_on_other_cell);
                      for (unsigned int i=0; i<total_dofs; ++i)
-                       {
-                         for (unsigned int j=0; j<total_dofs; ++j)
-                           {
-                             sparsity.add (dofs_on_this_cell[i],
-                                           dofs_on_other_cell[j]);
-                             sparsity.add (dofs_on_other_cell[i],
-                                           dofs_on_this_cell[j]);
-                           }
-                       }
+                        for (unsigned int j=0; j<total_dofs; ++j)
+                          {
+                            sparsity.add (dofs_on_this_cell[i],
+                                          dofs_on_other_cell[j]);
+                            sparsity.add (dofs_on_other_cell[i],
+                                          dofs_on_this_cell[j]);
+                          }
+                       
                      sub_neighbor->face(neighbor_face)->set_user_flag ();
                    }
-               } else {
+               }
+              else
+                {
                  neighbor->get_dof_indices (dofs_on_other_cell);
                  for (unsigned int i=0; i<total_dofs; ++i)
-                   {
-                     for (unsigned int j=0; j<total_dofs; ++j)
-                       {
-                         sparsity.add (dofs_on_this_cell[i],
-                                       dofs_on_other_cell[j]);
-                         sparsity.add (dofs_on_other_cell[i],
-                                       dofs_on_this_cell[j]);
-                       }
-                   }
+                    for (unsigned int j=0; j<total_dofs; ++j)
+                      {
+                        sparsity.add (dofs_on_this_cell[i],
+                                      dofs_on_other_cell[j]);
+                        sparsity.add (dofs_on_other_cell[i],
+                                      dofs_on_this_cell[j]);
+                      }
+               
                  neighbor->face(neighbor_face)->set_user_flag (); 
                }
            } 
@@ -575,7 +576,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
           face < GeometryInfo<dim>::faces_per_cell;
           ++face)
        {
-         typename DoFHandler<dim>::face_iterator cell_face = cell->face(face);
+         const typename DoFHandler<dim>::face_iterator
+            cell_face = cell->face(face);
          if (cell_face->user_flag_set ())
            continue;
 
@@ -589,16 +591,12 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                      const bool j_non_zero_i = support_on_face (j, face);
                      
                      if (flux_dof_mask(i,j) == 'f')
-                       {
-                         sparsity.add (dofs_on_this_cell[i],
-                                       dofs_on_this_cell[j]);
-                       }
+                        sparsity.add (dofs_on_this_cell[i],
+                                      dofs_on_this_cell[j]);
                      if (flux_dof_mask(i,j) == 'a')
-                       {
-                         if (i_non_zero_i && j_non_zero_i)
-                           sparsity.add (dofs_on_this_cell[i],
-                                         dofs_on_this_cell[j]);
-                       }
+                        if (i_non_zero_i && j_non_zero_i)
+                          sparsity.add (dofs_on_this_cell[i],
+                                        dofs_on_this_cell[j]);
                    }
                }
            }
@@ -611,7 +609,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
              if (neighbor->level() < cell->level())
                continue;
 
-             unsigned int neighbor_face = cell->neighbor_of_neighbor(face);
+             const unsigned int
+                neighbor_face = cell->neighbor_of_neighbor(face);
 
              if (neighbor->has_children())
                {
@@ -619,9 +618,9 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                       sub_nr != GeometryInfo<dim>::subfaces_per_face;
                       ++sub_nr)
                    {
-                     typename DoFHandler<dim>::cell_iterator sub_neighbor
-                       = neighbor->
-                       child(GeometryInfo<dim>::child_cell_on_face(neighbor_face, sub_nr));
+                     const typename DoFHandler<dim>::cell_iterator
+                        sub_neighbor
+                       = cell->neighbor_child_on_subface (face, sub_nr);
 
                      sub_neighbor->get_dof_indices (dofs_on_other_cell);
                      for (unsigned int i=0; i<total_dofs; ++i)
@@ -693,7 +692,9 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                        }
                      sub_neighbor->face(neighbor_face)->set_user_flag ();
                    }
-               } else {
+               }
+              else
+                {
                  neighbor->get_dof_indices (dofs_on_other_cell);
                  for (unsigned int i=0; i<total_dofs; ++i)
                    {
index 0c772f15445419e5ab2616f2ccae08d82aaeb860..5106e32986e69ff30eb08680255a8d0fd6318faf 100644 (file)
@@ -18,7 +18,6 @@
 #include <fe/fe.h>
 #include <fe/mapping.h>
 #include <fe/fe_dgp.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_values.h>
 
 #ifdef HAVE_STD_STRINGSTREAM
@@ -352,7 +351,7 @@ FE_DGP<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
   
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -377,10 +376,10 @@ FE_DGP<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
@@ -428,10 +427,10 @@ FE_DGP<dim>::fill_fe_subface_values (const Mapping<dim>                   &mappi
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // sub-faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face), 
+                 quadrature.n_quadrature_points));
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
index 3accfb98c43942ce04aafa0ec4c4482f51a877c7..79ace133ba9c0ff137df4f2a7a5a0f986d0dc45a 100644 (file)
@@ -18,7 +18,6 @@
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 #include <fe/fe.h>
-#include <fe/mapping_q1.h>
 #include <fe/mapping.h>
 #include <fe/fe_dgq.h>
 #include <fe/fe_values.h>
@@ -766,7 +765,7 @@ FE_DGQ<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
   
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -791,10 +790,10 @@ FE_DGQ<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
@@ -842,10 +841,10 @@ FE_DGQ<dim>::fill_fe_subface_values (const Mapping<dim>                   &mappi
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // sub-faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face),
+                 quadrature.n_quadrature_points));
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
index 03a21ead75e48b440388ed6acd42baaa928b4758..de12d5955cc6d97e467c5651d1c1ddf28d4186dc 100644 (file)
@@ -18,7 +18,6 @@
 #include <dofs/dof_accessor.h>
 #include <fe/fe.h>
 #include <fe/mapping.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_nedelec.h>
 #include <fe/fe_values.h>
 
@@ -1098,7 +1097,7 @@ FE_Nedelec<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
 
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -1123,10 +1122,10 @@ FE_Nedelec<dim>::fill_fe_face_values (const Mapping<dim>                   &mapp
                                    // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
 
                                   // get the flags indicating the
                                   // fields that have to be filled
@@ -1273,10 +1272,10 @@ FE_Nedelec<dim>::fill_fe_subface_values (const Mapping<dim>                   &m
                                    // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face),
+                 quadrature.n_quadrature_points));
 
                                   // get the flags indicating the
                                   // fields that have to be filled
index 572b02e9a44280e6266360b910c1fe8b10311ee4..f0a3aaa2ea28c565b863f8ae8dd2ed735ddac19d 100644 (file)
@@ -19,7 +19,6 @@
 #include <dofs/dof_accessor.h>
 #include <fe/fe.h>
 #include <fe/mapping.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
 
@@ -1689,7 +1688,7 @@ FE_Q<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
 
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -1714,10 +1713,10 @@ FE_Q<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
@@ -1765,10 +1764,10 @@ FE_Q<dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // sub-faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face),
+                 quadrature.n_quadrature_points));
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
index 39cb5bd1e905137013164a614d5f57c862611792..5a4ff175998276bbfaddd155ad874b5523d27a20 100644 (file)
@@ -19,7 +19,6 @@
 #include <dofs/dof_accessor.h>
 #include <fe/fe.h>
 #include <fe/mapping.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_q_hierarchical.h>
 #include <fe/fe_values.h>
 
@@ -1061,7 +1060,7 @@ FE_Q_Hierarchical<dim>::fill_fe_values (
 
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -1087,10 +1086,10 @@ FE_Q_Hierarchical<dim>::fill_fe_face_values (
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
@@ -1139,10 +1138,10 @@ FE_Q_Hierarchical<dim>::fill_fe_subface_values (
                                   // offset determines which data set
                                   // to take (all data sets for all
                                   // sub-faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face),
+                 quadrature.n_quadrature_points));
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
index 28d837df1407ba26703da46d0a45cdebd13bda07..b63f5b854b2d40922a2097157035fdae2a115dd7 100644 (file)
@@ -18,7 +18,6 @@
 #include <dofs/dof_accessor.h>
 #include <fe/fe.h>
 #include <fe/mapping.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_raviart_thomas.h>
 #include <fe/fe_values.h>
 
@@ -1611,7 +1610,7 @@ FE_RaviartThomas<dim>::fill_fe_values (const Mapping<dim>                   &map
 
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell,
-                       internal::DataSetDescriptor<dim>::cell().offset(),
+                       QProjector<dim>::DataSetDescriptor::cell(),
                        mapping_data, fe_data, data);
 }
 
@@ -1636,10 +1635,10 @@ FE_RaviartThomas<dim>::fill_fe_face_values (const Mapping<dim>
                                    // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       face (cell, face,
-             quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       face (face, cell->face_orientation(face),
+             quadrature.n_quadrature_points));
 
                                   // get the flags indicating the
                                   // fields that have to be filled
@@ -1778,10 +1777,10 @@ FE_RaviartThomas<dim>::fill_fe_subface_values (const Mapping<dim>
                                    // offset determines which data set
                                   // to take (all data sets for all
                                   // faces are stored contiguously)
-  const unsigned int offset
-    = (internal::DataSetDescriptor<dim>::
-       sub_face (cell, face, subface,
-                 quadrature.n_quadrature_points)).offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (QProjector<dim>::DataSetDescriptor::
+       sub_face (face, subface, cell->face_orientation(face),
+                 quadrature.n_quadrature_points));
 
                                   // get the flags indicating the
                                   // fields that have to be filled
index 461ffb1fde698d0db9917f3ea07de19a038f47f6..7b9efaeb96548d81ea56a9a3cce4628eb4d3a77e 100644 (file)
@@ -29,7 +29,6 @@
 #include <fe/fe_raviart_thomas.h>
 #include <fe/fe_system.h>
 #include <fe/fe_values.h>
-#include <fe/mapping_q1.h>
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
index 0ea0cdda921ec4d8f21db4f59538f6ef8a2973f6..6ae722cd5406470c0c5d9b515c26b69f69435527 100644 (file)
@@ -865,6 +865,10 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
          ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
   Assert (subface_no < GeometryInfo<dim>::subfaces_per_face,
          ExcIndexRange (subface_no, 0, GeometryInfo<dim>::subfaces_per_face));
+  Assert (cell->has_children() == false,
+          ExcMessage ("You can't use subface data for cells that are "
+                      "already refined. Iterate over their children "
+                      "instead in these cases."));
   
   this->present_cell  = cell;
   this->my_orientation = this->orientation_table[face_no];
index c7ebc882903f2cc0e9c4d328cc97cf65daf889f9..5bfcb69184d59db26ab246deebe7e070d45fce12 100644 (file)
@@ -19,7 +19,6 @@
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 #include <fe/mapping_cartesian.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_values.h>
 
 #include <cmath>
@@ -135,13 +134,29 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
                                     std::vector<Point<dim> > &quadrature_points,
                                     std::vector<Point<dim> > &normal_vectors) const
 {
-  UpdateFlags update_flags(data.current_update_flags());
+  const UpdateFlags update_flags(data.current_update_flags());
 
   const unsigned int npts = quadrature_points.size ();
 
-  typedef typename internal::DataSetDescriptor<dim> DataSetDescriptor;
-  unsigned int offset = DataSetDescriptor::cell().offset();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = (face_no == invalid_face_number
+       ?
+       QProjector<dim>::DataSetDescriptor::cell()
+       :
+       (sub_no == invalid_face_number
+        ?
+                                        // called from FEFaceValues
+        QProjector<dim>::DataSetDescriptor::face (face_no,
+                                                  cell->face_orientation(face_no),
+                                                  quadrature_points.size())
+        :
+                                        // called from FESubfaceValues
+        QProjector<dim>::DataSetDescriptor::sub_face (face_no, sub_no,
+                                                      cell->face_orientation(face_no),
+                                                      quadrature_points.size())
+        ));
   
+                                   // some more sanity checks
   if (face_no != invalid_face_number)
     {
                                       // Add 1 on both sides of
@@ -151,23 +166,11 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
       Assert (face_no+1 < GeometryInfo<dim>::faces_per_cell+1,
              ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
 
-      if (sub_no == invalid_face_number)
-                                        // called from FEFaceValues
-       offset = DataSetDescriptor::face (cell, face_no,
-                                          quadrature_points.size()).offset();
-      else
-       {
-                                          // called from
-                                           // FESubfaceValue (do the
-                                           // +1 trick to avoid a
-                                           // warning about comparison
-                                           // always being false in
-                                           // 1d)
-         Assert (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
-                 ExcIndexRange (sub_no, 0, GeometryInfo<dim>::subfaces_per_face));
-         offset = DataSetDescriptor::sub_face (cell, face_no, sub_no,
-                                                quadrature_points.size()).offset();
-       }
+      Assert ((sub_no == invalid_face_number)
+              ||
+              (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1),
+              ExcIndexRange (sub_no, 0,
+                             GeometryInfo<dim>::subfaces_per_face));
     }
   else
                                     // invalid face number, so
@@ -365,19 +368,16 @@ MappingCartesian<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell
       J *= data.length[d];
   
   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);
-    }
+    for (unsigned int i=0; i<JxW_values.size();++i)
+      JxW_values[i] = J * q.weight(i);
 
   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];
-    }
+    for (unsigned int i=0; i<boundary_forms.size();++i)
+      boundary_forms[i] = J * normal_vectors[i];
 }
 
 
+
 template <int dim>
 void
 MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
@@ -411,16 +411,12 @@ MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::c
       J *= data.length[d];
   
   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) / GeometryInfo<dim>::subfaces_per_face;
-    }
+    for (unsigned int i=0; i<JxW_values.size();++i)
+      JxW_values[i] = J * q.weight(i) / GeometryInfo<dim>::subfaces_per_face;
 
   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];
-    }
+    for (unsigned int i=0; i<boundary_forms.size();++i)
+      boundary_forms[i] = J * normal_vectors[i];
 }
 
 
@@ -476,8 +472,8 @@ MappingCartesian<dim>::transform_covariant (Tensor<1,dim>       *begin,
     {
       for (unsigned int d=0;d<dim;++d)
        (*begin)[d] = (*src)[d]/data.length[d];
-      begin++;
-      src++;
+      ++begin;
+      ++src;
     }
 }
 
@@ -502,8 +498,8 @@ MappingCartesian<dim>::transform_covariant (Tensor<2,dim>       *begin,
       for (unsigned int d=0; d<dim; ++d)
         for (unsigned int p=0; p<dim; ++p)
        (*begin)[d][p] = (*src)[d][p] / data.length[p];
-      begin++;
-      src++;
+      ++begin;
+      ++src;
     }
 }
 
@@ -561,8 +557,8 @@ MappingCartesian<dim>::transform_contravariant (Tensor<2,dim>       *begin,
       for (unsigned int d=0; d<dim; ++d)
         for (unsigned int p=0; p<dim; ++p)
        (*begin)[d][p] = data.length[d] * (*src)[d][p];
-      begin++;
-      src++;
+      ++begin;
+      ++src;
     }
 }
 
@@ -596,7 +592,7 @@ MappingCartesian<dim>::transform_unit_to_real_cell (
 
   Point<dim> p_real = cell->vertex(0);
   for (unsigned int d=0; d<dim; ++d)
-    p_real(d) +=length[d]*p(d);
+    p_real(d) += length[d]*p(d);
 
   return p_real;
 }
index 0b8542ac44a2e4f1a333c9cc293bbab21a983d6c..32c89cc490a9520e6d4273616936b04563caf975 100644 (file)
@@ -22,7 +22,6 @@
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary.h>
 #include <dofs/dof_accessor.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_tools.h>
 
 #include <numeric>
@@ -339,8 +338,9 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
   const unsigned int n_q_points=q.n_quadrature_points;
   this->compute_fill_face (cell, face_no, false,
                            n_q_points,
-                           MappingQ1<dim>::DataSetDescriptor::
-                           face (cell, face_no, n_q_points),
+                           QProjector<dim>::DataSetDescriptor::
+                           face (face_no, cell->face_orientation(face_no),
+                                 n_q_points),
                            q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
@@ -395,8 +395,10 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
   const unsigned int n_q_points=q.n_quadrature_points;
   this->compute_fill_face (cell, face_no, true,
                            n_q_points,
-                           MappingQ1<dim>::DataSetDescriptor::
-                           sub_face (cell, face_no, sub_no, n_q_points),
+                           QProjector<dim>::DataSetDescriptor::
+                           sub_face (face_no, sub_no,
+                                     cell->face_orientation(face_no),
+                                     n_q_points),
                            q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
index aaae80db024fa0dc0467490bac7613f14565febd..0c322fe21d99d341711473964bd7de50b008e79c 100644 (file)
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
-#include <fe/mapping_q1.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
 
 #include <cmath>
 #include <algorithm>
 
 
-namespace internal
-{
-  template <int dim>
-  DataSetDescriptor<dim>
-  DataSetDescriptor<dim>::cell ()
-  {
-    return 0;
-  }
-
-
-  template <int dim>
-  DataSetDescriptor<dim>
-  DataSetDescriptor<dim>::
-  face (const typename DoFHandler<dim>::cell_iterator &cell,
-        const unsigned int face_no,
-        const unsigned int n_quadrature_points)
-  {
-    Assert (dim != 1, ExcInternalError());
-    Assert (face_no < GeometryInfo<dim>::faces_per_cell,
-            ExcInternalError());
-    Assert (n_quadrature_points > 0, ExcInternalError());
-    
-    switch (dim)
-      {
-        case 1:
-        case 2:
-              return face_no * n_quadrature_points;
-
-                                               // in 3d, we have to
-                                               // account for faces
-                                               // that have reverse
-                                               // orientation. thus,
-                                               // we have to store
-                                               // _two_ data sets per
-                                               // face or subface
-        case 3:
-              return ((face_no +
-                       (cell->face_orientation(face_no) == true ?
-                        0 : GeometryInfo<dim>::faces_per_cell))
-                      * n_quadrature_points);
-
-        default:
-              Assert (false, ExcInternalError());
-      }
-    return static_cast<unsigned int>(-1);
-  }
-
-
-
-  template <int dim>
-  DataSetDescriptor<dim>
-  DataSetDescriptor<dim>::
-  sub_face (const typename DoFHandler<dim>::cell_iterator &cell,
-            const unsigned int face_no,
-            const unsigned int subface_no,
-            const unsigned int n_quadrature_points)
-  {
-    Assert (dim != 1, ExcInternalError());
-    Assert (face_no < GeometryInfo<dim>::faces_per_cell,
-            ExcInternalError());
-                                     // the trick with +1 prevents
-                                     // that we get a warning in 1d
-    Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
-            ExcInternalError());
-    Assert (n_quadrature_points > 0, ExcInternalError());
-
-    Assert (cell->has_children() == false,
-            ExcMessage ("You can't use subface data for cells that are "
-                        "already refined. Iterate over their children "
-                        "instead in these cases."));
-    
-    switch (dim)
-      {
-        case 1:
-        case 2:
-              return ((face_no * GeometryInfo<dim>::subfaces_per_face +
-                       subface_no)
-                      * n_quadrature_points);
-
-                                               // for 3d, same as above:
-        case 3:
-              return (((face_no * GeometryInfo<dim>::subfaces_per_face +
-                        subface_no)
-                       + (cell->face_orientation(face_no) == true ?
-                          0 :
-                          GeometryInfo<dim>::faces_per_cell *
-                          GeometryInfo<dim>::subfaces_per_face)
-                       )
-                      * n_quadrature_points);
-        default:
-              Assert (false, ExcInternalError());
-      }
-    return static_cast<unsigned int>(-1);              
-  }
-
-
-  template <int dim>
-  unsigned int
-  DataSetDescriptor<dim>::offset () const
-  {
-    return dataset_offset;
-  }
-
-
-
-  template <int dim>
-  DataSetDescriptor<dim>::
-  DataSetDescriptor (const unsigned int dataset_offset)
-                  :
-                  dataset_offset (dataset_offset)
-  {}
-
-
-  template <int dim>
-  DataSetDescriptor<dim>::
-  DataSetDescriptor ()
-                  :
-                  dataset_offset (static_cast<unsigned int>(-1))
-  {}
-}
-
-  
-  
 
 
 template <int dim>
@@ -625,7 +502,7 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
     for (unsigned int point=0; point<n_q_points; ++point)
       for (unsigned int k=0; k<data.n_shape_functions; ++k)
         quadrature_points[point]
-          += data.shape(point+data_set.offset(),k)
+          += data.shape(point+data_set,k)
              * data.mapping_support_points[k];
   
                                    // then Jacobians
@@ -635,7 +512,7 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
         for (unsigned int i=0; i<dim; ++i)
           for (unsigned int j=0; j<dim; ++j)
             data.contravariant[point][i][j]
-              += (data.derivative(point+data_set.offset(), k)[j]
+              += (data.derivative(point+data_set, k)[j]
                   *
                   data.mapping_support_points[k][i]);
 
@@ -811,7 +688,9 @@ MappingQ1<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterat
   
   compute_fill_face (cell, face_no, false,
                     n_q_points,
-                    DataSetDescriptor::face (cell, face_no, n_q_points),
+                    DataSetDescriptor::face (face_no,
+                                              cell->face_orientation(face_no),
+                                              n_q_points),
                     q.get_weights(),
                     data,
                     quadrature_points,
@@ -841,7 +720,9 @@ MappingQ1<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_ite
   
   compute_fill_face (cell, face_no, true,
                     n_q_points,
-                    DataSetDescriptor::sub_face (cell, face_no, sub_no, n_q_points),
+                    DataSetDescriptor::sub_face (face_no, sub_no,
+                                                  cell->face_orientation(face_no),
+                                                  n_q_points),
                     q.get_weights(),
                     data,
                     quadrature_points,
@@ -1163,7 +1044,3 @@ MappingQ1<dim>::transform_real_to_unit_cell_internal (
 //----------------------------------------------------------------------//
 
 template class MappingQ1<deal_II_dimension>;
-namespace internal
-{
-  template class DataSetDescriptor<deal_II_dimension>;
-}
index 907dc6edb85afc0505a24b2357d9ee77ebf4c030..aee011d7c9458b735791eef7700e04a4260c5af1 100644 (file)
@@ -802,11 +802,13 @@ integrate_over_regular_face (const DoFHandler<dim>               &dof_handler,
       
       const active_cell_iterator neighbor = cell->neighbor(face_no);
       
-                                      // find which number the current
-                                      // face has relative to the neighboring
-                                      // cell
-      const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-      Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
+                                      // find which number the
+                                      // current face has relative to
+                                      // the neighboring cell
+      const unsigned int neighbor_neighbor
+        = cell->neighbor_of_neighbor (face_no);
+      Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
+              ExcInternalError());
       
                                       // get restriction of finite element
                                       // function of @p{neighbor} to the
@@ -984,38 +986,24 @@ integrate_over_irregular_face (const DoFHandler<dim>               &dof_handler,
                                   // index the number of the
                                   // quadrature point
   
-                                  // store which number @p{cell} has in the
-                                  // list of neighbors of @p{neighbor}
-  const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
-  Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
+                                  // store which number @p{cell} has
+                                  // in the list of neighbors of
+                                  // @p{neighbor}
+  const unsigned int neighbor_neighbor
+    = cell->neighbor_of_neighbor (face_no);
+  Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell,
+          ExcInternalError());
   
                                   // loop over all subfaces
-  for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
+  for (unsigned int subface_no=0;
+       subface_no<GeometryInfo<dim>::subfaces_per_face;
        ++subface_no)
     {
                                       // get an iterator pointing to the
                                       // cell behind the present subface
-      static const unsigned int subface_translation[4]
-        = { 0, 3, 2, 1 };
-                                       // see whether face and
-                                       // the neighbor's
-                                       // counterface share the
-                                       // same indexing of
-                                       // children. if not so,
-                                       // translate child
-                                       // indices
-      const bool face_orientations_match
-        = (neighbor->face_orientation(neighbor_neighbor) ==
-           cell->face_orientation(face_no));
-      const unsigned int neighbor_child_index
-        = (GeometryInfo<dim>::
-           child_cell_on_face(neighbor_neighbor,
-                              (face_orientations_match ?
-                               subface_no :
-                               subface_translation[subface_no])));
       const active_cell_iterator neighbor_child
-       = neighbor->child(neighbor_child_index);
-      Assert (!neighbor->child(neighbor_child_index)->has_children(),
+        = cell->neighbor_child_on_subface (face_no, subface_no);
+      Assert (!neighbor_child->has_children(),
              ExcInternalError());
       
                                       // restrict the finite element
index 27c1df1674cf8bc30d5b7b9b16c741d43d218c45..9699185020b15a11be6eabf9a8c18fa2e30ad5e3 100644 (file)
@@ -32,15 +32,26 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3 style="color:red">Incompatibilities</h3>
 
 <ol>
+  <li> <p>
+       Changed: The <code class="class">QProjector</class> has functions that
+       project a lower-dimensional quadrature formula onto all faces or
+       subfaces of a cell. In 3d, it now does this but also adds projections of
+       these quadrature formula onto the faces from the other side. We need
+       this due to the fact that we now support faces in 3d that have a normal
+       vector opposite to the standard direction.
+       <br> 
+       (WB 2003/10/17)
+       </p>
+
   <li> <p> Moved and changed: The header file
-  <tt>include/numerics/dof_renumbering.h</tt> has been moved to the
-  directory <tt>include/dofs</tt>, where it logically
-  belongs. Furthermore, the sorting parameter of the function <code
-  class="class">DoFRenumbering</code>::<code
-  class="member">component_wise</code> has changed its meaning. See
-  the reference documentation for details.
-  <br>
-  (GK 2003/07/03)
+       <tt>include/numerics/dof_renumbering.h</tt> has been moved to the 
+       directory <tt>include/dofs</tt>, where it logically
+       belongs. Furthermore, the sorting parameter of the function <code
+       class="class">DoFRenumbering</code>::<code
+       class="member">component_wise</code> has changed its meaning. See
+       the reference documentation for details.
+       <br>
+       (GK 2003/07/03)
 </ol>
 
 
@@ -48,6 +59,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>General</h3>
 
 <ol>
+  <li> <p>
+       Augmented: The <code
+       class="member">GeometryInfo::child_cell_on_face</code> 
+       results in a result that might not be what you expect in 3d in some
+       cases. A warning has been added to the documentation, as well as a
+       reference to the function to call instead.
+       <br> 
+       (WB 2003/10/17)
+       </p>
+
   <li> <p> Fixed: The step-14 program had a bug in the rare case that
        there are more CPUs in a machine than there are cells. This is
        now fixed.
@@ -235,6 +256,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: There is now a function <code
+       class="member">CellAccessor::neighbor_child_on_subface</code> 
+       that should be called instead of using <code
+       class="member">GeometryInfo::child_cell_on_face</code> in most cases.
+       <br> 
+       (WB 2003/10/17)
+       </p>
+
   <li> <p>
        New: <code class="class">GridGenerator</code> has a new
        function <code class="member">subdivided_hyper_rectangle</code> 
index f6e8c1916aa826559a333aa4cad48a93432a869d..20c6f5fe0eb14fec5b6021b856c72db45f8ced4d 100644 (file)
@@ -959,9 +959,9 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // `behind' the
                                                       // current
                                                       // subface.
-                     typename DoFHandler<dim>::active_cell_iterator neighbor_child=
-                       neighbor->child(GeometryInfo<dim>::
-                                       child_cell_on_face(neighbor2,subface_no));
+                     typename DoFHandler<dim>::active_cell_iterator
+                        neighbor_child
+                        = cell->neighbor_child_on_subface (face_no, subface_no);
                      
                                                       // As these are
                                                       // quite
@@ -1123,9 +1123,10 @@ void DGMethod<dim>::assemble_system1 ()
                      const unsigned int neighbor_face_no=faceno_subfaceno.first,
                                      neighbor_subface_no=faceno_subfaceno.second;
 
-                     Assert (neighbor->neighbor(neighbor_face_no)
-                             ->child(GeometryInfo<dim>::child_cell_on_face(
-                               face_no,neighbor_subface_no)) == cell, ExcInternalError());
+                     Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+                                                                   neighbor_subface_no)
+                              == cell,
+                              ExcInternalError());
 
                                                       // Reinit the
                                                       // appropriate
@@ -1310,9 +1311,8 @@ void DGMethod<dim>::assemble_system2 ()
                       subface_no<GeometryInfo<dim>::subfaces_per_face;
                       ++subface_no)
                    {
-                     typename DoFHandler<dim>::cell_iterator neighbor_child=
-                       neighbor->child(GeometryInfo<dim>::child_cell_on_face(
-                         neighbor2,subface_no));
+                     typename DoFHandler<dim>::cell_iterator neighbor_child
+                        = cell->neighbor_child_on_subface (face_no, subface_no);
                      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
                              ExcInternalError());
                      Assert (!neighbor_child->has_children(), ExcInternalError());
index 2ca8e13af35a972654739f0bab05246ecec3f7fe..11a158a71431844769d980f75acc779a8290215e 100644 (file)
@@ -3557,9 +3557,7 @@ namespace LaplaceSolver
                                         // assertion will be removed
                                         // anyway.
        const active_cell_iterator neighbor_child
-         = neighbor->child(GeometryInfo<dim>::
-                           child_cell_on_face(neighbor_neighbor,
-                                              subface_no));
+          = cell->neighbor_child_on_subface (face_no, subface_no);
        Assert (neighbor_child->face(neighbor_neighbor) ==
                cell->face(face_no)->child(subface_no),
                ExcInternalError());

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.