]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove MappingQ1::compute_fill and simplify MappingQ1::compute_fill_face.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 27 Jul 2015 22:06:51 +0000 (17:06 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Jul 2015 02:43:53 +0000 (21:43 -0500)
Also move around a few functions in the .cc file to make the flow easier
to read. The old style used rather generic function names that don't
really give away what concretely a function does.

include/deal.II/fe/mapping_q1.h
source/fe/mapping_q.cc
source/fe/mapping_q1.cc

index f2b6f0524250f33989725b01b90f779ccc462ac4..81071b09a710cadb45c71df108f4abdb9b3332f5 100644 (file)
@@ -386,25 +386,14 @@ public:
                           const unsigned int n_orig_q_points,
                           InternalData &data) const;
 
-  /**
-   * Do the computation for the <tt>fill_*</tt> functions.
-   */
-  void compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                     const unsigned int      npts,
-                     const DataSetDescriptor data_set,
-                     const CellSimilarity::Similarity cell_similarity,
-                     const InternalData           &data,
-                     std::vector<Point<spacedim> > &quadrature_points) const;
-
   /**
    * Do the computation for the <tt>fill_*</tt> functions.
    */
   void compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                           const unsigned int                face_no,
                           const unsigned int                subface_no,
-                          const unsigned int                npts,
                           const DataSetDescriptor           data_set,
-                          const std::vector<double>        &weights,
+                          const Quadrature<dim-1>          &quadrature,
                           const InternalData                     &internal_data,
                           FEValuesData<dim,spacedim>                                &output_data) const;
 
index 034c0540c83c1813054cd7b7c6c646a105679e9f..edb1827ae4277ff1f6eb3c7867f6364b2b81471d 100644 (file)
@@ -334,14 +334,13 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
 
   const unsigned int n_q_points = quadrature.size();
   this->compute_fill_face (cell, face_no, numbers::invalid_unsigned_int,
-                           n_q_points,
                            QProjector<dim>::DataSetDescriptor::
                            face (face_no,
                                  cell->face_orientation(face_no),
                                  cell->face_flip(face_no),
                                  cell->face_rotation(face_no),
                                  n_q_points),
-                           quadrature.get_weights(),
+                           quadrature,
                            *p_data,
                            output_data);
 }
@@ -383,7 +382,6 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 
   const unsigned int n_q_points = quadrature.size();
   this->compute_fill_face (cell, face_no, subface_no,
-                           n_q_points,
                            QProjector<dim>::DataSetDescriptor::
                            subface (face_no, subface_no,
                                     cell->face_orientation(face_no),
@@ -391,7 +389,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
                                     cell->face_rotation(face_no),
                                     n_q_points,
                                     cell->subface_case(face_no)),
-                           quadrature.get_weights(),
+                           quadrature,
                            *p_data,
                            output_data);
 }
