]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split FEValuesData into two base classes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 21 Jul 2015 14:06:28 +0000 (09:06 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 31 Jul 2015 15:21:49 +0000 (10:21 -0500)
There is now one base class that stores Mapping-related data, and one
that stores FiniteElement-related data. These are declared in
fe_update_flags() so that they can be referenced from the finite
element and mapping classes without having to know about fe_values.h.

include/deal.II/fe/fe_update_flags.h
include/deal.II/fe/fe_values.h
include/deal.II/matrix_free/mapping_data_on_the_fly.h
include/deal.II/matrix_free/mapping_info.h
include/deal.II/matrix_free/mapping_info.templates.h
source/fe/fe_values.cc

index 0a934f287f6f306d63f713f49fe65d7c596b4608..4622f309a3a1b148998a3f4ce6bad213b234e7e9 100644 (file)
 
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <vector>
+
 
 DEAL_II_NAMESPACE_OPEN
 
+template <int,int> class FiniteElement;
+
+
 /*!@addtogroup feaccess */
 /*@{*/
 
@@ -336,6 +346,187 @@ namespace CellSimilarity
 }
 
 
+namespace internal
+{
+  namespace FEValues
+  {
+    /**
+     * A class that stores all of the mapping related data used in
+     * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues
+     * objects. Objects of this kind will be given
+     * as <i>output</i> argument when dealii::FEValues::reinit()
+     * calls Mapping::fill_fe_values() for a given cell, face, or subface.
+     *
+     * The data herein will then be provided as <i>input</i> argument in
+     * the following call to FiniteElement::fill_fe_values().
+     *
+     * @ingroup feaccess
+     */
+    template <int dim, int spacedim=dim>
+    class MappingRelatedData
+    {
+    public:
+      /**
+       * Initialize all vectors to correct size.
+       */
+      void initialize (const unsigned int n_quadrature_points,
+                       const UpdateFlags         flags);
+
+      /**
+       * Store an array of weights times the Jacobi determinant at the quadrature
+       * points. This function is reset each time reinit() is called. The Jacobi
+       * determinant is actually the reciprocal value of the Jacobi matrices
+       * stored in this class, see the general documentation of this class for
+       * more information.
+       *
+       * However, if this object refers to an FEFaceValues or FESubfaceValues
+       * object, then the JxW_values correspond to the Jacobian of the
+       * transformation of the face, not the cell, i.e. the dimensionality is that
+       * of a surface measure, not of a volume measure. In this case, it is
+       * computed from the boundary forms, rather than the Jacobian matrix.
+       */
+      std::vector<double>       JxW_values;
+
+      /**
+       * Array of the Jacobian matrices at the quadrature points.
+       */
+      std::vector< DerivativeForm<1,dim,spacedim> > jacobians;
+
+      /**
+       * Array of the derivatives of the Jacobian matrices at the quadrature
+       * points.
+       */
+      std::vector<DerivativeForm<2,dim,spacedim> >  jacobian_grads;
+
+      /**
+       * Array of the inverse Jacobian matrices at the quadrature points.
+       */
+      std::vector<DerivativeForm<1,spacedim,dim> > inverse_jacobians;
+
+      /**
+       * Array of quadrature points. This array is set up upon calling reinit()
+       * and contains the quadrature points on the real element, rather than on
+       * the reference element.
+       */
+      std::vector<Point<spacedim> >  quadrature_points;
+
+      /**
+       * List of outward normal vectors at the quadrature points. This field is
+       * filled in by the finite element class.
+       */
+      std::vector<Point<spacedim> >  normal_vectors;
+
+      /**
+       * List of boundary forms at the quadrature points. This field is filled in
+       * by the finite element class.
+       */
+      std::vector<Tensor<1,spacedim> >  boundary_forms;
+    };
+
+
+    /**
+     * A class that stores all of the shape function related data used in
+     * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues
+     * objects. Objects of this kind will be given
+     * as <i>output</i> argument when dealii::FEValues::reinit()
+     * calls FiniteElement::fill_fe_values().
+     *
+     * @ingroup feaccess
+     */
+    template <int dim, int spacedim=dim>
+    class FiniteElementRelatedData
+    {
+    public:
+      /**
+       * Initialize all vectors to correct size.
+       */
+      void initialize (const unsigned int        n_quadrature_points,
+                       const FiniteElement<dim,spacedim> &fe,
+                       const UpdateFlags         flags);
+
+      /**
+       * Storage type for shape values. Each row in the matrix denotes the values
+       * of a single shape function at the different points, columns are for a
+       * single point with the different shape functions.
+       *
+       * If a shape function has more than one non-zero component (in deal.II
+       * diction: it is non-primitive), then we allocate one row per non-zero
+       * component, and shift subsequent rows backward.  Lookup of the correct row
+       * for a shape function is thus simple in case the entire finite element is
+       * primitive (i.e. all shape functions are primitive), since then the shape
+       * function number equals the row number. Otherwise, use the
+       * #shape_function_to_row_table array to get at the first row that belongs
+       * to this particular shape function, and navigate among all the rows for
+       * this shape function using the FiniteElement::get_nonzero_components()
+       * function which tells us which components are non-zero and thus have a row
+       * in the array presently under discussion.
+       */
+      typedef dealii::Table<2,double> ShapeVector;
+
+      /**
+       * Storage type for gradients. The layout of data is the same as for the
+       * #ShapeVector data type.
+       */
+      typedef std::vector<std::vector<Tensor<1,spacedim> > > GradientVector;
+
+      /**
+       * Likewise for second order derivatives.
+       */
+      typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
+
+      /**
+       * Store the values of the shape functions at the quadrature points. See the
+       * description of the data type for the layout of the data in this field.
+       */
+      ShapeVector shape_values;
+
+      /**
+       * Store the gradients of the shape functions at the quadrature points. See
+       * the description of the data type for the layout of the data in this
+       * field.
+       */
+      GradientVector shape_gradients;
+
+      /**
+       * Store the 2nd derivatives of the shape functions at the quadrature
+       * points.  See the description of the data type for the layout of the data
+       * in this field.
+       */
+      HessianVector shape_hessians;
+
+      /**
+       * When asked for the value (or gradient, or Hessian) of shape function i's
+       * c-th vector component, we need to look it up in the #shape_values,
+       * #shape_gradients and #shape_hessians arrays.  The question is where in
+       * this array does the data for shape function i, component c reside. This
+       * is what this table answers.
+       *
+       * The format of the table is as follows: - It has dofs_per_cell times
+       * n_components entries. - The entry that corresponds to shape function i,
+       * component c is <code>i * n_components + c</code>. - The value stored at
+       * this position indicates the row in #shape_values and the other tables
+       * where the corresponding datum is stored for all the quadrature points.
+       *
+       * In the general, vector-valued context, the number of components is larger
+       * than one, but for a given shape function, not all vector components may
+       * be nonzero (e.g., if a shape function is primitive, then exactly one
+       * vector component is non-zero, while the others are all zero). For such
+       * zero components, #shape_values and friends do not have a row.
+       * Consequently, for vector components for which shape function i is zero,
+       * the entry in the current table is numbers::invalid_unsigned_int.
+       *
+       * On the other hand, the table is guaranteed to have at least one valid
+       * index for each shape function. In particular, for a primitive finite
+       * element, each shape function has exactly one nonzero component and so for
+       * each i, there is exactly one valid index within the range
+       * <code>[i*n_components, (i+1)*n_components)</code>.
+       */
+      std::vector<unsigned int> shape_function_to_row_table;
+    };
+  }
+}
+
+
 /*@}*/
 
 
