]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement a faster scheme to compute the FEValues by storing the values and gradients...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 1998 16:08:17 +0000 (16:08 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 1998 16:08:17 +0000 (16:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@447 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_values.cc

index cfe3502f39de81b64710464b92598a75ca85eee1..87c4b073e38e4c83ba37068cffeb51087db902f7 100644 (file)
@@ -53,28 +53,34 @@ struct FiniteElementData<1> {
                                      */
     const unsigned int total_dofs;
 
+                                    /**
+                                     * Number of basis functions used for the
+                                     * transformation from unit cell to real
+                                     * cell. For a linear mapping, this number
+                                     * equals the number of vertices.
+                                     */
+    const unsigned int n_transform_functions;
+    
                                     /**
                                      * Default constructor. Constructs an element
                                      * which is not so useful. Checking
                                      * #total_dofs# is therefore a good way to
                                      * check if something went wrong.
                                      */
-    FiniteElementData () :
-                   dofs_per_vertex(0),
-                   dofs_per_line(0),
-                   dofs_per_face(0),
-                   total_dofs(0) {};
+    FiniteElementData ();
 
                                     /**
                                      * A more useful version to construct
                                      * an object of this type.
                                      */
     FiniteElementData (const unsigned int dofs_per_vertex,
-                      const unsigned int dofs_per_line) :
+                      const unsigned int dofs_per_line,
+                      const unsigned int n_transform_functions) :
                    dofs_per_vertex(dofs_per_vertex),
                    dofs_per_line(dofs_per_line),
                    dofs_per_face(dofs_per_vertex),
-                   total_dofs (2*dofs_per_vertex+dofs_per_line) {};
+                   total_dofs (2*dofs_per_vertex+dofs_per_line),
+                   n_transform_functions(n_transform_functions) {};
 
                                     /**
                                      * Comparison operator. It is not clear to
@@ -82,6 +88,11 @@ struct FiniteElementData<1> {
                                      * this one explicitely.
                                      */
     bool operator == (const FiniteElementData<1> &) const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInternalError);
 };
 
 
@@ -128,18 +139,21 @@ struct FiniteElementData<2> {
                                      */
     const unsigned int total_dofs;
 
+                                    /**
+                                     * Number of basis functions used for the
+                                     * transformation from unit cell to real
+                                     * cell. For a linear mapping, this number
+                                     * equals the number of vertices.
+                                     */
+    const unsigned int n_transform_functions;
+
                                     /**
                                      * Default constructor. Constructs an element
                                      * which is not so useful. Checking
                                      * #total_dofs# is therefore a good way to
                                      * check if something went wrong.
                                      */
-    FiniteElementData () :
-                   dofs_per_vertex(0),
-                   dofs_per_line(0),
-                   dofs_per_quad(0),
-                   dofs_per_face(0),
-                   total_dofs(0) {};
+    FiniteElementData ();
 
                                     /**
                                      * A more useful version to construct
@@ -147,7 +161,8 @@ struct FiniteElementData<2> {
                                      */
     FiniteElementData (const unsigned int dofs_per_vertex,
                       const unsigned int dofs_per_line,
-                      const unsigned int dofs_per_quad) :
+                      const unsigned int dofs_per_quad,
+                      const unsigned int n_transform_functions) :
                    dofs_per_vertex(dofs_per_vertex),
                    dofs_per_line(dofs_per_line),
                    dofs_per_quad(dofs_per_quad),
@@ -155,7 +170,8 @@ struct FiniteElementData<2> {
                                  dofs_per_line),
                    total_dofs (4*dofs_per_vertex+
                                4*dofs_per_line+
-                               dofs_per_quad) {};
+                               dofs_per_quad),
+                   n_transform_functions (n_transform_functions) {};
 
                                     /**
                                      * Comparison operator. It is not clear to
@@ -163,6 +179,11 @@ struct FiniteElementData<2> {
                                      * this one explicitely.
                                      */
     bool operator == (const FiniteElementData<2> &) const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInternalError);
 };
 
     