index b5bc80917d55d5ced5d39c1e3f975adda03373ed..7e090cfb0826796ab24373b738d80923a7e3ebd7 100644 (file)
@@ -603,121 +603,122 @@ MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
 
 
 
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                                       const unsigned int  n_q_points,
-                                       const DataSetDescriptor  data_set,
-                                       const CellSimilarity::Similarity cell_similarity,
-                                       const InternalData  &data,
-                                       std::vector<Point<spacedim> > &quadrature_points) const
+namespace internal
 {
-  const UpdateFlags update_flags(data.current_update_flags());
-
-  // if necessary, recompute the
-  // support points of the
-  // transformation of this cell
-  // (note that we need to first
-  // check the triangulation pointer,
-  // since otherwise the second test
-  // might trigger an exception if
-  // the triangulations are not the
-  // same)
-  if ((data.mapping_support_points.size() == 0)
-      ||
-      (&cell->get_triangulation() !=
-       &data.cell_of_current_support_points->get_triangulation())
-      ||
-      (cell != data.cell_of_current_support_points))
-    {
-      compute_mapping_support_points(cell, data.mapping_support_points);
-      data.cell_of_current_support_points = cell;
-    }
-
-  // first compute quadrature points
-  if (update_flags & update_quadrature_points)
-    {
-      AssertDimension (quadrature_points.size(), n_q_points);
-
-      for (unsigned int point=0; point<n_q_points; ++point)
-        {
-          const double *shape = &data.shape(point+data_set,0);
-          Point<spacedim> result = (shape[0] *
-                                    data.mapping_support_points[0]);
-          for (unsigned int k=1; k<data.n_shape_functions; ++k)
-            for (unsigned int i=0; i<spacedim; ++i)
-              result[i] += shape[k] * data.mapping_support_points[k][i];
-          quadrature_points[point] = result;
-        }
-    }
-
-  // then Jacobians
-  if (update_flags & update_contravariant_transformation)
+  namespace
+  {
+    /**
+     * Compute the locations of quadrature points on the object described by
+     * the first argument (and the cell for which the mapping support points
+     * have already been set), but only if the update_flags of the @p data
+     * argument indicate so.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_compute_q_points (const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor  data_set,
+                            const typename dealii::MappingQ1<dim,spacedim>::InternalData      &data,
+                            std::vector<Point<spacedim> >                                     &quadrature_points)
     {
-      AssertDimension (data.contravariant.size(), n_q_points);
+      const UpdateFlags update_flags(data.current_update_flags());
 
-      // if the current cell is just a
-      // translation of the previous one, no
-      // need to recompute jacobians...
-      if (cell_similarity != CellSimilarity::translation)
+      if (update_flags & update_quadrature_points)
         {
-          std::fill(data.contravariant.begin(), data.contravariant.end(),
-                    DerivativeForm<1,dim,spacedim>());
-
-          Assert (data.n_shape_functions > 0, ExcInternalError());
-          const Tensor<1,spacedim> *supp_pts =
-            &data.mapping_support_points[0];
-
-          for (unsigned int point=0; point<n_q_points; ++point)
+          for (unsigned int point=0; point<quadrature_points.size(); ++point)
             {
-              const Tensor<1,dim> *data_derv =
-                &data.derivative(point+data_set, 0);
-
-              double result [spacedim][dim];
-
-              // peel away part of sum to avoid zeroing the
-              // entries and adding for the first time
-              for (unsigned int i=0; i<spacedim; ++i)
-                for (unsigned int j=0; j<dim; ++j)
-                  result[i][j] = data_derv[0][j] * supp_pts[0][i];
+              const double *shape = &data.shape(point+data_set,0);
+              Point<spacedim> result = (shape[0] *
+                                        data.mapping_support_points[0]);
               for (unsigned int k=1; k<data.n_shape_functions; ++k)
                 for (unsigned int i=0; i<spacedim; ++i)
-                  for (unsigned int j=0; j<dim; ++j)
-                    result[i][j] += data_derv[k][j] * supp_pts[k][i];
-
-              // write result into contravariant data. for
-              // j=dim in the case dim<spacedim, there will
-              // never be any nonzero data that arrives in
-              // here, so it is ok anyway because it was
-              // initialized to zero at the initialization
-              for (unsigned int i=0; i<spacedim; ++i)
-                for (unsigned int j=0; j<dim; ++j)
-                  data.contravariant[point][i][j] = result[i][j];
+                  result[i] += shape[k] * data.mapping_support_points[k][i];
+              quadrature_points[point] = result;
             }
         }
     }
 
-  if (update_flags & update_covariant_transformation)
-    {
-      AssertDimension (data.covariant.size(), n_q_points);
-      if (cell_similarity != CellSimilarity::translation)
-        for (unsigned int point=0; point<n_q_points; ++point)
-          {
-            data.covariant[point] = (data.contravariant[point]).covariant_form();
-          }
-    }
 
-  if (update_flags & update_volume_elements)
-    if (cell_similarity != CellSimilarity::translation)
-      for (unsigned int point=0; point<n_q_points; ++point)
-        data.volume_elements[point] = data.contravariant[point].determinant();
+    /**
+     * Update the co- and contravariant matrices as well as their determinant, for the cell
+     * described stored in the data object, but only if the update_flags of the @p data
+     * argument indicate so.
+     *
+     * Skip the computation if possible as indicated by the first argument.
+     */
+    template <int dim, int spacedim>
+    void
+    maybe_update_Jacobians (const CellSimilarity::Similarity                                   cell_similarity,
+                            const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor  data_set,
+                            const typename dealii::MappingQ1<dim,spacedim>::InternalData      &data)
+    {
+      const UpdateFlags update_flags(data.current_update_flags());
 
+      if (update_flags & update_contravariant_transformation)
+       // if the current cell is just a
+       // translation of the previous one, no
+       // need to recompute jacobians...
+       if (cell_similarity != CellSimilarity::translation)
+         {
+           const unsigned int n_q_points = data.contravariant.size();
+
+           std::fill(data.contravariant.begin(), data.contravariant.end(),
+                     DerivativeForm<1,dim,spacedim>());
+
+           Assert (data.n_shape_functions > 0, ExcInternalError());
+           const Tensor<1,spacedim> *supp_pts =
+             &data.mapping_support_points[0];
+
+           for (unsigned int point=0; point<n_q_points; ++point)
+             {
+               const Tensor<1,dim> *data_derv =
+                 &data.derivative(point+data_set, 0);
+
+               double result [spacedim][dim];
+
+               // peel away part of sum to avoid zeroing the
+               // entries and adding for the first time
+               for (unsigned int i=0; i<spacedim; ++i)
+                 for (unsigned int j=0; j<dim; ++j)
+                   result[i][j] = data_derv[0][j] * supp_pts[0][i];
+               for (unsigned int k=1; k<data.n_shape_functions; ++k)
+                 for (unsigned int i=0; i<spacedim; ++i)
+                   for (unsigned int j=0; j<dim; ++j)
+                     result[i][j] += data_derv[k][j] * supp_pts[k][i];
+
+               // write result into contravariant data. for
+               // j=dim in the case dim<spacedim, there will
+               // never be any nonzero data that arrives in
+               // here, so it is ok anyway because it was
+               // initialized to zero at the initialization
+               for (unsigned int i=0; i<spacedim; ++i)
+                 for (unsigned int j=0; j<dim; ++j)
+                   data.contravariant[point][i][j] = result[i][j];
+             }
+         }
+
+      if (update_flags & update_covariant_transformation)
+       if (cell_similarity != CellSimilarity::translation)
+         {
+           const unsigned int n_q_points = data.contravariant.size();
+            for (unsigned int point=0; point<n_q_points; ++point)
+              {
+                data.covariant[point] = (data.contravariant[point]).covariant_form();
+              }
+        }
 
+      if (update_flags & update_volume_elements)
+       if (cell_similarity != CellSimilarity::translation)
+         {
+           const unsigned int n_q_points = data.contravariant.size();    
+           for (unsigned int point=0; point<n_q_points; ++point)
+             data.volume_elements[point] = data.contravariant[point].determinant();
+       }
+      
+    }
+  }
 }
 
 
 