index dd8afbea03305a60e3bf069a3e95be3a32c124e7..f65a767a3a2ab486cde0328c6f736949f76b6100 100644 (file)
@@ -1243,38 +1243,19 @@ namespace internal
 
 
 //TODO: Add access to mapping values to FEValuesBase
-//TODO: Several FEValuesBase of a system should share Mapping
 
 /**
  * A class that contains all data vectors for FEValues, FEFaceValues, and
  * FESubfaceValues.
  *
- * This class has been extracted from FEValuesBase to encapsulate in one
- * place all of the data, independent of the functions thater later
- * access this data in the public interfaces of the FEValues and related
- * classes. Consequently, this base class is protected in FEValuesBase.
- *
- * The second reason is because in FEValuesBase::reinit, we first need to
- * call Mapping::fill_fe_values() to compute mapping related data, and later
- * call FiniteElement::fill_fe_values() to compute shape function related
- * data. In the first step, Mapping::fill_fe_values() gets a pointer to
- * its own internal data structure and a pointer to the FEValuesData base
- * object of FEValuesBase, and the mapping then places the computed data
- * into the data fields that pertain to the mapping below. In the second
- * step, the finite element receives a pointer to its own internal object,
- * and to the current object, and from both of these computes the shape
- * function related information and, again, places it into the current
- * FEValuesData object.
- *
  * More information can be found on the page on
  * @ref UpdateFlagsEssay.
  *
  * @ingroup feaccess
- * @author Guido Kanschat
- * @date 2000
  */
 template <int dim, int spacedim=dim>
