]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Separate the linear mapping from the rest of the finite element stuff. This is a...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 22 Jul 1998 08:26:16 +0000 (08:26 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 22 Jul 1998 08:26:16 +0000 (08:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@452 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/include/fe/fe_linear_mapping.h [new file with mode: 0644]
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_linear_mapping.cc [new file with mode: 0644]

index 0dcc0001f46c7c60f79c77363700dbfe6de69e85..fbe50a070f1845d3320a7e6ca971bbbb29c37140 100644 (file)
@@ -1041,143 +1041,6 @@ 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,
-                                       GeometryInfo<dim>::vertices_per_cell) {};
-
-                                    /**
-                                     * Return the value of the #i#th shape
-                                     * function at point #p# on the unit cell.
-                                     * Here, the (bi-)linear basis functions
-                                     * are meant, which are used for the
-                                     * computation of the transformation from
-                                     * unit cell to real space cell.
-                                     */
-    virtual double shape_value_transform (const unsigned int i,
-                                         const Point<dim> &p) const;
-
-                                    /**
-                                     * Return the gradient of the #i#th shape
-                                     * function at point #p# on the unit cell.
-                                     * Here, the (bi-)linear basis functions
-                                     * are meant, which are used for the
-                                     * computation of the transformation from
-                                     * unit cell to real space cell.
-                                     */
-    virtual Point<dim> shape_grad_transform (const unsigned int i,
-                                            const Point<dim> &p) const;
-
-                                    /**
-                                     * Refer to the base class for detailed
-                                     * 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 dFMatrix      &shape_values_transform,
-                                const vector<vector<Point<dim> > > &shape_grad_transform,
-                                const Boundary<dim> &boundary) const;
-};
-
-
-
   
 /*----------------------------   fe.h     ---------------------------*/
 /* end of #ifndef __fe_H */
index 608dd3345f852c535e9b632ca41232c3028d1507..22da61caf836a490d116866c060d6b53f38307aa 100644 (file)
@@ -5,7 +5,7 @@
 /*----------------------------   fe_lib.h     ---------------------------*/
 
 
-#include <fe/fe.h>
+#include <fe/fe_linear_mapping.h>
 
 
 
diff --git a/deal.II/deal.II/include/fe/fe_linear_mapping.h b/deal.II/deal.II/include/fe/fe_linear_mapping.h
new file mode 100644 (file)
index 0000000..86f2935
--- /dev/null
@@ -0,0 +1,154 @@
+/*----------------------------   fe_linear_mapping.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __fe_linear_mapping_H
+#define __fe_linear_mapping_H
+/*----------------------------   fe_linear_mapping.h     ---------------------------*/
+
+
+#include <fe/fe.h>
+
+
+
+
+/**
+ * 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,
+                                       GeometryInfo<dim>::vertices_per_cell) {};
+
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at point #p# on the unit cell.
+                                     * Here, the (bi-)linear basis functions
+                                     * are meant, which are used for the
+                                     * computation of the transformation from
+                                     * unit cell to real space cell.
+                                     */
+    virtual double shape_value_transform (const unsigned int i,
+                                         const Point<dim> &p) const;
+
+                                    /**
+                                     * Return the gradient of the #i#th shape
+                                     * function at point #p# on the unit cell.
+                                     * Here, the (bi-)linear basis functions
+                                     * are meant, which are used for the
+                                     * computation of the transformation from
+                                     * unit cell to real space cell.
+                                     */
+    virtual Point<dim> shape_grad_transform (const unsigned int i,
+                                            const Point<dim> &p) const;
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * 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 dFMatrix      &shape_values_transform,
+                                const vector<vector<Point<dim> > > &shape_grad_transform,
+                                const Boundary<dim> &boundary) const;
+};
+
+
+
+
+/*----------------------------   fe_linear_mapping.h     ---------------------------*/
+/* end of #ifndef __fe_linear_mapping_H */
+#endif
+/*----------------------------   fe_linear_mapping.h     ---------------------------*/
index ece34091c5f15c95bf360a826a233bf57b2520e9..449de51e7c30df641709648bc1e3bdcdb59d3309 100644 (file)
@@ -434,368 +434,11 @@ void FiniteElement<dim>::get_ansatz_points (const DoFHandler<dim>::cell_iterator
 
 
 
-/*---------------------------- FELinearMapping ----------------------------------*/
-
-#if deal_II_dimension == 1
-
-template <>
-inline
-double
-FELinearMapping<1>::shape_value_transform (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>::shape_grad_transform(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 FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                           const Boundary<1>         &,
-                                           const vector<Point<0> > &,
-                                           vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                              const unsigned int           ,
-                                              const vector<Point<0> > &,
-                                              vector<double>      &) const {
-  Assert (false, ExcInternalError());
-};
-
-
-
-template <>
-void FELinearMapping<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 FELinearMapping<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 dFMatrix      &shape_values_transform,
-                                        const vector<vector<Point<dim> > > &shape_gradients_transform,
-                                        const Boundary<1> &boundary) const {
-                                  // simply pass down
-  FiniteElement<1>::fill_fe_values (cell, unit_points,
-                                   jacobians, compute_jacobians,
-                                   ansatz_points, compute_ansatz_points,
-                                   q_points, compute_q_points,
-                                   shape_values_transform, shape_gradients_transform,
-                                   boundary);
-};
-
-
-
-#endif
-
-
-
-#if deal_II_dimension == 2
-
-template <>
-inline
-double
-FELinearMapping<2>::shape_value_transform (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>::shape_grad_transform (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 dFMatrix      &shape_values_transform,
-                                          const vector<vector<Point<dim> > > &shape_grad_transform,
-                                          const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
-         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  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] * shape_values_transform(j, 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 = shape_grad_transform[s][l];
-             for (unsigned int i=0; i<dim; ++i)
-               for (unsigned int j=0; j<dim; ++j)
-                 M(i,j) += vertices[s](i) * gradient(j);
-           };
-         jacobians[l].invert(M);
-       };
-    };
-
-  if (compute_ansatz_points)
-    get_ansatz_points (cell, boundary, ansatz_points);
-};
-
-
 
 /*------------------------------- Explicit Instantiations -------------*/
 
 template class FiniteElementBase<deal_II_dimension>;
 template class FiniteElement<deal_II_dimension>;