-
 template<int dim, int spacedim>
 void
 MappingQ1<dim,spacedim>::compute_mapping_support_points(
@@ -748,10 +749,27 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 
   const unsigned int n_q_points=quadrature.size();
 
-  compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity,
-                data,
-                output_data.quadrature_points);
+  // if necessary, recompute the support points of the transformation of this cell
+  // (note that we need to first check the triangulation pointer, since otherwise
+  // the second test might trigger an exception if the triangulations are not the
+  // same)
+  if ((data.mapping_support_points.size() == 0)
+      ||
+      (&cell->get_triangulation() !=
+       &data.cell_of_current_support_points->get_triangulation())
+      ||
+      (cell != data.cell_of_current_support_points))
+    {
+      compute_mapping_support_points(cell, data.mapping_support_points);
+      data.cell_of_current_support_points = cell;
+    }
 
+  internal::maybe_compute_q_points<dim,spacedim> (DataSetDescriptor::cell (),
+                                                  data,
+                                                  output_data.quadrature_points);
+  internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
+                                                  DataSetDescriptor::cell (),
+                                                  data);
 
   const UpdateFlags update_flags(data.current_update_flags());
   const std::vector<double> &weights=quadrature.get_weights();
@@ -915,16 +933,25 @@ namespace internal
 {
   namespace
   {
+    /**
+     * Depending on what information is called for in the update flags of the
+     * @p data object, compute the various pieces of information that is required
+     * by the fill_fe_face_values() and fill_fe_subface_values() functions.
+     * This function simply unifies the work that would be done by
+     * those two functions.
+     *
+     * The resulting data is put into the @p output_data argument.
+     */
     template <int dim, int spacedim>
     void
-    compute_fill_face (const dealii::MappingQ1<dim,spacedim> &mapping,
-                       const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
-                       const unsigned int               face_no,
-                       const unsigned int               subface_no,
-                       const unsigned int               n_q_points,
-                       const std::vector<double>        &weights,
-                       const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
-                       FEValuesData<dim,spacedim>                                   &output_data)
+    maybe_compute_face_data (const dealii::MappingQ1<dim,spacedim> &mapping,
+                             const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
+                             const unsigned int               face_no,
+                             const unsigned int               subface_no,
+                             const unsigned int               n_q_points,
+                             const std::vector<double>        &weights,
+                             const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
+                             FEValuesData<dim,spacedim>                                   &output_data)
     {
       const UpdateFlags update_flags(data.current_update_flags());
 
@@ -1060,21 +1087,38 @@ template<int dim, int spacedim>
 void
 MappingQ1<dim,spacedim>::
 compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                   const unsigned int               face_no,
-                   const unsigned int               subface_no,
-                   const unsigned int               n_q_points,
-                   const DataSetDescriptor          data_set,
-                   const std::vector<double>        &weights,
+                   const unsigned int                face_no,
+                   const unsigned int                subface_no,
+                   const DataSetDescriptor           data_set,
+                   const Quadrature<dim-1>           &quadrature,
                    const InternalData                     &internal_data,
                    FEValuesData<dim,spacedim>                                &output_data) const
 {
-  compute_fill (cell, n_q_points, data_set, CellSimilarity::none,
-                internal_data,
-                output_data.quadrature_points);
-  internal::compute_fill_face (*this,
-                               cell, face_no, subface_no, n_q_points,
-                               weights, internal_data,
-                               output_data);
+  // if necessary, recompute the support points of the transformation of this cell
+  // (note that we need to first check the triangulation pointer, since otherwise
+  // the second test might trigger an exception if the triangulations are not the
+  // same)
+  if ((internal_data.mapping_support_points.size() == 0)
+      ||
+      (&cell->get_triangulation() !=
+       &internal_data.cell_of_current_support_points->get_triangulation())
+      ||
+      (cell != internal_data.cell_of_current_support_points))
+    {
+      compute_mapping_support_points(cell, internal_data.mapping_support_points);
+      internal_data.cell_of_current_support_points = cell;
+    }
+
+  internal::maybe_compute_q_points<dim,spacedim> (data_set,
+                                                  internal_data,
+                                                  output_data.quadrature_points);
+  internal::maybe_update_Jacobians<dim,spacedim> (CellSimilarity::none,
+                                                  data_set,
+                                                  internal_data);
+  internal::maybe_compute_face_data (*this,
+                                     cell, face_no, subface_no, quadrature.size(),
+                                     quadrature.get_weights(), internal_data,
+                                     output_data);
 }
 
 
@@ -1096,13 +1140,12 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
   const unsigned int n_q_points = quadrature.size();
 
   compute_fill_face (cell, face_no, numbers::invalid_unsigned_int,
-                     n_q_points,
                      DataSetDescriptor::face (face_no,
                                               cell->face_orientation(face_no),
                                               cell->face_flip(face_no),
                                               cell->face_rotation(face_no),
                                               n_q_points),
-                     quadrature.get_weights(),
+                     quadrature,
                      data,
                      output_data);
 }
@@ -1128,14 +1171,13 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
   const unsigned int n_q_points = quadrature.size();
 
   compute_fill_face (cell, face_no, subface_no,
-                     n_q_points,
                      DataSetDescriptor::subface (face_no, subface_no,
                                                  cell->face_orientation(face_no),
                                                  cell->face_flip(face_no),
                                                  cell->face_rotation(face_no),
                                                  n_q_points,
                                                  cell->subface_case(face_no)),
-                     quadrature.get_weights(),
+                     quadrature,
                      data,
                      output_data);
 }

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.