]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Start of implementation of (bi-)quadratic elements.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 25 Jun 1998 11:20:23 +0000 (11:20 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 25 Jun 1998 11:20:23 +0000 (11:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@411 0785d39b-7218-0410-832d-ea1e28bc413d

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

index dda7a467590710cf156ac40188a91e34906b539e..1b3f40346319f9fd5d28da0b6b742190a8eee71e 100644 (file)
@@ -165,16 +165,16 @@ class FELinear : public FiniteElement<dim> {
 
 /**
  * Define a (bi-, tri-, etc)quadratic finite element in #dim# space dimensions.
- * In one space dimension, a linear (subparametric) mapping from the unit cell
+ * A linear (subparametric) mapping from the unit cell
  * to the real cell is implemented.
  */
 template <int dim>
-class FEQuadratic : public FiniteElement<dim> {
+class FEQuadraticSub : public FiniteElement<dim> {
   public:
                                     /**
                                      * Constructor
                                      */
-    FEQuadratic ();
+    FEQuadraticSub ();
 
                                     /**
                                      * Return the value of the #i#th shape
index 0dc692538b0e74a9b639f7aa766794be42c53b16..cfdf0bb20454fa99f5e510c064d91fad3fc32a8a 100644 (file)
@@ -47,7 +47,7 @@ FELinear<1>::FELinear () :
 template <>
 double
 FELinear<1>::shape_value(const unsigned int i,
-                        const Point<1>p) const
+                        const Point<1>     &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   switch (i)
@@ -647,13 +647,13 @@ void FELinear<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face
 #if deal_II_dimension == 1
 
 template <>
-FEQuadratic<1>::FEQuadratic () :
+FEQuadraticSub<1>::FEQuadraticSub () :
                FiniteElement<1> (1, 1) {};
 
 
 
 template <>
-void FEQuadratic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
+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,
@@ -672,16 +672,53 @@ void FEQuadratic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
 
 
 template <>
-void FEQuadratic<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
-                                       const Boundary<1>  &boundary,
-                                       vector<Point<1> >  &ansatz_points) const {
+double
+FEQuadraticSub<1>::shape_value(const unsigned int i,
+                              const Point<1>     &p) const
+{
+  Assert((i<total_dofs), ExcInvalidIndex(i));
+  const double xi = p(0);
+  switch (i)
+    {
+      case 0: return (1-xi)*(1-2*xi);
+      case 1: return xi*(2*xi-1);
+      case 2: return 4*xi*(1-xi);
+    }
+  return 0.;
+}
+
+
+
+template <>
+inline
+Point<1>
+FEQuadraticSub<1>::shape_grad(const unsigned int i,
+                             const Point<1>    &p) const
+{
+  Assert((i<total_dofs), ExcInvalidIndex(i));
+  const double xi = p(0);
+  switch (i)
+    {
+      case 0: return Point<1>(-3+4*xi);
+      case 1: return Point<1>(4*xi-1);
+      case 2: return Point<1>(4-8*xi);
+    }
+  return Point<1>();
+};
+
+
+
+template <>
+void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell,
+                                          const Boundary<1>  &boundary,
+                                          vector<Point<1> >  &ansatz_points) const {
   FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points);
 };
 
 
 
 template <>
-void FEQuadratic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
+void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &,
                                             const Boundary<1>  &,
                                             vector<Point<1> >  &) const {
   Assert (false, ExcInternalError());
@@ -690,42 +727,42 @@ void FEQuadratic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_
 
 
 template <>
-void FEQuadratic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                        const Boundary<1>         &,
-                                        const vector<Point<0> > &,
-                                        vector<double>      &) const {
+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 FEQuadratic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                           const unsigned int           ,
-                                           const vector<Point<0> > &,
-                                           vector<double>      &) const {
+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 FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                        const unsigned int,
-                                        const Boundary<1> &,
-                                        const vector<Point<0> > &,
-                                        vector<Point<1> > &) const {
+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 FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                        const unsigned int,
-                                        const unsigned int,
-                                        const vector<Point<0> > &,
-                                        vector<Point<1> > &) const {
+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());
 };
 
@@ -735,7 +772,7 @@ void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 #if deal_II_dimension == 2
 
 template <>
-FEQuadratic<2>::FEQuadratic () :
+FEQuadraticSub<2>::FEQuadraticSub () :
                FiniteElement<2> (1, 1, 1)
 {
   interface_constraints(0,2) = 1.0;
@@ -751,123 +788,257 @@ FEQuadratic<2>::FEQuadratic () :
   Assert (false, ExcNotImplemented());
 };
 
-#endif
-
-
 
-template <int dim>
+template <>
 double
-FEQuadratic<dim>::shape_value (const unsigned int i,
-                              const Point<dim> &) const
+FEQuadraticSub<2>::shape_value (const unsigned int i,
+                               const Point<2>    &p) const
 {
-  Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, ExcNotImplemented());
-  return 0.;
+  Assert (i<total_dofs, ExcInvalidIndex(i));
+
+  const double xi = p(0),
+              eta= p(1);
+  switch (i)
+    {
+      case 0: return (1.-xi)*( 2*xi-1) * (1.-eta)*( 2*eta-1);
+      case 1: return     xi *(-2*xi+1) * (1.-eta)*( 2*eta-1);
+      case 2: return     xi *(-2*xi+1) *     eta *(-2*eta+1);
+      case 3: return (1.-xi)*( 2*xi-1) *     eta *(-2*eta+1);
+      case 4: return 4 * (1-xi)*xi * (1-eta)*(1-2*eta);
+      case 5: return 4 *    xi *(-1+2*xi) * (1-eta)*eta;
+      case 6: return 4 * (1-xi)*xi *    eta *(-1+2*eta);
+      case 7: return 4 * (1.-xi)*(1-2*xi) * (1-eta)*eta;
+      case 8: return 16 * xi*(1-xi) * eta*(1-eta);
+    };
+  return 0;
 };
 
 
 
-template <int dim>
-Point<dim>
-FEQuadratic<dim>::shape_grad (const unsigned int i,
-                             const Point<dim> &) const
+template <>
+Point<2>
+FEQuadraticSub<2>::shape_grad (const unsigned int i,
+                              const Point<2>    &p) const
 {
-  Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
-  Assert (false, ExcNotImplemented());
-  return Point<dim> ();
-};
+  Assert (i<total_dofs, ExcInvalidIndex(i));
 
+  const double xi = p(0),
+              eta= p(1);
+  switch (i)
+    {
+      case 0: return Point<2>((-4*xi+3) * (1.-eta)*( 2*eta-1),
+                             (1.-xi)*( 2*xi-1) * (-4*eta+3));
+      case 1: return Point<2>((-4*xi+1) * (1.-eta)*( 2*eta-1) ,
+                             xi *(-2*xi+1) * (-4*eta+3));
+      case 2: return Point<2>((-4*xi+1) *     eta *(-2*eta+1),
+                             xi *(-2*xi+1) * (-4*eta+1));
+      case 3: return Point<2>((-4*xi+3) *     eta *(-2*eta+1),
+                             (1.-xi)*( 2*xi-1) * (-4*eta+1));
+      case 4: return Point<2>(4 * (1-2*xi) * (1-eta)*(1-2*eta),
+                             4 * (1-xi)*xi * (4*eta-3));
+      case 5: return Point<2>(4 * (4*xi-1) * (1-eta)*eta,
+                             4 *    xi *(-1+2*xi) * (1-2*eta));
+      case 6: return Point<2>(4 * (1-2*xi) *    eta *(-1+2*eta),
+                             4 * (1-xi)*xi * (4*eta-1));
+      case 7: return Point<2>(4 * (4*xi-3) * (1-eta)*eta,
+                             4 * (1.-xi)*(1-2*xi) * (1-2*eta));
+      case 8: return Point<2>(16 * (1-2*xi) * eta*(1-eta),
+                             16 * xi*(1-xi) * (1-2*eta));
+    };
+  return Point<2> ();
+};
 
 
-template <int dim>
-void FEQuadratic<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &,
-                                      const vector<Point<dim> >            &unit_points,
-                                      vector<dFMatrix>  &jacobians,
-                                      const bool,
-                                      vector<Point<dim> > &ansatz_points,
-                                      const bool,
-                                      vector<Point<dim> > &q_points,
-                                      const bool,
-                                      const Boundary<dim> &) 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()));
+template <>
+void FEQuadraticSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell,
+                                          const Boundary<2>&,
+                                          vector<Point<2> >  &ansatz_points) const {
   Assert (ansatz_points.size() == total_dofs,
-         ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+         ExcWrongFieldDimension (ansatz_points.size(), total_dofs));
+  
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    ansatz_points[vertex] = cell->vertex(vertex);
 
-  Assert (false, ExcNotImplemented());
+                                  // for the bilinear mapping, the centers
+                                  // of the face on the unit cell are mapped
+                                  // to the mean coordinates of the vertices
+  for (unsigned int line=0; line<4; ++line)
+    ansatz_points[4+line] = (cell->line(line)->vertex(0) +
+                            cell->line(line)->vertex(1)) / 2;
+                                  // same for the center of the square:
+                                  // since all four linear basis functions
+                                  // take on the value 1/4 at the center,
+                                  // the center is mapped to the mean
+                                  // coordinates of the four vertices
+  ansatz_points[8] = (ansatz_points[0] +
+                     ansatz_points[1] +
+                     ansatz_points[2] +
+                     ansatz_points[3]) / 4;
 };
 
 
 
-template <int dim>
-void FEQuadratic<dim>::get_ansatz_points (const typename DoFHandler<dim>::cell_iterator &,
-                                         const Boundary<dim>  &,
-                                         vector<Point<dim> >  &) const {
-  Assert (false, ExcNotImplemented());
-};
-
-
+template <>
+void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face,
+                                               const Boundary<2>  &,
+                                               vector<Point<2> >  &ansatz_points) const {
+  Assert (ansatz_points.size() == dofs_per_face,
+         ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face));
 
-template <int dim>
-void FEQuadratic<dim>::get_face_ansatz_points (const typename DoFHandler<dim>::face_iterator &,
-                                              const Boundary<dim>  &,
-                                              vector<Point<dim> >  &) const {
-  Assert (false, ExcNotImplemented());
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    ansatz_points[vertex] = face->vertex(vertex);
+  ansatz_points[2] = (ansatz_points[0] + ansatz_points[1]) / 1;
 };
 
 
 
-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 {
-  Assert (false, ExcNotImplemented());
+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 <int dim>
-void FEQuadratic<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+template <>
+void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face,
                                              const unsigned int           ,
-                                             const vector<Point<dim-1> > &,
-                                             vector<double>      &) const {
+                                             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 ());
 
-  Assert (false, ExcNotImplemented());
+                                  // 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 <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, ExcNotImplemented());
+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 <int dim>
-void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
-                                          const unsigned int           face_no,
+template <>
+void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell,
+                                          const unsigned int       face_no,
                                           const unsigned int,
-                                          const vector<Point<dim-1> > &,
-                                          vector<Point<dim> > &) const {
+                                          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 &,
+                                      const vector<Point<dim> >            &unit_points,
+                                      vector<dFMatrix>  &jacobians,
+                                      const bool,
+                                      vector<Point<dim> > &ansatz_points,
+                                      const bool,
+                                      vector<Point<dim> > &q_points,
+                                      const bool,
+                                      const Boundary<dim> &) 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));
+
   Assert (false, ExcNotImplemented());
 };
 
 
 
+
+
+
+
 template <int dim>
-void FEQuadratic<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
+void FEQuadraticSub<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
                                              const Boundary<dim> &,
                                              dFMatrix &) const {
   Assert (false, ExcNotImplemented());
@@ -1102,6 +1273,6 @@ void FECubic<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &
 // explicit instantiations
 
 template class FELinear<deal_II_dimension>;
-template class FEQuadratic<deal_II_dimension>;
+template class FEQuadraticSub<deal_II_dimension>;
 template class FECubic<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.