]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Still changes for the support of integration along faces.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 27 Apr 1998 18:03:45 +0000 (18:03 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 27 Apr 1998 18:03:45 +0000 (18:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@205 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_update_flags.h
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/fe/fe_values.cc

index 188a130e5cbd83eea2c4f95087276424503cdaa4..142b1b7cd4118ac11aab1eb6f7b02b001a28d613 100644 (file)
@@ -546,15 +546,16 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * performed by the specific finite element
                                      * class.
                                      *
-                                     * Two fields remain to be finite element
+                                     * Three fields remain to be finite element
                                      * specific in this standard implementation:
                                      * The jacobi determinants of the
                                      * transformation from the unit face to the
-                                     * real face and the ansatz points. For
-                                     * these two fields, there exist pure
-                                     * virtual functions, #get_face_jacobians#
-                                     * and the #get_face_ansatz_points#
-                                     * function.
+                                     * real face, the ansatz points
+                                     * and the outward normal vectors. For
+                                     * these fields, there exist pure
+                                     * virtual functions, #get_face_jacobians#,
+                                     * #get_face_ansatz_points# and
+                                     * #get_normal_vectors#.
                                      *
                                      * Though there is a standard
                                      * implementation, there
@@ -587,6 +588,8 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      const bool           compute_q_points,
                                      vector<double>      &face_jacobi_determinants,
                                      const bool           compute_face_jacobians,
+                                     vector<Point<dim> > &normal_vectors,
+                                     const bool           compute_normal_vectors,
                                      const Boundary<dim> &boundary) const;
     
                                     /**
@@ -644,6 +647,25 @@ class FiniteElement : public FiniteElementBase<dim> {
                                     const Boundary<dim>         &boundary,
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<double>      &face_jacobi_determinants) const =0;
+
+                                    /**
+                                     * Compute the normal vectors to the cell
+                                     * at the quadrature points. See the
+                                     * documentation for the #fill_fe_face_values#
+                                     * function for more details. The function
+                                     * must guarantee that the length of the
+                                     * vectors be one.
+                                     *
+                                     * Since any implementation for one
+                                     * dimension would be senseless, all
+                                     * derived classes should throw an error
+                                     * when called with #dim==1#.
+                                     */
+    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 =0;
     
                                     /**
                                      * Exception
index ab435d964e1bb1baca1a4b528dfa486162cfb735..b340857559608e505f879519c018ced731a7ee55 100644 (file)
@@ -83,6 +83,23 @@ class FELinear : public FiniteElement<dim> {
                                     const Boundary<dim>         &boundary,
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<double>      &face_jacobi_determinants) const;
+
+                                    /**
+                                     * 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;    
 };
 
 
@@ -149,6 +166,16 @@ class FEQuadratic : public FiniteElement<dim> {
                                     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_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;    
 };
 
 
@@ -215,6 +242,16 @@ class FECubic : public FiniteElement<dim> {
                                     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_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;    
 };
 
 
index c6562c858ec2a607a114c792403212f5fa767282..1a464f28783182218cbf4bde637c480226708600 100644 (file)
@@ -47,12 +47,18 @@ enum UpdateFlags {
                                        * on which the ansatz functions are
                                        * located.
                                        */
-      update_ansatz_points = 16
+      update_ansatz_points = 16,
+                                      /**
+                                       * Update the outward normal vectors
+                                       * to the face relative to this cell.
+                                       * This flag is only evaluated by
+                                       * the #FEFaceValues# class.
+                                       */
+      update_normal_vectors = 32
 };
 
 
 
-
 /*----------------------------   fe_update_flags.h     ---------------------------*/
 /* end of #ifndef __fe_update_flags_H */
 #endif
index e4946a8299b247f38a4f93458c145f3dfd22cf14..0715c95ad82a49de2dff042e9d4cbd748053cc51 100644 (file)
@@ -206,7 +206,7 @@ class FEValues {
                                      * segments, but higher order elements
                                      * may use other ways.)
                                      */