-class FEValuesData
+class FEValuesData : public internal::FEValues::MappingRelatedData<dim,spacedim>,
+  public internal::FEValues::FiniteElementRelatedData<dim,spacedim>
 {
 public:
   /**
@@ -1284,153 +1265,6 @@ public:
                    const FiniteElement<dim,spacedim> &fe,
                    const UpdateFlags         flags);
 
-  /**
-   * @name Fields filled by the finite element
-   * @{
-   */
-
-  /**
-   * Storage type for shape values. Each row in the matrix denotes the values
-   * of a single shape function at the different points, columns are for a
-   * single point with the different shape functions.
-   *
-   * If a shape function has more than one non-zero component (in deal.II
-   * diction: it is non-primitive), then we allocate one row per non-zero
-   * component, and shift subsequent rows backward.  Lookup of the correct row
-   * for a shape function is thus simple in case the entire finite element is
-   * primitive (i.e. all shape functions are primitive), since then the shape
-   * function number equals the row number. Otherwise, use the
-   * #shape_function_to_row_table array to get at the first row that belongs
-   * to this particular shape function, and navigate among all the rows for
-   * this shape function using the FiniteElement::get_nonzero_components()
-   * function which tells us which components are non-zero and thus have a row
-   * in the array presently under discussion.
-   */
-  typedef Table<2,double> ShapeVector;
-
-  /**
-   * Storage type for gradients. The layout of data is the same as for the
-   * #ShapeVector data type.
-   */
-  typedef std::vector<std::vector<Tensor<1,spacedim> > > GradientVector;
-
-  /**
-   * Likewise for second order derivatives.
-   */
-  typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
-
-  /**
-   * Store the values of the shape functions at the quadrature points. See the
-   * description of the data type for the layout of the data in this field.
-   */
-  ShapeVector shape_values;
-
-  /**
-   * Store the gradients of the shape functions at the quadrature points. See
-   * the description of the data type for the layout of the data in this
-   * field.
-   */
-  GradientVector shape_gradients;
-
-  /**
-   * Store the 2nd derivatives of the shape functions at the quadrature
-   * points.  See the description of the data type for the layout of the data
-   * in this field.
-   */
-  HessianVector shape_hessians;
-
-  /**
-   * @}
-   */
-
-  /**
-   * @name Fields filled by the mapping
-   * @{
-   */
-
-  /**
-   * Store an array of weights times the Jacobi determinant at the quadrature
-   * points. This function is reset each time reinit() is called. The Jacobi
-   * determinant is actually the reciprocal value of the Jacobi matrices
-   * stored in this class, see the general documentation of this class for
-   * more information.
-   *
-   * However, if this object refers to an FEFaceValues or FESubfaceValues
-   * object, then the JxW_values correspond to the Jacobian of the
-   * transformation of the face, not the cell, i.e. the dimensionality is that
-   * of a surface measure, not of a volume measure. In this case, it is
-   * computed from the boundary forms, rather than the Jacobian matrix.
-   */
-  std::vector<double>       JxW_values;
-
-  /**
-   * Array of the Jacobian matrices at the quadrature points.
-   */
-  std::vector< DerivativeForm<1,dim,spacedim> > jacobians;
-
-  /**
-   * Array of the derivatives of the Jacobian matrices at the quadrature
-   * points.
-   */
-  std::vector<DerivativeForm<2,dim,spacedim> >  jacobian_grads;
-
-  /**
-   * Array of the inverse Jacobian matrices at the quadrature points.
-   */
-  std::vector<DerivativeForm<1,spacedim,dim> > inverse_jacobians;
-
-  /**
-   * Array of quadrature points. This array is set up upon calling reinit()
-   * and contains the quadrature points on the real element, rather than on
-   * the reference element.
-   */
-  std::vector<Point<spacedim> >  quadrature_points;
-
-  /**
-   * List of outward normal vectors at the quadrature points. This field is
-   * filled in by the finite element class.
-   */
-  std::vector<Point<spacedim> >  normal_vectors;
-
-  /**
-   * List of boundary forms at the quadrature points. This field is filled in
-   * by the finite element class.
-   */
-  std::vector<Tensor<1,spacedim> >  boundary_forms;
-
-  /**
-   * @}
-   */
-
-  /**
-   * When asked for the value (or gradient, or Hessian) of shape function i's
-   * c-th vector component, we need to look it up in the #shape_values,
-   * #shape_gradients and #shape_hessians arrays.  The question is where in
-   * this array does the data for shape function i, component c reside. This
-   * is what this table answers.
-   *
-   * The format of the table is as follows: - It has dofs_per_cell times
-   * n_components entries. - The entry that corresponds to shape function i,
-   * component c is <code>i * n_components + c</code>. - The value stored at
-   * this position indicates the row in #shape_values and the other tables
-   * where the corresponding datum is stored for all the quadrature points.
-   *
-   * In the general, vector-valued context, the number of components is larger
-   * than one, but for a given shape function, not all vector components may
-   * be nonzero (e.g., if a shape function is primitive, then exactly one
-   * vector component is non-zero, while the others are all zero). For such
-   * zero components, #shape_values and friends do not have a row.
-   * Consequently, for vector components for which shape function i is zero,
-   * the entry in the current table is numbers::invalid_unsigned_int.
-   *
-   * On the other hand, the table is guaranteed to have at least one valid
-   * index for each shape function. In particular, for a primitive finite
-   * element, each shape function has exactly one nonzero component and so for
-   * each i, there is exactly one valid index within the range
-   * <code>[i*n_components, (i+1)*n_components)</code>.
-   */
-  std::vector<unsigned int> shape_function_to_row_table;
-
   /**
    * Original update flags handed to the constructor of FEValues.
    */
index 1f812b41703b46d7795cd5fbe4f90659f1689ea3..12b187c14d588f321027bb19a4ac02859281e0b0 100644 (file)
@@ -96,7 +96,7 @@ namespace internal
        * mapped quadrature points are accessible, as no finite element data is
        * actually used).
        */
