]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Overhaul the way we compute offsets. This was a big mess and unintelligible. Now...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Sep 2003 14:48:29 +0000 (14:48 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Sep 2003 14:48:29 +0000 (14:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@8072 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/mapping_q1.h
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 c6e24aa1292a6fdd17440237ca7908b1eb892cc9..fc7080e53a850529b3a21488afff5f65b11a90c4 100644 (file)
 /*!@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);
+  };
+}
+
+
+
 /**
  * Mapping of general quadrilateral/hexahedra by d-linear shape
  * functions.
@@ -351,6 +385,12 @@ class MappingQ1 : public Mapping<dim>
     DeclException0 (ExcAccessToUninitializedField);
 
   protected:
+                                     /**
+                                      * Typedef the right data set
+                                      * descriptor to a local type to
+                                      * make use somewhat simpler.
+                                      */
+    typedef internal::DataSetDescriptor<dim> DataSetDescriptor;    
     
                                     /**
                                      * Implementation of the interface in
@@ -442,9 +482,9 @@ class MappingQ1 : public Mapping<dim>
                                      * @p{fill_*} functions.
                                      */
     void compute_fill (const typename DoFHandler<dim>::cell_iterator &cell,
-                      const unsigned int   npts,
-                      const unsigned int   offset,
-                      InternalData        &data,
+                      const unsigned int      npts,
+                      const DataSetDescriptor data_set,
+                      InternalData           &data,
                       std::vector<Point<dim> > &quadrature_points) const;
     
                                     /**
@@ -455,7 +495,7 @@ class MappingQ1 : public Mapping<dim>
                            const unsigned int      face_no,
                            const bool              is_subface,
                            const unsigned int      npts,
-                           const unsigned int      offset,
+                           const DataSetDescriptor data_set,
                            const std::vector<double>   &weights,
                            InternalData           &mapping_data,
                            std::vector<Point<dim> >    &quadrature_points,
@@ -683,7 +723,7 @@ template <> void MappingQ1<1>::compute_fill_face (
   const unsigned int,
   const bool,
   const unsigned int,
-  const unsigned int,
+  const DataSetDescriptor,
   const std::vector<double> &,
   InternalData &,
   std::vector<Point<1> > &,
index fcfc88ae3ab35711bd78cd4e401340451b5e0bfc..c7ebc882903f2cc0e9c4d328cc97cf65daf889f9 100644 (file)
@@ -19,6 +19,7 @@
 #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>
@@ -137,7 +138,9 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
   UpdateFlags update_flags(data.current_update_flags());
 
   const unsigned int npts = quadrature_points.size ();
-  unsigned int offset = 0;
+
+  typedef typename internal::DataSetDescriptor<dim> DataSetDescriptor;
+  unsigned int offset = DataSetDescriptor::cell().offset();
   
   if (face_no != invalid_face_number)
     {
@@ -150,14 +153,20 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
 
       if (sub_no == invalid_face_number)
                                         // called from FEFaceValues
-       offset = face_no * quadrature_points.size();
+       offset = DataSetDescriptor::face (cell, face_no,
+                                          quadrature_points.size()).offset();
       else
        {
-                                          // called from FESubfaceValue
+                                          // 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 = (face_no * GeometryInfo<dim>::subfaces_per_face + sub_no)
-                  * quadrature_points.size();
+         offset = DataSetDescriptor::sub_face (cell, face_no, sub_no,
+                                                quadrature_points.size()).offset();
        }
     }
   else
index 47543ddd12ee7d5b8cd3d3a60096c86a0d86edf7..0b8542ac44a2e4f1a333c9cc293bbab21a983d6c 100644 (file)
@@ -226,7 +226,7 @@ MappingQ<dim>::get_face_data (const UpdateFlags update_flags,
                              const Quadrature<dim-1>& quadrature) const
 {
   InternalData *data = new InternalData(n_shape_functions);
-  Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature));
+  const Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature));
   this->compute_face_data (update_flags, q,
                            quadrature.n_quadrature_points, *data);
   if (!use_mapping_q_on_all_cells)
@@ -244,7 +244,7 @@ MappingQ<dim>::get_subface_data (const UpdateFlags update_flags,
                                 const Quadrature<dim-1>& quadrature) const
 {
   InternalData *data = new InternalData(n_shape_functions);
-  Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature));
+  const Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature));
   this->compute_face_data (update_flags, q,
                            quadrature.n_quadrature_points, *data);
   if (!use_mapping_q_on_all_cells)
@@ -327,9 +327,6 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                        || cell->has_boundary_lines());
 
-  const unsigned int npts=q.n_quadrature_points;
-  const unsigned int offset=face_no*npts;
-
                                   // depending on this result, use
                                   // this or the other data object
                                   // for the mapping
@@ -339,8 +336,12 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
   else
     p_data=&data;
 
+  const unsigned int n_q_points=q.n_quadrature_points;
   this->compute_fill_face (cell, face_no, false,
-                           npts, offset, q.get_weights(),
+                           n_q_points,
+                           MappingQ1<dim>::DataSetDescriptor::
+                           face (cell, face_no, n_q_points),
+                           q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
                            exterior_forms, normal_vectors);
@@ -382,10 +383,6 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                        || cell->has_boundary_lines());
 
-  const unsigned int npts=q.n_quadrature_points;
-  const unsigned int offset=
-    (face_no*GeometryInfo<dim>::subfaces_per_face + sub_no)*npts;
-
                                   // depending on this result, use
                                   // this or the other data object
                                   // for the mapping
@@ -395,8 +392,12 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
   else
     p_data=&data;
 
+  const unsigned int n_q_points=q.n_quadrature_points;
   this->compute_fill_face (cell, face_no, true,
-                           npts, offset, q.get_weights(),
+                           n_q_points,
+                           MappingQ1<dim>::DataSetDescriptor::
+                           sub_face (cell, face_no, sub_no, n_q_points),
+                           q.get_weights(),
                            *p_data,
                            quadrature_points, JxW_values,
                            exterior_forms, normal_vectors);
index c342fd7f8acb3f28a4e28492482aa3eeacf1ff0d..0140474a9753b8c4886e3ec6f4cf4bdac512ce2f 100644 (file)
 #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)
+  {
+    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->get_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)
+  {
+    switch (dim)
+      {
+        case 1:
+        case 2:
+              return ((face_no * GeometryInfo<dim>::faces_per_cell +
+                       subface_no)
+                      * n_quadrature_points);
+
+                                               // for 3d, same as above:
+        case 3:
+              return (((face_no * GeometryInfo<dim>::faces_per_cell +
+                        subface_no)
+                       + (cell->get_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>
 const unsigned int MappingQ1<dim>::n_shape_functions;
@@ -291,7 +395,7 @@ MappingQ1<dim>::compute_data (const UpdateFlags      update_flags,
                              const unsigned int     n_original_q_points,
                              InternalData          &data) const
 {
-  const unsigned int npts = q.n_quadrature_points;
+  const unsigned int n_q_points = q.n_quadrature_points;
 
   data.update_once = update_once(update_flags);
   data.update_each = update_each(update_flags);
@@ -300,10 +404,10 @@ MappingQ1<dim>::compute_data (const UpdateFlags      update_flags,
   const UpdateFlags flags(data.update_flags);
   
   if (flags & update_transformation_values)
-    data.shape_values.resize(data.n_shape_functions * npts);
+    data.shape_values.resize(data.n_shape_functions * n_q_points);
 
   if (flags & update_transformation_gradients)
-    data.shape_derivatives.resize(data.n_shape_functions * npts);
+    data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
 
   if (flags & update_covariant_transformation)
     data.covariant.resize(n_original_q_points);
@@ -438,17 +542,17 @@ MappingQ1<dim>::get_subface_data (const UpdateFlags update_flags,
 template <int dim>
 void
 MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cell,
-                             const unsigned int   npts,
-                             const unsigned int   offset,
-                             InternalData        &data,
+                             const unsigned int       n_q_points,
+                             const DataSetDescriptor  data_set,
+                             InternalData            &data,
                              std::vector<Point<dim> > &quadrature_points) const
 {
   const UpdateFlags update_flags(data.current_update_flags());
 
   if (update_flags & update_q_points)
     {
-      Assert (quadrature_points.size() == npts,
-             ExcDimensionMismatch(quadrature_points.size(), npts));
+      Assert (quadrature_points.size() == n_q_points,
+             ExcDimensionMismatch(quadrature_points.size(), n_q_points));
       std::fill(quadrature_points.begin(),
                quadrature_points.end(),
                Point<dim>());
@@ -456,14 +560,14 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
 
   if (update_flags & update_covariant_transformation)
     {
-      Assert (data.covariant.size() == npts,
-             ExcDimensionMismatch(data.covariant.size(), npts));
+      Assert (data.covariant.size() == n_q_points,
+             ExcDimensionMismatch(data.covariant.size(), n_q_points));
     }
   
   if (update_flags & update_contravariant_transformation)
     {
-      Assert (data.contravariant.size() == npts,
-             ExcDimensionMismatch(data.contravariant.size(), npts));
+      Assert (data.contravariant.size() == n_q_points,
+             ExcDimensionMismatch(data.contravariant.size(), n_q_points));
       std::fill(data.contravariant.begin(),
                data.contravariant.end(),
                Tensor<2,dim>());
@@ -472,8 +576,8 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
   if (update_flags & update_jacobian_grads)
     {
       Assert(false, ExcNotImplemented());
-//      Assert (covariant_grads.size () == npts,
-//           ExcDimensionMismatch(covariant_grads.size(), npts));
+//      Assert (covariant_grads.size () == n_q_points,
+//           ExcDimensionMismatch(covariant_grads.size(), n_q_points));
     }
   
                                   // if necessary, recompute the
@@ -498,19 +602,20 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
 
                                    // first compute quadrature points
   if (update_flags & update_q_points)
-    for (unsigned int point=0; point<npts; ++point)
+    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+offset,k) * data.mapping_support_points[k];
+          += data.shape(point+data_set.offset(),k)
+             * data.mapping_support_points[k];
   
                                    // then Jacobians
   if (update_flags & update_contravariant_transformation)
-    for (unsigned int point=0; point<npts; ++point)
+    for (unsigned int point=0; point<n_q_points; ++point)
       for (unsigned int k=0; k<data.n_shape_functions; ++k)
         for (unsigned int i=0; i<dim; ++i)
           for (unsigned int j=0; j<dim; ++j)
             data.contravariant[point][i][j]
-              += (data.derivative(point+offset, k)[j]
+              += (data.derivative(point+data_set.offset(), k)[j]
                   *
                   data.mapping_support_points[k][i]);
 
@@ -518,7 +623,7 @@ MappingQ1<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cel
                                       // covariant transformation
                                       // matrices
   if (update_flags & update_covariant_transformation)
-    for (unsigned int point=0; point<npts; ++point)
+    for (unsigned int point=0; point<n_q_points; ++point)
       data.covariant[point] = invert(data.contravariant[point]);
 }
 
@@ -549,13 +654,10 @@ MappingQ1<dim>::fill_fe_values (const typename DoFHandler<dim>::cell_iterator &c
   Assert(data_ptr!=0, ExcInternalError());
   InternalData &data=*data_ptr;
 
-  const unsigned int npts=q.n_quadrature_points;
+  const unsigned int n_q_points=q.n_quadrature_points;
   
-  compute_fill (cell,
-               npts,
-               0,
-               data,
-               quadrature_points);
+  compute_fill (cell, n_q_points, DataSetDescriptor::cell (),
+                data, quadrature_points);
 
   
   const UpdateFlags update_flags(data.current_update_flags());
@@ -565,9 +667,9 @@ MappingQ1<dim>::fill_fe_values (const typename DoFHandler<dim>::cell_iterator &c
                                   // by Jacobian determinants
   if (update_flags & update_JxW_values)
     {      
-      Assert (JxW_values.size() == npts,
-             ExcDimensionMismatch(JxW_values.size(), npts));
-      for (unsigned int point=0; point<npts; ++point)
+      Assert (JxW_values.size() == n_q_points,
+             ExcDimensionMismatch(JxW_values.size(), n_q_points));
+      for (unsigned int point=0; point<n_q_points; ++point)
        JxW_values[point]
          = determinant(data.contravariant[point])*weights[point];
     }
@@ -580,8 +682,8 @@ void
 MappingQ1<dim>::compute_fill_face (const typename DoFHandler<dim>::cell_iterator &cell,
                                   const unsigned int      face_no,
                                   const bool              is_subface,
-                                  const unsigned int      npts,
-                                  const unsigned int      offset,
+                                  const unsigned int      n_q_points,
+                                  const DataSetDescriptor data_set,
                                   const std::vector<double>   &weights,
                                   InternalData           &data,
                                   std::vector<Point<dim> >    &quadrature_points,
@@ -589,24 +691,20 @@ MappingQ1<dim>::compute_fill_face (const typename DoFHandler<dim>::cell_iterator
                                   std::vector<Tensor<1,dim> > &boundary_forms,
                                   std::vector<Point<dim> >    &normal_vectors) const
 {
-  compute_fill (cell,
-               npts,
-               offset,
-               data,
-               quadrature_points);
+  compute_fill (cell, n_q_points, data_set, data, quadrature_points);
 
   const UpdateFlags update_flags(data.current_update_flags());
   
   if (update_flags & update_boundary_forms)
     {
-      Assert (boundary_forms.size()==npts,
-             ExcDimensionMismatch(boundary_forms.size(), npts));
+      Assert (boundary_forms.size()==n_q_points,
+             ExcDimensionMismatch(boundary_forms.size(), n_q_points));
       if (update_flags & update_normal_vectors)
-       Assert (normal_vectors.size()==npts,
-               ExcDimensionMismatch(normal_vectors.size(), npts));
+       Assert (normal_vectors.size()==n_q_points,
+               ExcDimensionMismatch(normal_vectors.size(), n_q_points));
       if (update_flags & update_JxW_values)
-       Assert (JxW_values.size() == npts,
-               ExcDimensionMismatch(JxW_values.size(), npts));
+       Assert (JxW_values.size() == n_q_points,
+               ExcDimensionMismatch(JxW_values.size(), n_q_points));
       
       
       Assert (data.aux[0].size() <= data.unit_tangentials[face_no].size(),
@@ -665,7 +763,7 @@ MappingQ1<dim>::compute_fill_face (const typename DoFHandler<dim>::cell_iterator
              {
                JxW_values[i] = f * weights[i];
                if (is_subface)
-                 JxW_values[i]/=GeometryInfo<dim>::subfaces_per_face;
+                 JxW_values[i] /= GeometryInfo<dim>::subfaces_per_face;
              }
            if (update_flags & update_normal_vectors)
               normal_vectors[i] = boundary_forms[i] / f;
@@ -689,12 +787,11 @@ MappingQ1<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterat
   Assert(data_ptr!=0, ExcInternalError());
   InternalData &data=*data_ptr;
 
-  const unsigned int npts=q.n_quadrature_points;
-  const unsigned int offset=face_no*npts;
+  const unsigned int n_q_points = q.n_quadrature_points;
   
   compute_fill_face (cell, face_no, false,
-                    npts,
-                    offset,
+                    n_q_points,
+                    DataSetDescriptor::face (cell, face_no, n_q_points),
                     q.get_weights(),
                     data,
                     quadrature_points,
@@ -720,13 +817,11 @@ MappingQ1<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_ite
   Assert(data_ptr!=0, ExcInternalError());
   InternalData &data=*data_ptr;
 
-  const unsigned int npts=q.n_quadrature_points;
-  const unsigned int offset=
-    (face_no*GeometryInfo<dim>::subfaces_per_face + sub_no)*npts;
+  const unsigned int n_q_points = q.n_quadrature_points;
   
   compute_fill_face (cell, face_no, true,
-                    npts,
-                    offset,
+                    n_q_points,
+                    DataSetDescriptor::sub_face (cell, face_no, sub_no, n_q_points),
                     q.get_weights(),
                     data,
                     quadrature_points,
@@ -744,7 +839,7 @@ MappingQ1<1>::compute_fill_face (const DoFHandler<1>::cell_iterator &,
                                 const unsigned int,
                                 const bool,
                                 const unsigned int,
-                                const unsigned int,
+                                const DataSetDescriptor,
                                 const std::vector<double> &,
                                 InternalData &,
                                 std::vector<Point<1> > &,

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.