-    void reinit (const DoFHandler<dim>::cell_iterator &,
+    void reinit (const typename DoFHandler<dim>::cell_iterator &,
                 const FiniteElement<dim> &,
                 const Boundary<dim> &);
 
@@ -386,6 +386,14 @@ class FEValues {
   for the transformation of the unit face to the real face (needed to
   compute the weight factors for integration along faces). These two
   concepts have to be carefully separated.
+
+  Finally, we will often need the outward normal to a cell at the quadrature
+  points. While this could in principle be easily done using the Jacobi
+  matrices at the quadrature points and the normal vectors to the unit cell
+  (also easily derived, since they have an appealingly easy form for the unit
+  cell ;-), it is more efficiently done by the finite element class itself.
+  For example for (bi-, tri-)linear mappings the normal vector is readily
+  available without compicated matrix-vector-multiplications.
   */
 template <int dim>
 class FEFaceValues {
@@ -533,6 +541,21 @@ class FEFaceValues {
                                      * alike cells.
                                      */
     const vector<double> & get_JxW_values () const;
+
+                                    /**
+                                     * Return the outward normal vector to
+                                     * the cell at the #i#th quadrature
+                                     * point. The length of the vector
+                                     * is normalized to one.
+                                     */
+    const Point<dim> & normal_vector (const unsigned int i) const;
+    
+                                    /**
+                                     * Return the list of outward normal
+                                     * vectors to the cell at the
+                                     * quadrature points.
+                                     */
+    const vector<Point<dim> > & get_normal_vectors () const;
     
                                     /**
                                      * Reinitialize the gradients, Jacobi
@@ -553,7 +576,7 @@ class FEFaceValues {
                                      * segments, but higher order elements
                                      * may use other ways.)
                                      */
-    void reinit (const DoFHandler<dim>::cell_iterator &cell,
+    void reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                 const unsigned int                    face_no,
                 const FiniteElement<dim>             &fe,
                 const Boundary<dim>                  &boundary);
@@ -707,6 +730,13 @@ class FEFaceValues {
                                      * compute the JxW values.
                                      */
     vector<double>       face_jacobi_determinants;
+
+                                    /**
+                                     * List of outward normal vectors at the
+                                     * quadrature points. This field is filled
+                                     * in by the finite element class.
+                                     */
+    vector<Point<dim> >  normal_vectors;
     
                                     /**
                                      * Store which fields are to be updated by
@@ -833,6 +863,16 @@ FEFaceValues<dim>::get_JxW_values () const {
 
 
 
+template <int dim>
+inline
+const vector<Point<dim> > &
+FEFaceValues<dim>::get_normal_vectors () const {
+  Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+  return normal_vectors;
+};
+
+
+
 
 /*----------------------------   fe_values.h     ---------------------------*/
 /* end of #ifndef __fe_values_H */
index e08fda1ac9576a5e5fdcc98b8d2bc020866acfea..fea46fe84899bf817684055e1ad6b0956e3cb087 100644 (file)
@@ -196,6 +196,8 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &
                                            const bool               ,
                                            vector<double>          &,
                                            const bool              ,
+                                           vector<Point<1> >       &,
+                                           const bool,
                                            const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
 }
@@ -215,6 +217,8 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                                              const bool           compute_q_points,
                                              vector<double>      &face_jacobi_determinants,
                                              const bool           compute_face_jacobians,
+                                             vector<Point<dim> > &normal_vectors,
+                                             const bool           compute_normal_vectors,
                                              const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -232,17 +236,21 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
                  q_points, compute_q_points,
                  boundary);
   
-  cout << "Global unit points:\n";
-  for (unsigned int p=0; p<unit_points.size(); ++p)
-    cout << "    " << global_unit_points[p] << endl;
-
   if (compute_ansatz_points)
     get_face_ansatz_points (cell->face(face_no), boundary, ansatz_points);
 
   if (compute_face_jacobians)
     get_face_jacobians (cell->face(face_no), boundary,
                        unit_points, face_jacobi_determinants);
+
+  if (compute_normal_vectors)
+    get_normal_vectors (cell, face_no, boundary,
+                       unit_points, normal_vectors);
   
+  cout << "Global unit points:\n";
+  for (unsigned int p=0; p<unit_points.size(); ++p)
+    cout << "    " << global_unit_points[p] << endl;
+
   cout << "Global ansatz points:\n";
   for (unsigned int p=0; p<unit_points.size(); ++p)
     cout << "    " << ansatz_points[p] << endl;
index 45313af52c2abb58de76c80561510f78390c188a..739df3b6ade6fbd6952abe4216747ff6280d5f61 100644 (file)
@@ -106,6 +106,18 @@ void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+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());
+};
+
+
+
+
+
 FELinear<2>::FELinear () :
                FiniteElement<2> (1, 0, 0)
 {
@@ -357,6 +369,35 @@ void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
 
 
 
+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())));
+};
+
+
+
 
 
 
@@ -400,6 +441,18 @@ void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+void FEQuadratic<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());
+};
+
+
+
+
+
 FEQuadratic<2>::FEQuadratic () :
                FiniteElement<2> (1, 1, 1)
 {
@@ -483,6 +536,17 @@ void FEQuadratic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator
 
 
 
+template <int dim>
+void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &,
+                                          const unsigned int,
+                                          const Boundary<dim> &,
+                                          const vector<Point<dim-1> > &,
+                                          vector<Point<dim> > &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 
 
 
@@ -526,6 +590,16 @@ void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+void FECubic<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());
+};
+
+
+
 
 FECubic<2>::FECubic () :
                FiniteElement<2> (1, 2, 4) {};
@@ -598,6 +672,16 @@ void FECubic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
 
 
 
+template <int dim>
+void FECubic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &,
+                                      const unsigned int,
+                                      const Boundary<dim> &,
+                                      const vector<Point<dim-1> > &,
+                                      vector<Point<dim> > &) const {
+  Assert (false, ExcInternalError());
+};
+
+
 
 // explicit instantiations
 
index a332c178a432f47cd2e8d3d209b2dce9eba2d6c5..3fd49eb3262f9402c4faab5e8f3cd99262c944e4 100644 (file)
@@ -180,6 +180,7 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                ansatz_points (fe.total_dofs, Point<dim>()),
                jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)),
                face_jacobi_determinants (quadrature.n_quadrature_points,0),
+               normal_vectors (quadrature.n_quadrature_points,Point<dim>()),
                update_flags (update_flags),
                selected_face(0)
 {
@@ -298,6 +299,16 @@ const Point<dim> & FEFaceValues<dim>::ansatz_point (const unsigned int i) const
 
 
 
+template <int dim>
+const Point<dim> & FEFaceValues<dim>::normal_vector (const unsigned int i) const {
+  Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
+  Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+  
+  return normal_vectors[i];
+};
+
+
+
 template <int dim>
 double FEFaceValues<dim>::JxW (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
@@ -332,6 +343,8 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                            update_flags & update_q_points,
                            face_jacobi_determinants,
                            update_flags & update_JxW_values,
+                           normal_vectors,
+                           update_flags & update_normal_vectors,
                            boundary);
 
                                   // compute gradients on real element if

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.