@@ -202,7 +223,8 @@ struct FiniteElementBase : public FiniteElementData<dim> {
                                      */
     FiniteElementBase (const unsigned int dofs_per_vertex,
                       const unsigned int dofs_per_line,
-                      const unsigned int dofs_per_quad=0);
+                      const unsigned int dofs_per_quad,
+                      const unsigned int n_transform_functions);
     
                                     /**
                                      * Return a readonly reference to the
@@ -266,10 +288,6 @@ struct FiniteElementBase : public FiniteElementData<dim> {
                    << "The interface matrix has a size of " << arg1
                    << "x" << arg2
                    << ", which is not reasonable in the present dimension.");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
     
   protected:
                                     /**
@@ -474,10 +492,12 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      */
     FiniteElement (const unsigned int dofs_per_vertex,
                   const unsigned int dofs_per_line,
-                  const unsigned int dofs_per_quad=0) :
+                  const unsigned int dofs_per_quad,
+                  const unsigned int n_transform_functions) :
                    FiniteElementBase<dim> (dofs_per_vertex,
                                            dofs_per_line,
-                                           dofs_per_quad) {};
+                                           dofs_per_quad,
+                                           n_transform_functions) {};
 
                                     /**
                                      * Destructor. Only declared to have a
@@ -492,7 +512,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * #p# is a point on the reference element.
                                      */
     virtual double shape_value (const unsigned int i,
-                               const Point<dim>p) const = 0;
+                               const Point<dim> &p) const = 0;
 
                                     /**
                                      * Return the gradient of the #i#th shape
@@ -500,8 +520,28 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * #p# is a point on the reference element,
                                      */
     virtual Point<dim> shape_grad (const unsigned int i,
-                                  const Point<dim>p) const = 0;
+                                  const Point<dim> &p) const = 0;
 
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function of the transformation mapping
+                                     * from unit cell to real cell. For
+                                     * isoparametric elements, this function
+                                     * is the same as the ansatz functions,
+                                     * but for sublinear or other mappings,
+                                     * they differ.
+                                     */
+    virtual double shape_value_transform (const unsigned int i,
+                                         const Point<dim> &p) const = 0;
+
+                                    /**
+                                     * Same as above: return gradient of the
+                                     * #i#th shape function for the mapping
+                                     * from unit to real cell.
+                                     */
+    virtual Point<dim> shape_grad_transform (const unsigned int i,
+                                            const Point<dim> &p) const = 0;    
+    
                                     /**
                                      * Compute the Jacobian matrix and the
                                      * quadrature points as well as the ansatz
@@ -538,6 +578,17 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * unit cell to the real cell, which
                                      * makes things much more complicated.
                                      *
+                                     * The #shape_values/grads_transform#
+                                     * arrays store the values and gradients
+                                     * of the transformation basis functions.
+                                     * While this information is not necessary
+                                     * for the computation of the other fields,
+                                     * it allows for significant speedups, since
+                                     * the values and gradients of the transform
+                                     * functions at the quadrature points
+                                     * need not be recomputed each time this
+                                     * function is called.
+                                     *
                                      * The function assumes that the fields
                                      * already have the right number of
                                      * elements.
@@ -555,6 +606,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
+                                const dFMatrix      &shape_values_transform,
+                                const vector<vector<Point<dim> > > &shape_grads_transform,
                                 const Boundary<dim> &boundary) const;
 
                                     /**
@@ -675,6 +728,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      const bool           compute_face_jacobians,
                                      vector<Point<dim> > &normal_vectors,
                                      const bool           compute_normal_vectors,
+                                     const dFMatrix      &shape_values_transform,
+                                     const vector<vector<Point<dim> > > &shape_grads_transform,
                                      const Boundary<dim> &boundary) const;
 
                                     /**
@@ -718,6 +773,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                         const bool           compute_face_jacobians,
                                         vector<Point<dim> > &normal_vectors,
                                         const bool           compute_normal_vectors,
+                                        const dFMatrix      &shape_values_transform,
+                                        const vector<vector<Point<dim> > > &shape_grads_transform,
                                         const Boundary<dim> &boundary) const;
 
                                     /**
@@ -984,7 +1041,30 @@ class FELinearMapping : public FiniteElement<dim> {
                     const unsigned int dofs_per_quad=0) :
                    FiniteElement<dim> (dofs_per_vertex,
                                        dofs_per_line,
-                                       dofs_per_quad) {};
+                                       dofs_per_quad,
+                                       GeometryInfo<dim>::vertices_per_cell) {};
+
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at point #p# on the unit cell.
+                                     * Here, the (bi-)linear basis functions
+                                     * are meant, which are used for the
+                                     * computation of the transformation from
+                                     * unit cell to real space cell.
+                                     */
+    virtual double shape_value_transform (const unsigned int i,
+                                         const Point<dim> &p) const;
+
+                                    /**
+                                     * Return the gradient of the #i#th shape
+                                     * function at point #p# on the unit cell.
+                                     * Here, the (bi-)linear basis functions
+                                     * are meant, which are used for the
+                                     * computation of the transformation from
+                                     * unit cell to real space cell.
+                                     */
+    virtual Point<dim> shape_grad_transform (const unsigned int i,
+                                            const Point<dim> &p) const;
 
                                     /**
                                      * Refer to the base class for detailed
@@ -1073,30 +1153,9 @@ class FELinearMapping : public FiniteElement<dim> {
                                 const bool           compute_ansatz_points,
                                 vector<Point<dim> > &q_points,
                                 const bool           compute_q_points,
+                                const dFMatrix      &shape_values_transform,
+                                const vector<vector<Point<dim> > > &shape_grad_transform,
                                 const Boundary<dim> &boundary) const;
-
-  protected:
-                                    /**
-                                     * Return the value of the #i#th shape
-                                     * function at point #p# on the unit cell.
-                                     * Here, the (bi-)linear basis functions
-                                     * are meant, which are used for the
-                                     * computation of the transformation from
-                                     * unit cell to real space cell.
-                                     */
-    double linear_shape_value(const unsigned int i,
-                             const Point<dim>& p) const;
-
-                                    /**
-                                     * Return the gradient of the #i#th shape
-                                     * function at point #p# on the unit cell.
-                                     * Here, the (bi-)linear basis functions
-                                     * are meant, which are used for the
-                                     * computation of the transformation from
-                                     * unit cell to real space cell.
-                                     */
-    Point<dim> linear_shape_grad(const unsigned int i,
-                                const Point<dim>& p) const;
 };
 
 