-template class FELinearMapping<deal_II_dimension>;
+
 
 
diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc
new file mode 100644 (file)
index 0000000..4f86319
--- /dev/null
@@ -0,0 +1,374 @@
+/* $Id$ */
+
+#include <fe/fe_linear_mapping.h>
+#include <fe/quadrature.h>
+#include <grid/tria_iterator.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_boundary.h>
+
+
+
+
+
+/*---------------------------- FELinearMapping ----------------------------------*/
+
+#if deal_II_dimension == 1
+
+template <>
+inline
+double
+FELinearMapping<1>::shape_value_transform (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>::shape_grad_transform(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 FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+                                           const Boundary<1>         &,
+                                           const vector<Point<0> > &,
+                                           vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+                                              const unsigned int           ,
+                                              const vector<Point<0> > &,
+                                              vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <>
+void FELinearMapping<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 FELinearMapping<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 dFMatrix      &shape_values_transform,
+                                        const vector<vector<Point<dim> > > &shape_gradients_transform,
+                                        const Boundary<1> &boundary) const {
+                                  // simply pass down
+  FiniteElement<1>::fill_fe_values (cell, unit_points,
+                                   jacobians, compute_jacobians,
+                                   ansatz_points, compute_ansatz_points,
+                                   q_points, compute_q_points,
+                                   shape_values_transform, shape_gradients_transform,
+                                   boundary);
+};
+
+
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+inline
+double
+FELinearMapping<2>::shape_value_transform (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>::shape_grad_transform (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 dFMatrix      &shape_values_transform,
+                                          const vector<vector<Point<dim> > > &shape_grad_transform,
+                                          const Boundary<dim> &boundary) const {
+  Assert (jacobians.size() == unit_points.size(),
+         ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+  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] * shape_values_transform(j, 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 = shape_grad_transform[s][l];
+             for (unsigned int i=0; i<dim; ++i)
+               for (unsigned int j=0; j<dim; ++j)
+                 M(i,j) += vertices[s](i) * gradient(j);
+           };
+         jacobians[l].invert(M);
+       };
+    };
+
+  if (compute_ansatz_points)
+    get_ansatz_points (cell, boundary, ansatz_points);
+};
+
+
+
+
+/*------------------------------- Explicit Instantiations -------------*/
+
+template class FELinearMapping<deal_II_dimension>;

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.