]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add some new functions for the restriction of finite elements to subfaces. Fix a...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 5 May 1998 08:54:09 +0000 (08:54 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 5 May 1998 08:54:09 +0000 (08:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@248 0785d39b-7218-0410-832d-ea1e28bc413d

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 04cb987add6064f5eae3da10edeb57a6cc7e5a9c..602e91e364227aa1035d3d1b46d985a2baa0ac57 100644 (file)
@@ -200,7 +200,7 @@ void FiniteElement<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &
                                            const bool,
                                            const Boundary<1>       &) const {
   Assert (false, ExcNotImplemented());
-}
+};
 
 
 
@@ -250,6 +250,65 @@ void FiniteElement<dim>::fill_fe_face_values (const DoFHandler<dim>::cell_iterat
 
 
 
+
+void FiniteElement<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &,
+                                              const unsigned int       ,
+                                              const unsigned int       ,
+                                              const vector<Point<0> > &,
+                                              const vector<Point<1> > &,
+                                              vector<dFMatrix>        &,
+                                              const bool               ,
+                                              vector<Point<1> >       &,
+                                              const bool               ,
+                                              vector<double>          &,
+                                              const bool               ,
+                                              vector<Point<1> >       &,
+                                              const bool,
+                                              const Boundary<1>       &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FiniteElement<dim>::fill_fe_subface_values (const DoFHandler<dim>::cell_iterator &cell,
+                                                const unsigned int           face_no,
+                                                const unsigned int           subface_no,
+                                                const vector<Point<dim-1> > &unit_points,
+                                                const vector<Point<dim> > &global_unit_points,
+                                                vector<dFMatrix>    &jacobians,
+                                                const bool           compute_jacobians,
+                                                vector<Point<dim> > &q_points,
+                                                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()));
+  Assert (q_points.size() == unit_points.size(),
+         ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+  Assert (global_unit_points.size() == unit_points.size(),
+         ExcWrongFieldDimension(global_unit_points.size(), unit_points.size()));
+
+  static vector<Point<dim> > dummy(total_dofs);
+  fill_fe_values (cell, global_unit_points,
+                 jacobians, compute_jacobians,
+                 dummy, false,
+                 q_points, compute_q_points,
+                 boundary);
+  
+  if (compute_face_jacobians)
+    get_subface_jacobians (cell->face(face_no), subface_no,
+                          unit_points, face_jacobi_determinants);
+
+  if (compute_normal_vectors)
+    get_normal_vectors (cell, face_no, subface_no,
+                       unit_points, normal_vectors);
+};
+
+
 /*------------------------------- Explicit Instantiations -------------*/
 
 template class FiniteElementData<1>;
index 739df3b6ade6fbd6952abe4216747ff6280d5f61..e242af02055404a205ad6f47f25b12bb153b9222 100644 (file)
@@ -106,6 +106,15 @@ void FELinear<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+void FELinear<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+                                        const unsigned int           ,
+                                        const vector<Point<0> > &,
+                                        vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
                                      const unsigned int,
                                      const Boundary<1> &,
@@ -116,6 +125,16 @@ void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 
 
 
+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());
+};
+
+
+
 
 
 FELinear<2>::FELinear () :
@@ -369,6 +388,27 @@ void FELinear<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face,
 
 
 
+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);
+};
+
+
+
 void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
                                      const unsigned int       face_no,
                                      const Boundary<2> &,
@@ -398,6 +438,40 @@ void FELinear<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
 
 
 
+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())));
+};
+
+
+
 
 
 
@@ -441,6 +515,15 @@ void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+void FEQuadratic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+                                           const unsigned int           ,
+                                           const vector<Point<0> > &,
+                                           vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
                                         const unsigned int,
                                         const Boundary<1> &,
@@ -451,6 +534,15 @@ void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 
 
 
+void FEQuadratic<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());
+};
+
+
 
 
 FEQuadratic<2>::FEQuadratic () :
@@ -528,9 +620,22 @@ void FEQuadratic<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::f
 
 template <int dim>
 void FEQuadratic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
-                                     const Boundary<dim>         &,
-                                     const vector<Point<dim-1> > &,
-                                     vector<double>      &) const {
+                                          const Boundary<dim>         &,
+                                          const vector<Point<dim-1> > &,
+                                          vector<double>      &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                             const unsigned int           ,
+                                             const vector<Point<dim-1> > &,
+                                             vector<double>      &) const {
+  Assert (face->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+
   Assert (false, ExcNotImplemented());
 };
 
@@ -542,7 +647,20 @@ void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator
                                           const Boundary<dim> &,
                                           const vector<Point<dim-1> > &,
                                           vector<Point<dim> > &) const {
-  Assert (false, ExcInternalError());
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                          const unsigned int           face_no,
+                                          const unsigned int,
+                                          const vector<Point<dim-1> > &,
+                                          vector<Point<dim> > &) const {
+  Assert (cell->face(face_no)->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+  Assert (false, ExcNotImplemented());
 };
 
 
@@ -590,6 +708,15 @@ void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
 
 
 
+void FECubic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+                                       const unsigned int           ,
+                                       const vector<Point<0> > &,
+                                       vector<double>      &) const {
+  Assert (false, ExcInternalError());
+};
+
+
+
 void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
                                     const unsigned int,
                                     const Boundary<1> &,