-      const FEValues<dim> &get_fe_values () const;
+      const dealii::FEValues<dim> &get_fe_values () const;
 
       /**
        * Return a vector of inverse transpose Jacobians. For compatibility
@@ -152,7 +152,7 @@ namespace internal
       /**
        * An underlying FEValues object that performs the (scalar) evaluation.
        */
-      FEValues<dim> fe_values;
+      dealii::FEValues<dim> fe_values;
 
       /**
        * Get 1D quadrature formula to be used for reinitializing shape info.
@@ -272,7 +272,7 @@ namespace internal
 
     template <int dim, typename Number>
     inline
-    const FEValues<dim> &
+    const dealii::FEValues<dim> &
     MappingDataOnTheFly<dim,Number>::get_fe_values() const
     {
       return fe_values;
index a5293eab9d08c9e347b3c3b97875a94968fe4468..906e2ae6ad8126d2e6f6428ca1816b020e052be9 100644 (file)
@@ -328,7 +328,7 @@ namespace internal
                              const unsigned int  my_q,
                              CellType (&cell_t_prev)[n_vector_elements],
                              CellType (&cell_t)[n_vector_elements],
-                             FEValues<dim,dim> &fe_values,
+                             dealii::FEValues<dim,dim> &fe_values,
                              CellData          &cell_data) const;
     };
 
index d2f74ad0861ce5c25b633637d557d2cfbb4a0dd4..1c68ab104e133f76d9f8bf0c0fd242a3c7027f48 100644 (file)
@@ -240,7 +240,7 @@ namespace internal
           // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
           // select the correct finite element, so just hold a vector of
           // FEValues
-          std::vector<std_cxx11::shared_ptr<FEValues<dim> > >
+          std::vector<std_cxx11::shared_ptr<dealii::FEValues<dim> > >
           fe_values (current_data.quadrature.size());
           UpdateFlags update_flags_feval =
             (update_flags & update_inverse_jacobians ? update_jacobians : update_default) |
@@ -290,10 +290,10 @@ namespace internal
               const unsigned int n_q_points = current_data.n_q_points[fe_index];
               if (fe_values[fe_index].get() == 0)
                 fe_values[fe_index].reset
-                (new FEValues<dim> (mapping, dummy_fe,
-                                    current_data.quadrature[fe_index],
-                                    update_flags_feval));
-              FEValues<dim> &fe_val = *fe_values[fe_index];
+                (new dealii::FEValues<dim> (mapping, dummy_fe,
+                                            current_data.quadrature[fe_index],
+                                            update_flags_feval));
+              dealii::FEValues<dim> &fe_val = *fe_values[fe_index];
               data.resize (n_q_points);
 
               // if the fe index has changed from the previous cell, set the
@@ -608,7 +608,7 @@ namespace internal
                                                const unsigned int  my_q,
                                                CellType (&cell_t_prev)[n_vector_elements],
                                                CellType (&cell_t)[n_vector_elements],
-                                               FEValues<dim,dim> &fe_val,
+                                               dealii::FEValues<dim,dim> &fe_val,
                                                CellData          &data) const
     {
       const unsigned int n_q_points = fe_val.n_quadrature_points;
index 909f6d11d3bd1e5a672e16ee23d6c679cd4bd7d0..49792db89b8110d78df22141663f63bf1f561e23 100644 (file)
@@ -2073,67 +2073,93 @@ get_interpolated_dof_values (const IndexSet &,
 /* --------------------- FEValuesData ----------------- */
 
 
