]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Separate the linear mapping from the elements implemented as of yet.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 13 Jul 1998 14:01:58 +0000 (14:01 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 13 Jul 1998 14:01:58 +0000 (14:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@443 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.cubic.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/fe/fe_lib.quadratic.cc
deal.II/deal.II/source/fe/fe_lib.quartic.cc

index e04a12f35d8aca35852129e308a4ea4c6f98766f..cfe3502f39de81b64710464b92598a75ca85eee1 100644 (file)
@@ -966,6 +966,141 @@ class FiniteElement : public FiniteElementBase<dim> {
 
 
 
+/**
+ * Abstract base class for concrete finite elements which use a
+ * (bi-,tri-)linear mapping from the unit cell to the real cell. Some
+ * functions can be singled out from these elements and are collected
+ * in this one.
+ */
+template <int dim>
+class FELinearMapping : public FiniteElement<dim> {
+  public:
+                                    /**
+                                     * Constructor. Simply passes through
+                                     * its arguments to the base class.
+                                     */
+    FELinearMapping (const unsigned int dofs_per_vertex,
+                    const unsigned int dofs_per_line,
+                    const unsigned int dofs_per_quad=0) :
+                   FiniteElement<dim> (dofs_per_vertex,
+                                       dofs_per_line,
+                                       dofs_per_quad) {};
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     *
+                                     * In two spatial dimensions, this function
+                                     * simply returns the length of the face.
+                                     */
+    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<double>      &face_jacobi_determinants) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     *
+                                     * In two spatial dimensions, this function
+                                     * simply returns half the length of the
+                                     * whole face.
+                                     */
+    virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                       const unsigned int           subface_no,
+                                       const vector<Point<dim-1> > &unit_points,
+                                       vector<double>      &face_jacobi_determinants) const;
+
+                                    /**
+                                     * Return the normal vectors to the
+                                     * face with number #face_no# of #cell#.
+                                     *
+                                     * For linear finite elements, this function
+                                     * is particularly simple since all normal
+                                     * vectors are equal and can easiliy be
+                                     * computed from the direction of the face
+                                     * without using the transformation (Jacobi)
+                                     * matrix, at least for two dimensions.
+                                     *
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                    const unsigned int          face_no,
+                                    const Boundary<dim>         &boundary,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<Point<dim> >         &normal_vectors) const;
+
+                                    /**
+                                     * Return the normal vectors to the
+                                     * subface with number #subface_no# of
+                                     * the face with number #face_no# of #cell#.
+                                     *
+                                     * For linear finite elements, this function
+                                     * is particularly simple since all normal
+                                     * vectors are equal and can easiliy be
+                                     * computed from the direction of the face
+                                     * without using the transformation (Jacobi)
+                                     * matrix, at least for two dimensions.
+                                     *
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                    const unsigned int           face_no,
+                                    const unsigned int           subface_no,
+                                    const vector<Point<dim-1> > &unit_points,
+                                    vector<Point<dim> >         &normal_vectors) const;    
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     *
+                                     * For one dimensional elements, this
+                                     * function simply passes through to
+                                     * the one implemented in the base class.
+                                     * For higher dimensional finite elements
+                                     * we use linear mappings and therefore
+                                     * the boundary object is ignored since
+                                     * the boundary is approximated using
+                                     * piecewise multilinear boundary segments.
+                                     */
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
+                                vector<dFMatrix>    &jacobians,
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &ansatz_points,
+                                const bool           compute_ansatz_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points,
+                                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;
+};
+
+
+
   
 /*----------------------------   fe.h     ---------------------------*/
 /* end of #ifndef __fe_H */
index d11198b085050528e3373dbe7da6efd4378df9e3..45777de5c8deecda254fec71aef55fddb8326d9a 100644 (file)
@@ -25,7 +25,7 @@
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim>
-class FELinear : public FiniteElement<dim> {
+class FELinear : public FELinearMapping<dim> {
   public:
                                     /**
                                      * Constructor
@@ -46,29 +46,6 @@ class FELinear : public FiniteElement<dim> {
     virtual Point<dim> shape_grad(const unsigned int i,
                                  const Point<dim>& p) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * For one dimensional elements, this
-                                     * function simply passes through to
-                                     * the one implemented in the base class.
-                                     * For higher dimensional finite elements
-                                     * we use linear mappings and therefore
-                                     * the boundary object is ignored since
-                                     * the boundary is approximated using
-                                     * piecewise straight boundary segments.
-                                     */
-    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &ansatz_points,
-                                const bool           compute_ansatz_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const Boundary<dim> &boundary) const;
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -85,72 +62,6 @@ class FELinear : public FiniteElement<dim> {
                                         const Boundary<dim> &boundary,
                                         vector<Point<dim> > &ansatz_points) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * In two spatial dimensions, this function
-                                     * simply returns the length of the face.
-                                     */
-    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * In two spatial dimensions, this function
-                                     * simply returns half the length of the
-                                     * whole face.
-                                     */
-    virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                       const unsigned int           subface_no,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Return the normal vectors to the
-                                     * face with number #face_no# of #cell#.
-                                     *
-                                     * For linear finite elements, this function
-                                     * is particularly simple since all normal
-                                     * vectors are equal and can easiliy be
-                                     * computed from the direction of the face
-                                     * without using the transformation (Jacobi)
-                                     * matrix, at least for two dimensions.
-                                     *
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int          face_no,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
-                                    /**
-                                     * Return the normal vectors to the
-                                     * subface with number #subface_no# of
-                                     * the face with number #face_no# of #cell#.
-                                     *
-                                     * For linear finite elements, this function
-                                     * is particularly simple since all normal
-                                     * vectors are equal and can easiliy be
-                                     * computed from the direction of the face
-                                     * without using the transformation (Jacobi)
-                                     * matrix, at least for two dimensions.
-                                     *
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int           face_no,
-                                    const unsigned int           subface_no,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -169,7 +80,7 @@ class FELinear : public FiniteElement<dim> {
  * to the real cell is implemented.
  */
 template <int dim>
-class FEQuadraticSub : public FiniteElement<dim> {
+class FEQuadraticSub : public FELinearMapping<dim> {
   public:
                                     /**
                                      * Constructor
@@ -190,24 +101,6 @@ class FEQuadraticSub : public FiniteElement<dim> {
     virtual Point<dim> shape_grad(const unsigned int i,
                                  const Point<dim>& p) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * For one dimensional elements, this
-                                     * function simply passes through to
-                                     * the one implemented in the base class.
-                                     */
-    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &ansatz_points,
-                                const bool           compute_ansatz_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const Boundary<dim> &boundary) const;
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -224,44 +117,6 @@ class FEQuadraticSub : public FiniteElement<dim> {
                                         const Boundary<dim> &boundary,
                                         vector<Point<dim> > &ansatz_points) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                       const unsigned int           subface_no,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int          face_no,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int           subface_no,
-                                    const unsigned int           face_no,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -269,29 +124,6 @@ class FEQuadraticSub : public FiniteElement<dim> {
     virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
                                        const Boundary<dim> &boundary,
                                        dFMatrix &local_mass_matrix) const;
-
-  private:
-                                    /**
-                                     * 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;
 };
 
 
@@ -322,7 +154,7 @@ class FEQuadraticSub : public FiniteElement<dim> {
  * freedom.
  */
 template <int dim>
-class FECubicSub : public FiniteElement<dim> {
+class FECubicSub : public FELinearMapping<dim> {
   public:
                                     /**
                                      * Constructor
@@ -343,24 +175,6 @@ class FECubicSub : public FiniteElement<dim> {
     virtual Point<dim> shape_grad(const unsigned int i,
                                  const Point<dim>& p) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * For one dimensional elements, this
-                                     * function simply passes through to
-                                     * the one implemented in the base class.
-                                     */
-    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &ansatz_points,
-                                const bool           compute_ansatz_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const Boundary<dim> &boundary) const;
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -377,44 +191,6 @@ class FECubicSub : public FiniteElement<dim> {
                                         const Boundary<dim> &boundary,
                                         vector<Point<dim> > &ansatz_points) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                       const unsigned int           subface_no,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int          face_no,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int           subface_no,
-                                    const unsigned int           face_no,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -422,29 +198,6 @@ class FECubicSub : public FiniteElement<dim> {
     virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
                                        const Boundary<dim> &boundary,
                                        dFMatrix &local_mass_matrix) const;
-
-  private:
-                                    /**
-                                     * 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;
 };
 
 
@@ -476,7 +229,7 @@ class FECubicSub : public FiniteElement<dim> {
  * freedom.
  */
 template <int dim>
-class FEQuarticSub : public FiniteElement<dim> {
+class FEQuarticSub : public FELinearMapping<dim> {
   public:
                                     /**
                                      * Constructor
@@ -497,24 +250,6 @@ class FEQuarticSub : public FiniteElement<dim> {
     virtual Point<dim> shape_grad(const unsigned int i,
                                  const Point<dim>& p) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     *
-                                     * For one dimensional elements, this
-                                     * function simply passes through to
-                                     * the one implemented in the base class.
-                                     */
-    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                const vector<Point<dim> >            &unit_points,
-                                vector<dFMatrix>    &jacobians,
-                                const bool           compute_jacobians,
-                                vector<Point<dim> > &ansatz_points,
-                                const bool           compute_ansatz_points,
-                                vector<Point<dim> > &q_points,
-                                const bool           compute_q_points,
-                                const Boundary<dim> &boundary) const;
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -531,44 +266,6 @@ class FEQuarticSub : public FiniteElement<dim> {
                                         const Boundary<dim> &boundary,
                                         vector<Point<dim> > &ansatz_points) const;
 
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
-                                       const unsigned int           subface_no,
-                                       const vector<Point<dim-1> > &unit_points,
-                                       vector<double>      &face_jacobi_determinants) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int          face_no,
-                                    const Boundary<dim>         &boundary,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * information on this function.
-                                     */
-    virtual void get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int           subface_no,
-                                    const unsigned int           face_no,
-                                    const vector<Point<dim-1> > &unit_points,
-                                    vector<Point<dim> >         &normal_vectors) const;    
-
                                     /**
                                      * Refer to the base class for detailed
                                      * information on this function.
@@ -576,29 +273,6 @@ class FEQuarticSub : public FiniteElement<dim> {
     virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
                                        const Boundary<dim> &boundary,
                                        dFMatrix &local_mass_matrix) const;
-
-  private:
-                                    /**
-                                     * 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 e64d0b47105d5af0879ad9ae9a95c8a19f561a16..b8e20795d4fd3843a69c1f701c7f506ad3767df0 100644 (file)
@@ -362,10 +362,363 @@ void FiniteElement<dim>::get_ansatz_points (const DoFHandler<dim>::cell_iterator
 
 
 
+#if deal_II_dimension == 1
+
+template <>
+inline
+double
+FELinearMapping<1>::linear_shape_value(const unsigned int i,
+                                      const Point<1>     &p) const
+{
+  Assert((i<2), ExcInvalidIndex(i));
+  const double xi = p(0);
+  switch (i)
+    {
+      case 0: return 1.-xi;
+      case 1: return xi;
+    }
+  return 0.;
+};
+
+
+
+template <>
+inline
+Point<1>
+FELinearMapping<1>::linear_shape_grad(const unsigned int i,
+                                    const Point<1>&) const
+{
+  Assert((i<2), ExcInvalidIndex(i));
+  switch (i)
+    {
+    case 0: return Point<1>(-1.);
+    case 1: return Point<1>(1.);
+    }
+  return Point<1>();
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+                                           const Boundary<1>         &,
+                                           const vector<Point<0> > &,
+                                           vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+                                              const unsigned int           ,
+                                              const vector<Point<0> > &,
+                                              vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+                                           const unsigned int,
+                                           const Boundary<1> &,
+                                           const vector<Point<0> > &,
+                                           vector<Point<1> > &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+                                           const unsigned int,
+                                           const unsigned int,
+                                           const vector<Point<0> > &,
+                                           vector<Point<1> > &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+template <>
+void FELinearMapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+                                        const vector<Point<1> >            &unit_points,
+                                        vector<dFMatrix>  &jacobians,
+                                        const bool         compute_jacobians,
+                                        vector<Point<1> > &ansatz_points,
+                                        const bool         compute_ansatz_points,
+                                        vector<Point<1> > &q_points,
+                                        const bool         compute_q_points,
+                                        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);
+};
+
+
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+inline
+double
+FELinearMapping<2>::linear_shape_value (const unsigned int i,
+                                      const Point<2>& p) const
+{
+  Assert((i<4), ExcInvalidIndex(i));
+  switch (i)
+    {
+    case 0: return (1.-p(0)) * (1.-p(1));
+    case 1: return p(0) * (1.-p(1));
+    case 2: return p(0) * p(1);
+    case 3: return (1.-p(0)) * p(1);
+    }
+  return 0.;
+};
+
+
+
+template <>
+inline
+Point<2>
+FELinearMapping<2>::linear_shape_grad (const unsigned int i,
+                                     const Point<2>& p) const
+{
+  Assert((i<4), ExcInvalidIndex(i));
+  switch (i)
+    {
+    case 0: return Point<2> (p(1)-1., p(0)-1.);
+    case 1: return Point<2> (1.-p(1), -p(0));
+    case 2: return Point<2> (p(1), p(0));
+    case 3: return Point<2> (-p(1), 1.-p(0));
+    }
+  return Point<2> ();
+};
+
+
+
+template <>
+void FELinearMapping<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
+                                            const Boundary<2>         &,
+                                            const vector<Point<1> > &unit_points,
+                                            vector<double> &face_jacobians) const {
+                                  // more or less copied from the linear
+                                  // finite element
+  Assert (unit_points.size() == face_jacobians.size(),
+         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
+
+                                  // a linear mapping for a single line
+                                  // produces particularly simple
+                                  // expressions for the jacobi
+                                  // determinant :-)
+  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
+  fill_n (face_jacobians.begin(),
+         unit_points.size(),
+         h);  
+};
+
+
+
+template <>
+void FELinearMapping<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
+                                               const unsigned int           ,
+                                               const vector<Point<1> > &unit_points,
+                                               vector<double> &face_jacobians) const {
+                                  // more or less copied from the linear
+                                  // finite element
+  Assert (unit_points.size() == face_jacobians.size(),
+         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
+  Assert (face->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+
+                                  // a linear mapping for a single line
+                                  // produces particularly simple
+                                  // expressions for the jacobi
+                                  // determinant :-)
+  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
+  fill_n (face_jacobians.begin(),
+         unit_points.size(),
+         h/2);
+};
+
+
+
+template <>
+void FELinearMapping<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
+                                          const unsigned int       face_no,
+                                          const Boundary<2>       &,
+                                          const vector<Point<1> > &unit_points,
+                                          vector<Point<2> > &normal_vectors) const {
+                                  // more or less copied from the linear
+                                  // finite element
+  Assert (unit_points.size() == normal_vectors.size(),
+         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
+
+  const DoFHandler<2>::face_iterator face = cell->face(face_no);
+                                  // compute direction of line
+  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
+                                  // rotate to the right by 90 degrees
+  const Point<2> normal_direction(line_direction(1),
+                                 -line_direction(0));
+
+  if (face_no <= 1)
+                                    // for sides 0 and 1: return the correctly
+                                    // scaled vector
+    fill (normal_vectors.begin(), normal_vectors.end(),
+         normal_direction / sqrt(normal_direction.square()));
+  else
+                                    // for sides 2 and 3: scale and invert
+                                    // vector
+    fill (normal_vectors.begin(), normal_vectors.end(),
+         normal_direction / (-sqrt(normal_direction.square())));
+};
+
+
+
+template <>
+void FELinearMapping<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
+                                          const unsigned int       face_no,
+                                          const unsigned int,
+                                          const vector<Point<1> > &unit_points,
+                                          vector<Point<2> > &normal_vectors) const {
+                                  // more or less copied from the linear
+                                  // finite element
+                                  // note, that in 2D the normal vectors to the
+                                  // subface have the same direction as that
+                                  // for the face
+  Assert (unit_points.size() == normal_vectors.size(),
+         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
+  Assert (cell->face(face_no)->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+
+  const DoFHandler<2>::face_iterator face = cell->face(face_no);
+                                  // compute direction of line
+  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
+                                  // rotate to the right by 90 degrees
+  const Point<2> normal_direction(line_direction(1),
+                                 -line_direction(0));
+
+  if (face_no <= 1)
+                                    // for sides 0 and 1: return the correctly
+                                    // scaled vector
+    fill (normal_vectors.begin(), normal_vectors.end(),
+         normal_direction / sqrt(normal_direction.square()));
+  else
+                                    // for sides 2 and 3: scale and invert
+                                    // vector
+    fill (normal_vectors.begin(), normal_vectors.end(),
+         normal_direction / (-sqrt(normal_direction.square())));
+};
+
+#endif
+
+
+
+template <int dim>
+void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                          const vector<Point<dim> >            &unit_points,
+                                          vector<dFMatrix>    &jacobians,
+                                          const bool           compute_jacobians,
+                                          vector<Point<dim> > &ansatz_points,
+                                          const bool           compute_ansatz_points,
+                                          vector<Point<dim> > &q_points,
+                                          const bool           compute_q_points,
+                                          const Boundary<dim> &boundary) const {
+  Assert (jacobians.size() == unit_points.size(),
+         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+  Assert (q_points.size() == unit_points.size(),
+         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+  Assert (ansatz_points.size() == total_dofs,
+         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
+  
+  unsigned int n_points=unit_points.size();
+
+  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
+    vertices[l] = cell->vertex(l);
+  
+
+  if (compute_q_points) 
+    {
+                                      // initialize points to zero
+      for (unsigned int i=0; i<n_points; ++i)
+       q_points[i] = Point<dim> ();
+      
+                                      // note: let x_l be the vector of the
+                                      // lth quadrature point in real space and
+                                      // xi_l that on the unit cell, let further
+                                      // p_j be the vector of the jth vertex
+                                      // of the cell in real space and
+                                      // N_j(xi_l) be the value of the associated
+                                      // basis function at xi_l, then
+                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
+                                      //
+                                      // Here, N_j is the *linear* basis function,
+                                      // not that of the finite element, since we
+                                      // 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]);
+    };
+  
+
+/* jacobi matrices: compute d(x)/d(xi) and invert this
+   Let M(l) be the inverse of J at the quadrature point l, then
+     M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
+   where p_i(s) is the i-th coordinate of the s-th vertex vector,
+   N_s(l) is the value of the s-th vertex shape function at the
+   quadrature point l.
+
+   We could therefore write:
+   l=0..n_points-1
+     i=0..dim-1
+       j=0..dim-1
+         M_{ij}(l) = 0
+        s=0..n_vertices
+          M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
+
+  However, we rewrite the loops to only compute the gradient once for
+  each integration point and basis function.
+*/
+  if (compute_jacobians) 
+    {
+      dFMatrix M(dim,dim);
+      for (unsigned int l=0; l<n_points; ++l) 
+       {
+         M.clear ();
+         for (unsigned int s=0; s<GeometryInfo<dim>::vertices_per_cell; ++s)
+           {
+                                              // we want the linear transform,
+                                              // so use that function
+             const Point<dim> gradient = linear_shape_grad (s, unit_points[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);
+           };
+         jacobians[l].invert(M);
+       };
+    };
+
+                                  // compute ansatz points, which are
+                                  // the corners for linear elements
+  if (compute_ansatz_points)
+    get_ansatz_points (cell, boundary, ansatz_points);
+};
+
+
+
 /*------------------------------- 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 82730253c3f3efa3b00efc634abf19fce93a9897..2263f03d14d3c101e42d563aefbd5a64849488a9 100644 (file)
 
 template <>
 FECubicSub<1>::FECubicSub () :
-               FiniteElement<1> (1, 2) {
+               FELinearMapping<1> (1, 2) {
   prolongation[0](0,0) = 1.0;
   prolongation[0](0,1) = 0.0;
   prolongation[0](0,2) = 0.0;
@@ -340,25 +340,6 @@ FECubicSub<1>::FECubicSub () :
 
 
 
-template <>
-void FECubicSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
-                                    const vector<Point<1> >            &unit_points,
-                                    vector<dFMatrix>  &jacobians,
-                                    const bool         compute_jacobians,
-                                    vector<Point<1> > &ansatz_points,
-                                    const bool         compute_ansatz_points,
-                                    vector<Point<1> > &q_points,
-                                    const bool         compute_q_points,
-                                    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);
-};
-
-
-
 template <>
 double
 FECubicSub<1>::shape_value(const unsigned int i,
@@ -378,24 +359,6 @@ FECubicSub<1>::shape_value(const unsigned int i,
 
 
 
-template <>
-inline
-double
-FECubicSub<1>::linear_shape_value(const unsigned int i,
-                                     const Point<1>     &p) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  const double xi = p(0);
-  switch (i)
-    {
-      case 0: return 1.-xi;
-      case 1: return xi;
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<1>
 FECubicSub<1>::shape_grad(const unsigned int i,
@@ -415,23 +378,6 @@ FECubicSub<1>::shape_grad(const unsigned int i,
 
 
 
-template <>
-inline
-Point<1>
-FECubicSub<1>::linear_shape_grad(const unsigned int i,
-                                    const Point<1>&) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return Point<1>(-1.);
-    case 1: return Point<1>(1.);
-    }
-  return Point<1>();
-};
-
-
-
 template <>
 void FECubicSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
                                           const Boundary<1>  &boundary,
@@ -450,48 +396,6 @@ void FECubicSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_i
 
 
 
-template <>
-void FECubicSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                           const Boundary<1>         &,
-                                           const vector<Point<0> > &,
-                                           vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                              const unsigned int           ,
-                                              const vector<Point<0> > &,
-                                              vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                           const unsigned int,
-                                           const Boundary<1> &,
-                                           const vector<Point<0> > &,
-                                           vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FECubicSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                           const unsigned int,
-                                           const unsigned int,
-                                           const vector<Point<0> > &,
-                                           vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
 template <>
 void FECubicSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
                                               const Boundary<1> &,
@@ -537,7 +441,7 @@ void FECubicSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &c
 
 template <>
 FECubicSub<2>::FECubicSub () :
-               FiniteElement<2> (1, 2, 4)
+               FELinearMapping<2> (1, 2, 4)
 {
   interface_constraints(0,0) = -1.0/16.0;
   interface_constraints(0,1) = -1.0/16.0;
@@ -1023,25 +927,6 @@ xi*xi)*eta*eta*eta;
 
 
 
-template <>
-inline
-double
-FECubicSub<2>::linear_shape_value (const unsigned int i,
-                                      const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return (1.-p(0)) * (1.-p(1));
-    case 1: return p(0) * (1.-p(1));
-    case 2: return p(0) * p(1);
-    case 3: return (1.-p(0)) * p(1);
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<2>
 FECubicSub<2>::shape_grad (const unsigned int i,
@@ -1091,25 +976,6 @@ FECubicSub<2>::shape_grad (const unsigned int i,
 
 
 
-template <>
-inline
-Point<2>
-FECubicSub<2>::linear_shape_grad (const unsigned int i,
-                                     const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return Point<2> (p(1)-1., p(0)-1.);
-    case 1: return Point<2> (1.-p(1), -p(0));
-    case 2: return Point<2> (p(1), p(0));
-    case 3: return Point<2> (-p(1), 1.-p(0));
-    }
-  return Point<2> ();
-};
-
-
-
 template <>
 void FECubicSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
                                               const Boundary<2> &,
@@ -1763,220 +1629,12 @@ void FECubicSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_i
 
 
 
-template <>
-void FECubicSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
-                                           const Boundary<2>         &,
-                                           const vector<Point<1> > &unit_points,
-                                           vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h);  
-};
-
-
-
-template <>
-void FECubicSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
-                                             const unsigned int           ,
-                                             const vector<Point<1> > &unit_points,
-                                             vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-  Assert (face->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h/2);
-};
-
-
-
-template <>
-void FECubicSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                          const unsigned int       face_no,
-                                          const Boundary<2>       &,
-                                          const vector<Point<1> > &unit_points,
-                                          vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
-
-
-template <>
-void FECubicSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                          const unsigned int       face_no,
-                                          const unsigned int,
-                                          const vector<Point<1> > &unit_points,
-                                          vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-                                  // note, that in 2D the normal vectors to the
-                                  // subface have the same direction as that
-                                  // for the face
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-  Assert (cell->face(face_no)->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
 #endif
 
 
 
 
 
-template <int dim>
-void FECubicSub<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                         const vector<Point<dim> >            &unit_points,
-                                         vector<dFMatrix>    &jacobians,
-                                         const bool           compute_jacobians,
-                                         vector<Point<dim> > &ansatz_points,
-                                         const bool           compute_ansatz_points,
-                                         vector<Point<dim> > &q_points,
-                                         const bool           compute_q_points,
-                                         const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
-         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
-         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (ansatz_points.size() == total_dofs,
-         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
-
-  
-  unsigned int n_points=unit_points.size();
-
-  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
-  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
-    vertices[l] = cell->vertex(l);
-  
-
-  if (compute_q_points) 
-    {
-                                      // initialize points to zero
-      for (unsigned int i=0; i<n_points; ++i)
-       q_points[i] = Point<dim> ();
-      
-                                      // note: let x_l be the vector of the
-                                      // lth quadrature point in real space and
-                                      // xi_l that on the unit cell, let further
-                                      // p_j be the vector of the jth vertex
-                                      // of the cell in real space and
-                                      // N_j(xi_l) be the value of the associated
-                                      // basis function at xi_l, then
-                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
-                                      //
-                                      // Here, N_j is the *linear* basis function,
-                                      // not that of the finite element, since we
-                                      // 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]);
-    };
-  
-
-/* jacobi matrices: compute d(x)/d(xi) and invert this
-   Let M(l) be the inverse of J at the quadrature point l, then
-     M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
-   where p_i(s) is the i-th coordinate of the s-th vertex vector,
-   N_s(l) is the value of the s-th vertex shape function at the
-   quadrature point l (linear shape functions implied, as these
-   are used for the mapping).
-
-   We could therefore write:
-   l=0..n_points-1
-     i=0..dim-1
-       j=0..dim-1
-         M_{ij}(l) = 0
-        s=0..n_vertices
-          M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
-
-  However, we rewrite the loops to only compute the gradient once for
-  each integration point and basis function.
-*/
-  if (compute_jacobians) 
-    {
-      dFMatrix M(dim,dim);
-      for (unsigned int l=0; l<n_points; ++l) 
-       {
-         M.clear ();
-         for (unsigned int s=0; s<GeometryInfo<dim>::vertices_per_cell; ++s)
-           {
-                                              // we want the linear transform,
-                                              // so use that function
-             const Point<dim> gradient = linear_shape_grad (s, unit_points[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);
-           };
-         jacobians[l].invert(M);
-       };
-    };
-
-                                  // compute ansatz points, which are
-                                  // the corners for linear elements
-  if (compute_ansatz_points)
-    get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
 
 // explicit instantiations
 
index a58dd9cd1238ef391d4bd38150206dee7a33c0b4..f293ae8dc9b899772466b529f68b3d4a3c138243 100644 (file)
@@ -13,7 +13,7 @@
 
 template <>
 FELinear<1>::FELinear () :
-               FiniteElement<1> (1, 0)
+               FELinearMapping<1> (1, 0)
 {
                                   // for restriction and prolongation matrices:
                                   // note that we do not add up all the
@@ -79,25 +79,6 @@ FELinear<1>::shape_grad(const unsigned int i,
 
 
 
-template <>
-void FELinear<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
-                                 const vector<Point<1> >            &unit_points,
-                                 vector<dFMatrix>  &jacobians,
-                                 const bool         compute_jacobians,
-                                 vector<Point<1> > &ansatz_points,
-                                 const bool         compute_ansatz_points,
-                                 vector<Point<1> > &q_points,
-                                 const bool         compute_q_points,
-                                 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);
-};
-
-
-
 template <>
 void FELinear<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
                                     const Boundary<1>  &boundary,
@@ -116,48 +97,6 @@ void FELinear<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_ite
 
 
 
-template <>
-void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                     const Boundary<1>         &,
-                                     const vector<Point<0> > &,
-                                     vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                        const unsigned int           ,
-                                        const vector<Point<0> > &,
-                                        vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                     const unsigned int,
-                                     const Boundary<1> &,
-                                     const vector<Point<0> > &,
-                                     vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                     const unsigned int,
-                                     const unsigned int,
-                                     const vector<Point<0> > &,
-                                     vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
 template <>
 void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
                                         const Boundary<1> &,
@@ -183,7 +122,7 @@ void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cel
 
 template <>
 FELinear<2>::FELinear () :
-               FiniteElement<2> (1, 0, 0)
+               FELinearMapping<2> (1, 0, 0)
 {
   interface_constraints(0,0) = 1./2.;
   interface_constraints(0,1) = 1./2.;
@@ -413,209 +352,10 @@ void FELinear<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cel
 
 
 
-template <>
-void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
-                                     const Boundary<2>         &,
-                                     const vector<Point<1> > &unit_points,
-                                     vector<double>      &face_jacobians) const {
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h);  
-};
-
-
-
-template <>
-void FELinear<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
-                                        const unsigned int,
-                                        const vector<Point<1> > &unit_points,
-                                        vector<double>      &face_jacobians) const {
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-  Assert (face->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h/2);
-};
-
-
-
-template <>
-void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                     const unsigned int       face_no,
-                                     const Boundary<2> &,
-                                     const vector<Point<1> >  &unit_points,
-                                     vector<Point<2> >        &normal_vectors) const {
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
-
-
-template <>
-void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                     const unsigned int       face_no,
-                                     const unsigned int,
-                                     const vector<Point<1> >  &unit_points,
-                                     vector<Point<2> >        &normal_vectors) const {
-                                  // note, that in 2D the normal vectors to the
-                                  // subface have the same direction as that
-                                  // for the face
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-  Assert (cell->face(face_no)->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
 #endif
 
 
 
-template <int dim>
-void FELinear<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                   const vector<Point<dim> >            &unit_points,
-                                   vector<dFMatrix>    &jacobians,
-                                   const bool           compute_jacobians,
-                                   vector<Point<dim> > &ansatz_points,
-                                   const bool           compute_ansatz_points,
-                                   vector<Point<dim> > &q_points,
-                                   const bool           compute_q_points,
-                                   const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
-         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
-         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (ansatz_points.size() == total_dofs,
-         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
-  
-  unsigned int n_points=unit_points.size();
-
-  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
-  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
-    vertices[l] = cell->vertex(l);
-  
-
-  if (compute_q_points) 
-    {
-                                      // initialize points to zero
-      for (unsigned int i=0; i<n_points; ++i)
-       q_points[i] = Point<dim> ();
-      
-                                      // note: let x_l be the vector of the
-                                      // lth quadrature point in real space and
-                                      // xi_l that on the unit cell, let further
-                                      // p_j be the vector of the jth vertex
-                                      // of the cell in real space and
-                                      // N_j(xi_l) be the value of the associated
-                                      // basis function at xi_l, then
-                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
-      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] *
-                        FELinear<dim>::shape_value(j, unit_points[l]);
-    };
-  
-
-/* jacobi matrices: compute d(x)/d(xi) and invert this
-   Let M(l) be the inverse of J at the quadrature point l, then
-     M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
-   where p_i(s) is the i-th coordinate of the s-th vertex vector,
-   N_s(l) is the value of the s-th vertex shape function at the
-   quadrature point l.
-
-   We could therefore write:
-   l=0..n_points-1
-     i=0..dim-1
-       j=0..dim-1
-         M_{ij}(l) = 0
-        s=0..n_vertices
-          M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
-
-  However, we rewrite the loops to only compute the gradient once for
-  each integration point and basis function.
-*/
-  if (compute_jacobians) 
-    {
-      dFMatrix M(dim,dim);
-      for (unsigned int l=0; l<n_points; ++l) 
-       {
-         M.clear ();
-         for (unsigned int s=0; s<GeometryInfo<dim>::vertices_per_cell; ++s)
-           {
-                                              // we want a linear transform and
-                                              // if we prepend the class name in
-                                              // front of the #shape_grad#, we
-                                              // need not use virtual function
-                                              // calls.
-             const Point<dim> gradient
-               = FELinear<dim>::shape_grad (s, unit_points[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);
-           };
-         jacobians[l].invert(M);
-       };
-    };
-
-                                  // compute ansatz points, which are
-                                  // the corners for linear elements
-  if (compute_ansatz_points)
-    get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
 template <int dim>
 void FELinear<dim>::get_ansatz_points (const typename DoFHandler<dim>::cell_iterator &cell,
                                       const Boundary<dim>  &,
index 52884f1df1c81ab8191a4ffd879202ad3a8e720b..2f10e9d1a8b7a968c17bce6be83e62224b35f090 100644 (file)
@@ -15,7 +15,7 @@
 
 template <>
 FEQuadraticSub<1>::FEQuadraticSub () :
-               FiniteElement<1> (1, 1) {
+               FELinearMapping<1> (1, 1) {
 /*
   Get the prolongation matrices by the following little maple script:
 
@@ -62,25 +62,6 @@ FEQuadraticSub<1>::FEQuadraticSub () :
 
 
 
-template <>
-void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
-                                    const vector<Point<1> >            &unit_points,
-                                    vector<dFMatrix>  &jacobians,
-                                    const bool         compute_jacobians,
-                                    vector<Point<1> > &ansatz_points,
-                                    const bool         compute_ansatz_points,
-                                    vector<Point<1> > &q_points,
-                                    const bool         compute_q_points,
-                                    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);
-};
-
-
-
 template <>
 double
 FEQuadraticSub<1>::shape_value(const unsigned int i,
@@ -99,24 +80,6 @@ FEQuadraticSub<1>::shape_value(const unsigned int i,
 
 
 
-template <>
-inline
-double
-FEQuadraticSub<1>::linear_shape_value(const unsigned int i,
-                                     const Point<1>     &p) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  const double xi = p(0);
-  switch (i)
-    {
-      case 0: return 1.-xi;
-      case 1: return xi;
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<1>
 FEQuadraticSub<1>::shape_grad(const unsigned int i,
@@ -135,23 +98,6 @@ FEQuadraticSub<1>::shape_grad(const unsigned int i,
 
 
 
-template <>
-inline
-Point<1>
-FEQuadraticSub<1>::linear_shape_grad(const unsigned int i,
-                                    const Point<1>&) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return Point<1>(-1.);
-    case 1: return Point<1>(1.);
-    }
-  return Point<1>();
-};
-
-
-
 template <>
 void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
                                           const Boundary<1>  &boundary,
@@ -170,48 +116,6 @@ void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::fa
 
 
 
-template <>
-void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                           const Boundary<1>         &,
-                                           const vector<Point<0> > &,
-                                           vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                              const unsigned int           ,
-                                              const vector<Point<0> > &,
-                                              vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                           const unsigned int,
-                                           const Boundary<1> &,
-                                           const vector<Point<0> > &,
-                                           vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                           const unsigned int,
-                                           const unsigned int,
-                                           const vector<Point<0> > &,
-                                           vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
 template <>
 void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
                                               const Boundary<1> &,
@@ -239,7 +143,7 @@ void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterato
 
 template <>
 FEQuadraticSub<2>::FEQuadraticSub () :
-               FiniteElement<2> (1, 1, 1)
+               FELinearMapping<2> (1, 1, 1)
 {
   interface_constraints(0,2) = 1.0;
   interface_constraints(1,0) = 3./8.;
@@ -641,25 +545,6 @@ FEQuadraticSub<2>::shape_value (const unsigned int i,
 
 
 
-template <>
-inline
-double
-FEQuadraticSub<2>::linear_shape_value (const unsigned int i,
-                                      const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return (1.-p(0)) * (1.-p(1));
-    case 1: return p(0) * (1.-p(1));
-    case 2: return p(0) * p(1);
-    case 3: return (1.-p(0)) * p(1);
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<2>
 FEQuadraticSub<2>::shape_grad (const unsigned int i,
@@ -695,25 +580,6 @@ FEQuadraticSub<2>::shape_grad (const unsigned int i,
 
 
 
-template <>
-inline
-Point<2>
-FEQuadraticSub<2>::linear_shape_grad (const unsigned int i,
-                                     const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-    case 0: return Point<2> (p(1)-1., p(0)-1.);
-    case 1: return Point<2> (1.-p(1), -p(0));
-    case 2: return Point<2> (p(1), p(0));
-    case 3: return Point<2> (-p(1), 1.-p(0));
-    }
-  return Point<2> ();
-};
-
-
-
 template <>
 void FEQuadraticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
                                               const Boundary<2> &,
@@ -1033,220 +899,12 @@ void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::fa
 
 
 
-template <>
-void FEQuadraticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
-                                           const Boundary<2>         &,
-                                           const vector<Point<1> > &unit_points,
-                                           vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h);  
-};
-
-
-
-template <>
-void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
-                                             const unsigned int           ,
-                                             const vector<Point<1> > &unit_points,
-                                             vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-  Assert (face->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h/2);
-};
-
-
-
-template <>
-void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                          const unsigned int       face_no,
-                                          const Boundary<2>       &,
-                                          const vector<Point<1> > &unit_points,
-                                          vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
-
-
-template <>
-void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                          const unsigned int       face_no,
-                                          const unsigned int,
-                                          const vector<Point<1> > &unit_points,
-                                          vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-                                  // note, that in 2D the normal vectors to the
-                                  // subface have the same direction as that
-                                  // for the face
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-  Assert (cell->face(face_no)->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
 #endif
 
 
 
 
 
-template <int dim>
-void FEQuadraticSub<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                         const vector<Point<dim> >            &unit_points,
-                                         vector<dFMatrix>    &jacobians,
-                                         const bool           compute_jacobians,
-                                         vector<Point<dim> > &ansatz_points,
-                                         const bool           compute_ansatz_points,
-                                         vector<Point<dim> > &q_points,
-                                         const bool           compute_q_points,
-                                         const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
-         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
-         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (ansatz_points.size() == total_dofs,
-         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
-
-  
-  unsigned int n_points=unit_points.size();
-
-  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
-  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
-    vertices[l] = cell->vertex(l);
-  
-
-  if (compute_q_points) 
-    {
-                                      // initialize points to zero
-      for (unsigned int i=0; i<n_points; ++i)
-       q_points[i] = Point<dim> ();
-      
-                                      // note: let x_l be the vector of the
-                                      // lth quadrature point in real space and
-                                      // xi_l that on the unit cell, let further
-                                      // p_j be the vector of the jth vertex
-                                      // of the cell in real space and
-                                      // N_j(xi_l) be the value of the associated
-                                      // basis function at xi_l, then
-                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
-                                      //
-                                      // Here, N_j is the *linear* basis function,
-                                      // not that of the finite element, since we
-                                      // 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]);
-    };
-  
-
-/* jacobi matrices: compute d(x)/d(xi) and invert this
-   Let M(l) be the inverse of J at the quadrature point l, then
-     M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
-   where p_i(s) is the i-th coordinate of the s-th vertex vector,
-   N_s(l) is the value of the s-th vertex shape function at the
-   quadrature point l.
-
-   We could therefore write:
-   l=0..n_points-1
-     i=0..dim-1
-       j=0..dim-1
-         M_{ij}(l) = 0
-        s=0..n_vertices
-          M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
-
-  However, we rewrite the loops to only compute the gradient once for
-  each integration point and basis function.
-*/
-  if (compute_jacobians) 
-    {
-      dFMatrix M(dim,dim);
-      for (unsigned int l=0; l<n_points; ++l) 
-       {
-         M.clear ();
-         for (unsigned int s=0; s<GeometryInfo<dim>::vertices_per_cell; ++s)
-           {
-                                              // we want the linear transform,
-                                              // so use that function
-             const Point<dim> gradient = linear_shape_grad (s, unit_points[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);
-           };
-         jacobians[l].invert(M);
-       };
-    };
-
-                                  // compute ansatz points, which are
-                                  // the corners for linear elements
-  if (compute_ansatz_points)
-    get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
-
-
 // explicit instantiations
 
 template class FEQuadraticSub<deal_II_dimension>;
index fb875967895795681c24e46604e8ea50483c519f..71ad977592173ecc3e30c4ba2532ce9f52b810db 100644 (file)
 
 template <>
 FEQuarticSub<1>::FEQuarticSub () :
-               FiniteElement<1> (1, 3) {
+               FELinearMapping<1> (1, 3) {
   prolongation[0](0,0) = 1.0;
   prolongation[0](1,3) = 1.0;
   prolongation[0](2,0) = 35.0/128.0;
@@ -351,25 +351,6 @@ FEQuarticSub<1>::FEQuarticSub () :
 
 
 
-template <>
-void FEQuarticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
-                                     const vector<Point<1> >            &unit_points,
-                                     vector<dFMatrix>  &jacobians,
-                                     const bool         compute_jacobians,
-                                     vector<Point<1> > &ansatz_points,
-                                     const bool         compute_ansatz_points,
-                                     vector<Point<1> > &q_points,
-                                     const bool         compute_q_points,
-                                     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);
-};
-
-
-
 template <>
 double
 FEQuarticSub<1>::shape_value(const unsigned int i,
@@ -390,24 +371,6 @@ FEQuarticSub<1>::shape_value(const unsigned int i,
 
 
 
-template <>
-inline
-double
-FEQuarticSub<1>::linear_shape_value(const unsigned int i,
-                                   const Point<1>     &p) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  const double xi = p(0);
-  switch (i)
-    {
-      case 0: return 1.-xi;
-      case 1: return xi;
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<1>
 FEQuarticSub<1>::shape_grad(const unsigned int i,
@@ -428,23 +391,6 @@ FEQuarticSub<1>::shape_grad(const unsigned int i,
 
 
 
-template <>
-inline
-Point<1>
-FEQuarticSub<1>::linear_shape_grad(const unsigned int i,
-                                  const Point<1>&) const
-{
-  Assert((i<2), ExcInvalidIndex(i));
-  switch (i)
-    {
-      case 0: return Point<1>(-1.);
-      case 1: return Point<1>(1.);
-    }
-  return Point<1>();
-};
-
-
-
 template <>
 void FEQuarticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
                                         const Boundary<1>  &boundary,
@@ -463,48 +409,6 @@ void FEQuarticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face
 
 
 
-template <>
-void FEQuarticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                         const Boundary<1>         &,
-                                         const vector<Point<0> > &,
-                                         vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                            const unsigned int           ,
-                                            const vector<Point<0> > &,
-                                            vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                         const unsigned int,
-                                         const Boundary<1> &,
-                                         const vector<Point<0> > &,
-                                         vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FEQuarticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                         const unsigned int,
-                                         const unsigned int,
-                                         const vector<Point<0> > &,
-                                         vector<Point<1> > &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
 template <>
 void FEQuarticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
                                             const Boundary<1> &,
@@ -561,7 +465,7 @@ void FEQuarticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator
 
 template <>
 FEQuarticSub<2>::FEQuarticSub () :
-               FiniteElement<2> (1, 3, 9)
+               FELinearMapping<2> (1, 3, 9)
 {
   interface_constraints(0,3) = 1.0;
   interface_constraints(1,0) = 35.0/128.0;
@@ -1414,25 +1318,6 @@ FEQuarticSub<2>::shape_value (const unsigned int i,
 
 
 
-template <>
-inline
-double
-FEQuarticSub<2>::linear_shape_value (const unsigned int i,
-                                    const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-      case 0: return (1.-p(0)) * (1.-p(1));
-      case 1: return p(0) * (1.-p(1));
-      case 2: return p(0) * p(1);
-      case 3: return (1.-p(0)) * p(1);
-    }
-  return 0.;
-};
-
-
-
 template <>
 Point<2>
 FEQuarticSub<2>::shape_grad (const unsigned int i,
@@ -1500,25 +1385,6 @@ FEQuarticSub<2>::shape_grad (const unsigned int i,
 
 
 
-template <>
-inline
-Point<2>
-FEQuarticSub<2>::linear_shape_grad (const unsigned int i,
-                                   const Point<2>& p) const
-{
-  Assert((i<4), ExcInvalidIndex(i));
-  switch (i)
-    {
-      case 0: return Point<2> (p(1)-1., p(0)-1.);
-      case 1: return Point<2> (1.-p(1), -p(0));
-      case 2: return Point<2> (p(1), p(0));
-      case 3: return Point<2> (-p(1), 1.-p(0));
-    }
-  return Point<2> ();
-};
-
-
-
 template <>
 void FEQuarticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
                                             const Boundary<2> &,
@@ -2912,219 +2778,12 @@ void FEQuarticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face
 
 
 
-template <>
-void FEQuarticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
-                                         const Boundary<2>         &,
-                                         const vector<Point<1> > &unit_points,
-                                         vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h);  
-};
-
-
-
-template <>
-void FEQuarticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
-                                            const unsigned int           ,
-                                            const vector<Point<1> > &unit_points,
-                                            vector<double> &face_jacobians) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == face_jacobians.size(),
-         ExcWrongFieldDimension (unit_points.size(), face_jacobians.size()));
-  Assert (face->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-                                  // a linear mapping for a single line
-                                  // produces particularly simple
-                                  // expressions for the jacobi
-                                  // determinant :-)
-  const double h = sqrt((face->vertex(1) - face->vertex(0)).square());
-  fill_n (face_jacobians.begin(),
-         unit_points.size(),
-         h/2);
-};
-
-
-
-template <>
-void FEQuarticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                         const unsigned int       face_no,
-                                         const Boundary<2>       &,
-                                         const vector<Point<1> > &unit_points,
-                                         vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
-
-
-template <>
-void FEQuarticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
-                                         const unsigned int       face_no,
-                                         const unsigned int,
-                                         const vector<Point<1> > &unit_points,
-                                         vector<Point<2> > &normal_vectors) const {
-                                  // more or less copied from the linear
-                                  // finite element
-                                  // note, that in 2D the normal vectors to the
-                                  // subface have the same direction as that
-                                  // for the face
-  Assert (unit_points.size() == normal_vectors.size(),
-         ExcWrongFieldDimension (unit_points.size(), normal_vectors.size()));
-  Assert (cell->face(face_no)->at_boundary() == false,
-         ExcBoundaryFaceUsed ());
-
-  const DoFHandler<2>::face_iterator face = cell->face(face_no);
-                                  // compute direction of line
-  const Point<2> line_direction = (face->vertex(1) - face->vertex(0));
-                                  // rotate to the right by 90 degrees
-  const Point<2> normal_direction(line_direction(1),
-                                 -line_direction(0));
-
-  if (face_no <= 1)
-                                    // for sides 0 and 1: return the correctly
-                                    // scaled vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / sqrt(normal_direction.square()));
-  else
-                                    // for sides 2 and 3: scale and invert
-                                    // vector
-    fill (normal_vectors.begin(), normal_vectors.end(),
-         normal_direction / (-sqrt(normal_direction.square())));
-};
-
 #endif
 
 
 
 
 
-template <int dim>
-void FEQuarticSub<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
-                                       const vector<Point<dim> >            &unit_points,
-                                       vector<dFMatrix>    &jacobians,
-                                       const bool           compute_jacobians,
-                                       vector<Point<dim> > &ansatz_points,
-                                       const bool           compute_ansatz_points,
-                                       vector<Point<dim> > &q_points,
-                                       const bool           compute_q_points,
-                                       const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
-         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
-         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (ansatz_points.size() == total_dofs,
-         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
-
-  
-  unsigned int n_points=unit_points.size();
-
-  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
-  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
-    vertices[l] = cell->vertex(l);
-  
-
-  if (compute_q_points) 
-    {
-                                      // initialize points to zero
-      for (unsigned int i=0; i<n_points; ++i)
-       q_points[i] = Point<dim> ();
-      
-                                      // note: let x_l be the vector of the
-                                      // lth quadrature point in real space and
-                                      // xi_l that on the unit cell, let further
-                                      // p_j be the vector of the jth vertex
-                                      // of the cell in real space and
-                                      // N_j(xi_l) be the value of the associated
-                                      // basis function at xi_l, then
-                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
-                                      //
-                                      // Here, N_j is the *linear* basis function,
-                                      // not that of the finite element, since we
-                                      // 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]);
-    };
-  
-
-/* jacobi matrices: compute d(x)/d(xi) and invert this
-   Let M(l) be the inverse of J at the quadrature point l, then
-   M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
-   where p_i(s) is the i-th coordinate of the s-th vertex vector,
-   N_s(l) is the value of the s-th vertex shape function at the
-   quadrature point l (linear shape functions implied, as these
-   are used for the mapping).
-
-   We could therefore write:
-   l=0..n_points-1
-   i=0..dim-1
-   j=0..dim-1
-   M_{ij}(l) = 0
-   s=0..n_vertices
-   M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
-
-   However, we rewrite the loops to only compute the gradient once for
-   each integration point and basis function.
-*/
-  if (compute_jacobians) 
-    {
-      dFMatrix M(dim,dim);
-      for (unsigned int l=0; l<n_points; ++l) 
-       {
-         M.clear ();
-         for (unsigned int s=0; s<GeometryInfo<dim>::vertices_per_cell; ++s)
-           {
-                                              // we want the linear transform,
-                                              // so use that function
-             const Point<dim> gradient = linear_shape_grad (s, unit_points[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);
-           };
-         jacobians[l].invert(M);
-       };
-    };
-
-                                  // compute ansatz points, which are
-                                  // the corners for linear elements
-  if (compute_ansatz_points)
-    get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
 
 
 // explicit instantiations

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.