index 7b153329a32f3a52b9ea97eab5f3d2748fe76318..fc156f0f1fc175bb41a41a9e5525e128ad5cab40 100644 (file)
@@ -57,6 +57,32 @@ template <int dim> class Quadrature;
  *  the derived classes and the #get_values# function for the exact usage of
  *  this variable.
  *
+ *  For many of the actual computations done by the #fill_fe_*# functions of
+ *  the #FiniteElement# class and its decendants, the values and gradients of
+ *  the transformation functions are needed. For example, for the computation
+ *  of the real location of a quadrature point from the location on the unit
+ *  cell, the values are needed, while for the computation of the Jacobian
+ *  matrix the gradient is needed. While for linear elements the transformation
+ *  functions coincide with the ansatz functions, this does not hold for higher
+ *  order elements with subparametric mappings and for other types of elements
+ *  such as non-conforming ones, etc, such that the precomputed values and
+ *  gradients of the ansatz functions (#unit_shape_values# and
+ *  #unit_shape_grads#) cannot be used for the present purpose.
+ *  In principle, these values could be computed each time the #fill_fe_*#
+ *  function is called; however, this computation is highly redundant, since
+ *  only the values on the unit cell and only at the quadrature points are
+ *  needed, i.e. they are the same for each cell that #fill_fe_*# is called.
+ *  Therefore, two additional arrays, #unit_shape_values_transform# and
+ *  #unit_shape_grads_transform# are provided, which are filled upon construction
+ *  of an object of this type, which the actual finite element may or may not
+ *  use. Later on, the #fill_fe_*# functions are passed pointers to these
+ *  arrays, which they may use to extract the values and gradients of the
+ *  transform functions. If a concrete finite element choses not to use this
+ *  field, it shall set its field #n_transform_functions# to zero.
+ *
+ *  The #unit_shape_grads_transform# array is provided by the derived classes
+ *  to allow for the inclusion of multiple faces, etc.
+ *
  *
  *  \subsection{Definitions}
  *
@@ -141,7 +167,6 @@ template <int dim> class Quadrature;
  * not assume that someone who wants to get the #JxW_values# must set the
  * #update_jacobians# flag besides the #update_JxW_values# flag.
  *
- *
  * It is also forbidden that the constructor of a class set the
  * #update_jacobians# flag if the user specifies #update_JxW_values#. This is
  * since derived classes may be able to compute the #JxW_values# field without
@@ -169,6 +194,14 @@ class FEValuesBase {
                                      */
     const unsigned int total_dofs;
 
+                                    /**
+                                     * Number of basis functions for the
+                                     * transformation from the unit cell
+                                     * to the real cell. See the docs for
+                                     * more information on this field.
+                                     */
+    const unsigned int n_transform_functions;
+    
                                     /**
                                      * Constructor. Set up the array sizes
                                      * with #n_q_points# quadrature points,
@@ -186,6 +219,7 @@ class FEValuesBase {
     FEValuesBase (const unsigned int n_q_points,
                  const unsigned int n_ansatz_points,
                  const unsigned int n_dofs,
+                 const unsigned int n_transform_functions,
                  const unsigned int n_values_array,
                  const UpdateFlags update_flags);
     
@@ -470,6 +504,27 @@ class FEValuesBase {
                                      */
     vector<dFMatrix>     jacobi_matrices;
 
+                                    /**
+                                     * Store the values of the basis functions
+                                     * of the transformation from unit cell
+                                     * to real cell at the quadrature points.
+                                     *
+                                     * This field stores some data which is not
+                                     * really needed for the assemblage of
+                                     * matrices and vectors but makes that
+                                     * operation much faster. Each time the
+                                     * #FEValues::reinit# function calls
+                                     * the #FiniteElemenet::fill_fe_values#
+                                     * function, this and the next array are
+                                     * passed. The #fill_fe_values# function
+                                     * may or may not use this field.
+                                     *
+                                     * The element #(i,j)# denotes the value
+                                     * of the #i#th transfer basis function
+                                     * at the #j#th quadrature point.
+                                     */
+    vector<dFMatrix>     shape_values_transform;
+
                                     /**
                                      * Store which of the data sets in the
                                      * #shape_values# array is presently
@@ -484,7 +539,7 @@ class FEValuesBase {
                                      * Store which fields are to be updated by
                                      * the reinit function.
                                      */
-    UpdateFlags         update_flags;
+    UpdateFlags          update_flags;
 
                                     /**
                                      * Store the cell selected last time
@@ -575,7 +630,14 @@ class FEValues : public FEValuesBase<dim> {
                                      */
     vector<vector<Point<dim> > >   unit_shape_gradients;
     
-
+                                    /**
+                                     * Gradients of the basis
+                                     * functions of the transformation.
+                                     * Analogous to the #shape_values_transform#
+                                     * array of the base class.
+                                     */
+    vector<vector<Point<dim> > >   unit_shape_gradients_transform;
+    
                                     /**
                                      * Array of quadrature points in the unit
                                      * cell. This array is set up upon
@@ -680,6 +742,7 @@ class FEFaceValuesBase : public FEValuesBase<dim> {
     FEFaceValuesBase (const unsigned int n_q_points,
                      const unsigned int n_ansatz_points,
                      const unsigned int n_dofs,
+                     const unsigned int n_transform_functions,
                      const unsigned int n_faces_or_subfaces,
                      const UpdateFlags update_flags);
 
@@ -712,6 +775,14 @@ class FEFaceValuesBase : public FEValuesBase<dim> {
                                      * #unit_shape_gradients[face][dof][q_point]#
                                      */
     vector<vector<vector<Point<dim> > > > unit_shape_gradients;
+    
+                                    /**
+                                     * Gradients of the basis
+                                     * functions of the transformation.
+                                     * Analogous to the #shape_values_transform#
+                                     * array of the base class.
+                                     */
+    vector<vector<vector<Point<dim> > > > unit_shape_gradients_transform;
 
                                     /**
                                      * Array of quadrature points on the
index fefe9e5745a96f3315b65d75f009fccf7eaac504..cc8e4a672edd8e08c003e3958a72362a1a6e7cfb 100644 (file)
 
 #if deal_II_dimension == 1
 
+FiniteElementData<1>::FiniteElementData () :
+               dofs_per_vertex(0),
+               dofs_per_line(0),
+               dofs_per_face(0),
+               total_dofs(0),
+               n_transform_functions(0) {
+  Assert (false, ExcInternalError());
+};
+
+
+
 bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const {
   return ((dofs_per_vertex == f.dofs_per_vertex) &&
          (dofs_per_line == f.dofs_per_line) &&
@@ -26,6 +37,19 @@ bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const {
 
 #if deal_II_dimension == 2
 
+
+FiniteElementData<2>::FiniteElementData () :
+               dofs_per_vertex(0),
+               dofs_per_line(0),
+               dofs_per_quad(0),
+               dofs_per_face(0),
+               total_dofs(0),
+               n_transform_functions(0) {
+  Assert (false, ExcInternalError());
+};
+
+
+
 bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const {
   return ((dofs_per_vertex == f.dofs_per_vertex) &&
          (dofs_per_line == f.dofs_per_line) &&
@@ -45,9 +69,11 @@ bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const {
 template <>
 FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex,
                                         const unsigned int dofs_per_line,
-                                        const unsigned int dofs_per_quad) :
+                                        const unsigned int dofs_per_quad,
+                                        const unsigned int n_transform_funcs) :
                FiniteElementData<1> (dofs_per_vertex,
-                                     dofs_per_line)
+                                     dofs_per_line,
+                                     n_transform_functs)
 {
   Assert (dofs_per_quad==0, ExcInternalError());
 
@@ -69,10 +95,12 @@ FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex,
 template <>
 FiniteElementBase<2>::FiniteElementBase (const unsigned int dofs_per_vertex,
                                         const unsigned int dofs_per_line,
-                                        const unsigned int dofs_per_quad) :
+                                        const unsigned int dofs_per_quad,
+                                        const unsigned int n_transform_funcs) :
                FiniteElementData<2> (dofs_per_vertex,
                                      dofs_per_line,
-                                     dofs_per_quad)
+                                     dofs_per_quad,
+                                     n_transform_funcs)
 {
   const unsigned int dim=2;
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i) 
@@ -153,6 +181,8 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
                                       const bool         compute_ansatz_points,
                                       vector<Point<1> > &q_points,
                                       const bool         compute_q_points,
+                                      const dFMatrix      &shape_values_transform,
+                                      const vector<vector<Point<dim> > > &shape_grad_transform,
                                       const Boundary<1> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -198,6 +228,8 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &
                                            const bool              ,
                                            vector<Point<1> >       &,
                                            const bool,
+                                           const dFMatrix          &,
+                                           const vector<vector<Point<1> > > &,
                                            const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
 };
@@ -217,6 +249,8 @@ void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterato
                                               const bool               ,
                                               vector<Point<1> >       &,
                                               const bool,
+                                              const dFMatrix          &,
+                                              const vector<vector<Point<1> > > &,
                                               const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
 };
@@ -260,6 +294,8 @@ void FiniteElement<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
                                         const bool,
                                         vector<Point<dim> > &,
                                         const bool,
+                                        const dFMatrix      &,
+                                        const vector<vector<Point<dim> > > &,
                                         const Boundary<dim> &) const {
   Assert (false, ExcPureFunctionCalled());
 };
@@ -281,6 +317,8 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                                              const bool           compute_face_jacobians,
                                              vector<Point<dim> > &normal_vectors,
                                              const bool           compute_normal_vectors,
+                                             const dFMatrix      &shape_values_transform,
+                                             const vector<vector<Point<dim> > > &shape_gradients_transform,
                                              const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -296,6 +334,7 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                  jacobians, compute_jacobians,
                  dummy, false,
                  q_points, compute_q_points,
+                 shape_values_transform, shape_gradients_transform,
                  boundary);
   
   if (compute_ansatz_points)
@@ -327,6 +366,8 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
                                                 const bool           compute_face_jacobians,
                                                 vector<Point<dim> > &normal_vectors,
                                                 const bool           compute_normal_vectors,
+                                                const dFMatrix      &shape_values_transform,
+                                                const vector<vector<Point<dim> > > &shape_gradients_transform,
                                                 const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -340,6 +381,7 @@ void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_ite
                  jacobians, compute_jacobians,
                  dummy, false,
                  q_points, compute_q_points,
+                 shape_values_transform, shape_gradients_transform,
                  boundary);
   
   if (compute_face_jacobians)
@@ -362,13 +404,16 @@ void FiniteElement<dim>::get_ansatz_points (const DoFHandler<dim>::cell_iterator
 
 
 
+
+/*---------------------------- FELinearMapping ----------------------------------*/
+
 #if deal_II_dimension == 1
 
 template <>
 inline
 double
-FELinearMapping<1>::linear_shape_value(const unsigned int i,
-                                      const Point<1>     &p) const
+FELinearMapping<1>::shape_value_transform (const unsigned int i,
+                                          const Point<1>     &p) const
 {
   Assert((i<2), ExcInvalidIndex(i));
   const double xi = p(0);
@@ -385,8 +430,8 @@ FELinearMapping<1>::linear_shape_value(const unsigned int i,
 template <>
 inline
 Point<1>
-FELinearMapping<1>::linear_shape_grad(const unsigned int i,
-                                    const Point<1>&) const
+FELinearMapping<1>::shape_grad_transform(const unsigned int i,
+                                        const Point<1>&) const
 {
   Assert((i<2), ExcInvalidIndex(i));
   switch (i)
@@ -400,7 +445,7 @@ FELinearMapping<1>::linear_shape_grad(const unsigned int i,
 
 
 template <>
-void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+void FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
                                            const Boundary<1>         &,
                                            const vector<Point<0> > &,
                                            vector<double>      &) const {
@@ -410,7 +455,7 @@ void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &
 
 
 template <>
-void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
                                               const unsigned int           ,
                                               const vector<Point<0> > &,
                                               vector<double>      &) const {
@@ -420,7 +465,7 @@ void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterato
 
 
 template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
                                            const unsigned int,
                                            const Boundary<1> &,
                                            const vector<Point<0> > &,
@@ -431,7 +476,7 @@ void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &
 
 
 template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
                                            const unsigned int,
                                            const unsigned int,
                                            const vector<Point<0> > &,
@@ -449,12 +494,16 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel
                                         const bool         compute_ansatz_points,
                                         vector<Point<1> > &q_points,
                                         const bool         compute_q_points,
+                                        const dFMatrix      &shape_values_transform,
+                                        const vector<vector<Point<dim> > > &shape_gradients_transform,
                                         const Boundary<1> &boundary) const {
                                   // simply pass down
   FiniteElement<1>::fill_fe_values (cell, unit_points,
                                    jacobians, compute_jacobians,
                                    ansatz_points, compute_ansatz_points,
-                                   q_points, compute_q_points, boundary);
+                                   q_points, compute_q_points,
+                                   shape_values_transform, shape_gradients_transform,
+                                   boundary);
 };
 
 
@@ -468,8 +517,8 @@ void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cel
 template <>
 inline
 double
-FELinearMapping<2>::linear_shape_value (const unsigned int i,
-                                      const Point<2>& p) const
+FELinearMapping<2>::shape_value_transform (const unsigned int i,
+                                          const Point<2>& p) const
 {
   Assert((i<4), ExcInvalidIndex(i));
   switch (i)
@@ -487,8 +536,8 @@ FELinearMapping<2>::linear_shape_value (const unsigned int i,
 template <>
 inline
 Point<2>
-FELinearMapping<2>::linear_shape_grad (const unsigned int i,
-                                     const Point<2>& p) const
+FELinearMapping<2>::shape_grad_transform (const unsigned int i,
+                                         const Point<2>& p) const
 {
   Assert((i<4), ExcInvalidIndex(i));
   switch (i)
@@ -629,6 +678,8 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                           const bool           compute_ansatz_points,
                                           vector<Point<dim> > &q_points,
                                           const bool           compute_q_points,
+                                          const dFMatrix      &shape_values_transform,
+                                          const vector<vector<Point<dim> > > &shape_grad_transform,
                                           const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -665,7 +716,7 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                       // use a subparametric mapping
       for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j) 
        for (unsigned int l=0; l<n_points; ++l) 
-         q_points[l] += vertices[j] * linear_shape_value(j, unit_points[l]);
+         q_points[l] += vertices[j] * shape_values_transform(j, l);
     };
   
 
@@ -697,7 +748,7 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
            {
                                               // we want the linear transform,
                                               // so use that function
-             const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
+             const Point<dim> gradient = shape_grad_transform[s][l];
              for (unsigned int i=0; i<dim; ++i)
                for (unsigned int j=0; j<dim; ++j)
                  M(i,j) += vertices[s](i) * gradient(j);
@@ -714,7 +765,6 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
 
 /*------------------------------- Explicit Instantiations -------------*/
 
-template class FiniteElementData<deal_II_dimension>;
 template class FiniteElementBase<deal_II_dimension>;
 template class FiniteElement<deal_II_dimension>;
 template class FELinearMapping<deal_II_dimension>;
index afa2c707910cc6e852499979591bbb6fb491a6b3..0a8c85e5dedf5e8d54ccdb8a0b441aa16837db0b 100644 (file)
@@ -17,10 +17,12 @@ template <int dim>
 FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
                                 const unsigned int n_ansatz_points,
                                 const unsigned int n_dofs,
+                                const unsigned int n_transform_functions,
                                 const unsigned int n_values_arrays,
                                 const UpdateFlags update_flags) :
                n_quadrature_points (n_q_points),
                total_dofs (n_dofs),
+               n_transform_functions (n_transform_functions),
                shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
                shape_gradients (n_dofs, vector<Point<dim> >(n_q_points)),
                weights (n_q_points, 0),
@@ -28,6 +30,9 @@ FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
                quadrature_points (n_q_points, Point<dim>()),
                ansatz_points (n_ansatz_points, Point<dim>()),
                jacobi_matrices (n_q_points, dFMatrix(dim,dim)),
+               shape_values_transform (n_values_arrays,
+                                       dFMatrix(n_transform_functions,
+                                                n_quadrature_points)),
                selected_dataset (0),
                update_flags (update_flags) {};
 
@@ -154,10 +159,13 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                FEValuesBase<dim> (quadrature.n_quadrature_points,
                                   fe.total_dofs,
                                   fe.total_dofs,
+                                  fe.n_transform_functions,
                                   1,
                                   update_flags),
                unit_shape_gradients(fe.total_dofs,
                                     vector<Point<dim> >(quadrature.n_quadrature_points)),
+               unit_shape_gradients_transform(fe.n_transform_functions,
+                                              vector<Point<dim> >(quadrature.n_quadrature_points)),
                unit_quadrature_points(quadrature.get_quad_points())
 {
   Assert ((update_flags & update_normal_vectors) == false,
@@ -171,6 +179,15 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
          = fe.shape_grad(i, unit_quadrature_points[j]);
       };
 
+  for (unsigned int i=0; i<n_transform_functions; ++i)
+    for (unsigned int j=0; j<n_quadrature_points; ++j)
+      {
+       shape_values_transform[0] (i,j)
+         = fe.shape_value_transform (i, unit_quadrature_points[j]);
+       unit_shape_gradients_transform[i][j]
+         = fe.shape_grad_transform(i, unit_quadrature_points[j]);
+      };
+  
   weights = quadrature.get_weights ();
 };
 
@@ -202,6 +219,7 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                       update_flags & update_ansatz_points,
                       quadrature_points,
                       update_flags & update_q_points,
+                      shape_values_transform[0], unit_shape_gradients_transform,
                       boundary);
 
                                   // compute gradients on real element if
@@ -244,16 +262,21 @@ template <int dim>
 FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
                                         const unsigned int n_ansatz_points,
                                         const unsigned int n_dofs,
+                                        const unsigned int n_transform_functions,
                                         const unsigned int n_faces_or_subfaces,
                                         const UpdateFlags update_flags) :
                FEValuesBase<dim> (n_q_points,
                                   n_ansatz_points,
                                   n_dofs,
+                                  n_transform_functions,
                                   n_faces_or_subfaces,
                                   update_flags),
                unit_shape_gradients (n_faces_or_subfaces,
                                      vector<vector<Point<dim> > >(n_dofs,
                                                                   vector<Point<dim> >(n_q_points))),
+               unit_shape_gradients_transform (n_faces_or_subfaces,
+                                               vector<vector<Point<dim> > >(n_transform_functions,
+                                                                            vector<Point<dim> >(n_q_points))),
                unit_face_quadrature_points (n_q_points, Point<dim-1>()),
                unit_quadrature_points (n_faces_or_subfaces,
                                        vector<Point<dim> >(n_q_points, Point<dim>())),
@@ -287,7 +310,8 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
                                       fe.dofs_per_face,
                                       fe.total_dofs,
-                                      2*dim,
+                                      fe.n_transform_functions,
+                                      GeometryInfo<dim>::faces_per_cell,
                                       update_flags)
 {
   unit_face_quadrature_points = quadrature.get_quad_points();
@@ -298,10 +322,10 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                                   // of the space with #dim# dimensions.
                                   // the points are still on the unit
                                   // cell, not on the real cell.
-  for (unsigned int face=0; face<2*dim; ++face)
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
     QProjector<dim>::project_to_face (quadrature, face, unit_quadrature_points[face]);
 
-  for (unsigned int face=0; face<2*dim; ++face)
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
     for (unsigned int i=0; i<fe.total_dofs; ++i)
       for (unsigned int j=0; j<n_quadrature_points; ++j) 
        {
@@ -310,6 +334,16 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
          unit_shape_gradients[face][i][j]
            = fe.shape_grad(i, unit_quadrature_points[face][j]);
        };
+
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int i=0; i<n_transform_functions; ++i)
+      for (unsigned int j=0; j<n_quadrature_points; ++j)
+       {
+         shape_values_transform[face] (i,j)
+           = fe.shape_value_transform (i, unit_quadrature_points[face][j]);
+         unit_shape_gradients_transform[face][i][j]
+           = fe.shape_grad_transform(i, unit_quadrature_points[face][j]);
+      };
 };
 
 
@@ -345,6 +379,8 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                            update_flags & update_JxW_values,
                            normal_vectors,
                            update_flags & update_normal_vectors,
+                           shape_values_transform[face_no],
+                           unit_shape_gradients_transform[face_no],
                            boundary);
 
                                   // compute gradients on real element if