@@ -600,6 +727,15 @@ void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 
 
 
+void FECubic<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());
+};
+
+
 
 FECubic<2>::FECubic () :
                FiniteElement<2> (1, 2, 4) {};
@@ -672,13 +808,40 @@ void FECubic<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &,
 
 
 
+template <int dim>
+void FECubic<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                         const unsigned int           ,
+                                         const vector<Point<dim-1> > &,
+                                         vector<double>      &) const {
+  Assert (face->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+
+  Assert (false, ExcNotImplemented());
+};
+
+
+
 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());
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+void FECubic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                      const unsigned int           face_no,
+                                      const unsigned int           ,
+                                      const vector<Point<dim-1> > &,
+                                      vector<Point<dim> > &) const {
+  Assert (cell->face(face_no)->at_boundary() == false,
+         ExcBoundaryFaceUsed ());
+
+  Assert (false, ExcNotImplemented());
 };
 
 
index 462690cf5cf8c906e23d3e97b347e7cbfb25cfbb..eca2ef337d8f87b4c0ba8e8ef1c2ef48685c41eb 100644 (file)
@@ -160,7 +160,7 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                                     vector<Point<dim> >(quadrature.n_quadrature_points)),
                unit_quadrature_points(quadrature.get_quad_points())
 {
-  Assert ((update_flags | update_normal_vectors) == false,
+  Assert ((update_flags & update_normal_vectors) == false,
          ExcInvalidUpdateFlag());
 
   for (unsigned int i=0; i<fe.total_dofs; ++i)
@@ -254,20 +254,15 @@ FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
                                   n_dofs,
                                   n_faces_or_subfaces,
                                   update_flags),
+               unit_shape_gradients (n_faces_or_subfaces,
+                                     vector<vector<Point<dim> > >(n_dofs,
+                                                                  vector<Point<dim> >(n_q_points))),
                unit_face_quadrature_points (n_q_points, Point<dim-1>()),
                unit_quadrature_points (n_faces_or_subfaces,
                                        vector<Point<dim> >(n_q_points, Point<dim>())),
                face_jacobi_determinants (n_q_points, 0),
                normal_vectors (n_q_points)
-{
-  for (unsigned int i=0; i<n_faces_or_subfaces; ++i) 
-    {
-      unit_shape_gradients[i].resize (n_dofs,
-                                     vector<Point<dim> >(n_q_points));
-      unit_quadrature_points[i].resize (n_q_points,
-                                       Point<dim>());
-    };
-};
+{};
 
 
 
@@ -358,10 +353,11 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
       Assert (update_flags & update_jacobians, ExcCannotInitializeField());
       
       for (unsigned int i=0; i<fe.total_dofs; ++i)
-       for (unsigned int j=0; j<n_quadrature_points; ++j) 
-         {
-           shape_gradients[i][j] = Point<dim>();
-           
+       {
+         fill_n (shape_gradients[i].begin(),
+                 n_quadrature_points,
+                 Point<dim>());
+         for (unsigned int j=0; j<n_quadrature_points; ++j) 
            for (unsigned int s=0; s<dim; ++s)
                                               // (grad psi)_s =
                                               // (grad_{\xi\eta})_b J_{bs}
@@ -370,7 +366,7 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                shape_gradients[i][j](s)
                  += (unit_shape_gradients[face_no][i][j](b) *
                      jacobi_matrices[j](b,s));
-         };
+       };
     };
   
   
@@ -443,6 +439,9 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
                                   const unsigned int         subface_no,
                                   const FiniteElement<dim>  &fe,
                                   const Boundary<dim>       &boundary) {
+  Assert (cell->face(face_no)->at_boundary() == false,
+         ExcReinitCalledWithBoundaryFace());
+  
   present_cell  = cell;
   selected_dataset = face_no*(1<<(dim-1)) + subface_no;
                                   // fill jacobi matrices and real
@@ -471,11 +470,12 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
     {
       Assert (update_flags & update_jacobians, ExcCannotInitializeField());
       
-      for (unsigned int i=0; i<fe.total_dofs; ++i)
-       for (unsigned int j=0; j<n_quadrature_points; ++j) 
-         {
-           shape_gradients[i][j] = Point<dim>();
-           
+      for (unsigned int i=0; i<fe.total_dofs; ++i) 
+       {
+         fill_n (shape_gradients[i].begin(),
+                 n_quadrature_points,
+                 Point<dim>());
+         for (unsigned int j=0; j<n_quadrature_points; ++j) 
            for (unsigned int s=0; s<dim; ++s)
                                               // (grad psi)_s =
                                               // (grad_{\xi\eta})_b J_{bs}
@@ -484,7 +484,7 @@ void FESubfaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator
                shape_gradients[i][j](s)
                  += (unit_shape_gradients[selected_dataset][i][j](b) *
                      jacobi_matrices[j](b,s));
-         };
+       };
     };
   
   

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.