+namespace internal
+{
+  namespace FEValues
+  {
+    template <int dim, int spacedim>
+    void
+    MappingRelatedData<dim,spacedim>::initialize (const unsigned int        n_quadrature_points,
+                                                  const UpdateFlags         flags)
+    {
+      if (flags & update_quadrature_points)
+        this->quadrature_points.resize(n_quadrature_points);
+
+      if (flags & update_JxW_values)
+        this->JxW_values.resize(n_quadrature_points);
+
+      if (flags & update_jacobians)
+        this->jacobians.resize(n_quadrature_points);
+
+      if (flags & update_jacobian_grads)
+        this->jacobian_grads.resize(n_quadrature_points);
+
+      if (flags & update_inverse_jacobians)
+        this->inverse_jacobians.resize(n_quadrature_points);
+
+      if (flags & update_boundary_forms)
+        this->boundary_forms.resize(n_quadrature_points);
+
+      if (flags & update_normal_vectors)
+        this->normal_vectors.resize(n_quadrature_points);
+    }
+
+
+    template <int dim, int spacedim>
+    void
+    FiniteElementRelatedData<dim,spacedim>::initialize (const unsigned int        n_quadrature_points,
+                                                        const FiniteElement<dim,spacedim> &fe,
+                                                        const UpdateFlags         flags)
+    {
+
+      // initialize the table mapping
+      // from shape function number to
+      // the rows in the tables storing
+      // the data by shape function and
+      // nonzero component
+      this->shape_function_to_row_table
+        = make_shape_function_to_row_table (fe);
+
+      // count the total number of non-zero
+      // components accumulated over all shape
+      // functions
+      unsigned int n_nonzero_shape_components = 0;
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+        n_nonzero_shape_components += fe.n_nonzero_components (i);
+      Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
+              ExcInternalError());
+
+      // with the number of rows now
+      // known, initialize those fields
+      // that we will need to their
+      // correct size
+      if (flags & update_values)
+        this->shape_values.reinit(n_nonzero_shape_components,
+                                  n_quadrature_points);
+
+      if (flags & update_gradients)
+        this->shape_gradients.resize (n_nonzero_shape_components,
+                                      std::vector<Tensor<1,spacedim> > (n_quadrature_points));
+
+      if (flags & update_hessians)
+        this->shape_hessians.resize (n_nonzero_shape_components,
+                                     std::vector<Tensor<2,spacedim> > (n_quadrature_points));
+    }
+  }
+}
+
+
 template <int dim, int spacedim>
 void
 FEValuesData<dim,spacedim>::initialize (const unsigned int        n_quadrature_points,
                                         const FiniteElement<dim,spacedim> &fe,
                                         const UpdateFlags         flags)
 {
-  this->update_flags = flags;
+  // initialize the base classes
+  internal::FEValues::MappingRelatedData<dim,spacedim>::initialize(n_quadrature_points, flags);
+  internal::FEValues::FiniteElementRelatedData<dim,spacedim>::initialize(n_quadrature_points, fe, flags);
 
-  // initialize the table mapping
-  // from shape function number to
-  // the rows in the tables storing
-  // the data by shape function and
-  // nonzero component
-  this->shape_function_to_row_table
-    = make_shape_function_to_row_table (fe);
-
-  // count the total number of non-zero
-  // components accumulated over all shape
-  // functions
-  unsigned int n_nonzero_shape_components = 0;
-  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    n_nonzero_shape_components += fe.n_nonzero_components (i);
-  Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
-          ExcInternalError());
-
-  // with the number of rows now
-  // known, initialize those fields
-  // that we will need to their
-  // correct size
-  if (flags & update_values)
-    this->shape_values.reinit(n_nonzero_shape_components,
-                              n_quadrature_points);
-
-  if (flags & update_gradients)
-    this->shape_gradients.resize (n_nonzero_shape_components,
-                                  std::vector<Tensor<1,spacedim> > (n_quadrature_points));
-
-  if (flags & update_hessians)
-    this->shape_hessians.resize (n_nonzero_shape_components,
-                                 std::vector<Tensor<2,spacedim> > (n_quadrature_points));
-
-  if (flags & update_quadrature_points)
-    this->quadrature_points.resize(n_quadrature_points);
-
-  if (flags & update_JxW_values)
-    this->JxW_values.resize(n_quadrature_points);
-
-  if (flags & update_jacobians)
-    this->jacobians.resize(n_quadrature_points);
-
-  if (flags & update_jacobian_grads)
-    this->jacobian_grads.resize(n_quadrature_points);
-
-  if (flags & update_inverse_jacobians)
-    this->inverse_jacobians.resize(n_quadrature_points);
-
-  if (flags & update_boundary_forms)
-    this->boundary_forms.resize(n_quadrature_points);
-
-  if (flags & update_normal_vectors)
-    this->normal_vectors.resize(n_quadrature_points);
+  this->update_flags = flags;
 }
 
 

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.