@@ -392,7 +428,8 @@ FESubfaceValues<dim>::FESubfaceValues (const FiniteElement<dim> &fe,
                FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
                                       0,
                                       fe.total_dofs,
-                                      2*dim*(1<<(dim-1)),
+                                      fe.n_transform_functions,
+                                      GeometryInfo<dim>::faces_per_cell * GeometryInfo<dim>::subfaces_per_face,
                                       update_flags)
 {
   Assert ((update_flags & update_ansatz_points) == false,
@@ -406,21 +443,31 @@ FESubfaceValues<dim>::FESubfaceValues (const FiniteElement<dim> &fe,
                                   // of the space with #dim# dimensions.
                                   // the points are still on the unit
                                   // cell, not on the real cell.
-  for (unsigned int face=0; face<2*dim; ++face)
-    for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface)
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
       QProjector<dim>::project_to_subface (quadrature,
                                           face, subface,
                                           unit_quadrature_points[face*(1<<(dim-1))+subface]);
 
-  for (unsigned int face=0; face<2*dim; ++face)
-    for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface)
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
       for (unsigned int i=0; i<fe.total_dofs; ++i)
        for (unsigned int j=0; j<n_quadrature_points; ++j) 
          {
-           shape_values[face*(1<<(dim-1))+subface](i,j)
-             = fe.shape_value(i, unit_quadrature_points[face*(1<<(dim-1))+subface][j]);
-           unit_shape_gradients[face*(1<<(dim-1))+subface][i][j]
-             = fe.shape_grad(i, unit_quadrature_points[face*(1<<(dim-1))+subface][j]);
+           shape_values[face*GeometryInfo<dim>::subfaces_per_face+subface](i,j)
+             = fe.shape_value(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+           unit_shape_gradients[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+             = fe.shape_grad(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+         };
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
+      for (unsigned int i=0; i<n_transform_functions; ++i)
+       for (unsigned int j=0; j<n_quadrature_points; ++j)
+         {
+           shape_values_transform[face*GeometryInfo<dim>::subfaces_per_face+subface] (i,j)
+             = fe.shape_value_transform (i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+           unit_shape_gradients_transform[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+             = fe.shape_grad_transform(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
          };
 };
 
@@ -459,6 +506,8 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
                               update_flags & update_JxW_values,
                               normal_vectors,
                               update_flags & update_normal_vectors,
+                              shape_values_transform[selected_dataset],
+                              unit_shape_gradients_transform[selected_dataset],
                               boundary);
 
                                   // compute gradients on real element if

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.