]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
continuing on FEDG_Pk
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jun 2000 22:11:44 +0000 (22:11 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jun 2000 22:11:44 +0000 (22:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@3003 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_lib.dgp.h
deal.II/deal.II/source/fe/fe_lib.dgp2.cc
deal.II/deal.II/source/fe/fe_lib.dgp3.cc [new file with mode: 0644]

index 2e407a24fa09a064034f6a9d3d03e25ee80a71cc..cc3f4a637f0a19661a7698a351b0bd317a84eab4 100644 (file)
 
 #include <fe/q1_mapping.h>
 
+
 /**
- * Isoparametric Q1 finite element in #dim# space dimensions.
- *
- * The linear, isoparametric mapping from a point $\vec \xi$ on the unit cell
- * to a point $\vec x$ on the real cell is defined as
- * $$ \vec x(\vec \xi)  = \sum_j {\vec p_j} N_j(\xi) $$
- * where $\vec p_j$ is the vector to the $j$th corner point of the cell in
- * real space and $N_j(\vec \xi)$ is the value of the basis function associated
- * with the $j$th corner point, on the unit cell at point $\vec \xi$. The sum
- * over $j$ runs over all corner points.
+ * Discontinuous P1-element on hypercubes.
  *
- * The number of degrees of freedom equal the number of the respective vertex
- * of the cell
+ * This is the implementation of a linear (sic) polynomial space on a
+ * d-dimensional hypercube. The shape functions are the first @p{d+1}
+ * of @p{1,x,y,z}. Later on, these should be exchanged for mutually
+ * orthogonal, preferably by changing the unit cell to $[-1,1]^d$.
  *
- * @author Wolfgang Bangerth, 1998, 1999
+ * @author Guido Kanschat, 2000
  */
 template <int dim>
 class FEDG_P1 : public FEQ1Mapping<dim>
@@ -119,86 +114,15 @@ class FEDG_P1 : public FEQ1Mapping<dim>
 
 
 /**
- * Subparametric Q2 finite element in #dim# space dimensions.
- * A Q1 mapping from the unit cell
- * to the real cell is implemented.
- *
- * The numbering of the degrees of freedom is as follows:
- * \begin{itemize}
- * \item 1D case:
- *   \begin{verbatim}
- *      0---2---1
- *   \end{verbatim}
- *
- * \item 2D case:
- *   \begin{verbatim}
- *      3---6---2
- *      |       |
- *      7   8   5
- *      |       |
- *      0---4---1
- *   \end{verbatim}
+ * Discontinuous P2-element on hypercubes.
  *
- * \item 3D case:
- *   \begin{verbatim}
- *         7--14---6        7--14---6
- *        /|       |       /       /|
- *      19 |       13     19      1813
- *      /  15      |     /       /  |
- *     3   |       |    3---10--2   |
- *     |   4--12---5    |       |   5
- *     |  /       /     |       9  /
- *    11 16      17     11      | 17
- *     |/       /       |       |/
- *     0---8---1        0---8---1
+ * This is the implementation of a linear (sic) polynomial space on a
+ * d-dimensional hypercube. The shape functions are those of
+ * @p{1,x,y,z, x*x, x*y, x*z, y*y, y*z, z*z} applying to the space
+ * dimension. Later on, these should be exchanged for mutually
+ * orthogonal, preferably by changing the unit cell to $[-1,1]^d$.
  *
- *         *-------*        *-------*
- *        /|       |       /       /|
- *       / |  21   |      /  24   / |
- *      /  |       |     /       /  |
- *     *   |       |    *-------*   |
- *     |25 *-------*    |       |23 *
- *     |  /       /     |   20  |  /
- *     | /  22   /      |       | /
- *     |/       /       |       |/
- *     *-------*        *-------* 
- *   \end{verbatim}
- *   The center vertex has number 26.
- *
- *   The respective coordinate values of the support points of the degrees
- *   of freedom are as follows:
- *   \begin{itemize}
- *   \item Index 0: #[0, 0, 0]#;
- *   \item Index 1: #[1, 0, 0]#;
- *   \item Index 2: #[1, 0, 1]#;
- *   \item Index 3: #[0, 0, 1]#;
- *   \item Index 4: #[0, 1, 0]#;
- *   \item Index 5: #[1, 1, 0]#;
- *   \item Index 6: #[1, 1, 1]#;
- *   \item Index 7: #[0, 1, 1]#;
- *   \item Index 8: #[1/2, 0, 0]#;
- *   \item Index 9: #[1, 0, 1/2]#;
- *   \item Index 10: # [1/2, 0, 1]#;
- *   \item Index 11: # [0, 0, 1/2]#;
- *   \item Index 12: # [1/2, 1, 0]#;
- *   \item Index 13: # [1, 1, 1/2]#;
- *   \item Index 14: # [1/2, 1, 1]#;
- *   \item Index 15: # [0, 1, 1/2]#;
- *   \item Index 16: # [0, 1/2, 0]#;
- *   \item Index 17: # [1, 1/2, 0]#;
- *   \item Index 18: # [1, 1/2, 1]#;
- *   \item Index 19: # [0, 1/2, 1]#;
- *   \item Index 20: # [1/2, 0, 1/2]#;
- *   \item Index 21: # [1/2, 1, 1/2]#;
- *   \item Index 22: # [1/2, 1/2, 0]#;
- *   \item Index 23: # [1, 1/2, 1/2]#;
- *   \item Index 24: # [1/2, 1/2, 1]#;
- *   \item Index 25: # [0, 1/2, 1/2]#;
- *   \item Index 26: # [1/2, 1/2, 1/2]#; 
- *   \end{itemize}
- * \end{itemize}
- *
- * @author Wolfgang Bangerth, 1998, 1999
+ * @author Guido Kanschat, 2000
  */
 template <int dim>
 class FEDG_P2 : public FEQ1Mapping<dim>
@@ -282,30 +206,15 @@ class FEDG_P2 : public FEQ1Mapping<dim>
 
 
 /**
- * Subparametric Q3 finite element in #dim# space dimensions.
- * A Q1 mapping from the unit cell
- * to the real cell is implemented.
+ * Discontinuous P3-element on hypercubes.
  *
- * The numbering of degrees of freedom in one spatial dimension is as follows:
- * \begin{verbatim}
- *   0--2--3--1
- * \end{verbatim}
+ * This is the implementation of a linear (sic) polynomial space on a
+ * d-dimensional hypercube. The shape functions are the basis
+ * polynomials spanning the space of cubic polynomials. Later on,
+ * they should be exchanged for mutually orthogonal, preferably by
+ * changing the unit cell to $[-1,1]^d$.
  *
- * The numbering of degrees of freedom in two spatial dimension is as follows:
- * \begin{verbatim}
- *   3--8--9--2
- *   |        |
- *   11 15 14 7
- *   |        |
- *   10 12 13 6
- *   |        |
- *   0--4--5--1
- * \end{verbatim}
- * Note the reverse ordering of degrees of freedom on the left and upper
- * line and the counterclockwise numbering of the interior degrees of
- * freedom.
- *
- * @author Wolfgang Bangerth, 1998
+ * @author Guido Kanschat, 2000
  */
 template <int dim>
 class FEDG_P3 : public FEQ1Mapping<dim>
@@ -390,32 +299,15 @@ class FEDG_P3 : public FEQ1Mapping<dim>
 
 
 /**
- * Subparametric Q4 finite element in #dim# space dimensions.
- * A linear (subparametric) mapping from the unit cell
- * to the real cell is implemented.
- *
- * The numbering of degrees of freedom in one spatial dimension is as follows:
- * \begin{verbatim}
- *   0--2--3--4--1
- * \end{verbatim}
+ * Discontinuous P3-element on hypercubes.
  *
- * The numbering of degrees of freedom in two spatial dimension is as follows:
- * \begin{verbatim}
- *   3--10-11-12-2
- *   |           |
- *   15 19 22 18 9
- *   |           |
- *   14 23 24 21 8
- *   |           |
- *   13 16 20 17 7
- *   |           |
- *   0--4--5--6--1
- * \end{verbatim}
- * Note the reverse ordering of degrees of freedom on the left and upper
- * line and the numbering of the interior degrees of
- * freedom.
+ * This is the implementation of a linear (sic) polynomial space on a
+ * d-dimensional hypercube. The shape functions are the basis
+ * polynomials spanning the space of cubic polynomials. Later on,
+ * they should be exchanged for mutually orthogonal, preferably by
+ * changing the unit cell to $[-1,1]^d$.
  *
- * @author Wolfgang Bangerth, 1998
+ * @author Guido Kanschat, 2000
  */
 template <int dim>
 class FEDG_P4 : public FEQ1Mapping<dim>
index 3168fc4d5727d624f4996a04fa7c7e1989978de3..555f5fd5c723d472d4287c44066af6f3fc48c0c6 100644 (file)
@@ -237,6 +237,7 @@ FEDG_P2<3>::FEDG_P2 () :
                FEQ1Mapping<3> (0, 0, 0, 4, 1,
                                vector<bool> (1, true))
 {
+  Assert(false, ExcNotImplemented ());
 //  initialize_matrices ();
 };
 
diff --git a/deal.II/deal.II/source/fe/fe_lib.dgp3.cc b/deal.II/deal.II/source/fe/fe_lib.dgp3.cc
new file mode 100644 (file)
index 0000000..d4ca025
--- /dev/null
@@ -0,0 +1,375 @@
+//----------------------------  $RCSFile$  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  $RCSFile$  ---------------------------
+
+
+#include <fe/fe_lib.dgp.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <grid/geometry_info.h>
+#include <algorithm>
+
+
+// declare explicit specializations before use:
+template <> void FEDG_P3<deal_II_dimension>::initialize_matrices ();
+
+
+#if deal_II_dimension == 1
+
+template <>
+FEDG_P3<1>::FEDG_P3 () :
+               FEQ1Mapping<1> (0, 4, 0, 0, 1,
+                               vector<bool> (1, true))
+{
+//  initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P3<1>::initialize_matrices ()
+{
+  Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+double
+FEDG_P3<1>::shape_value(const unsigned int i,
+                        const Point<1>     &p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+  switch (i)
+    {
+      case 0: return 1.;
+      case 1: return p(0);
+      case 2: return p(0)*p(0);
+      case 2: return p(0)*p(0)*p(0);
+    }
+  return 0.;
+}
+
+
+template <>
+inline
+Tensor<1,1>
+FEDG_P3<1>::shape_grad(const unsigned int i,
+                       const Point<1>&p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
+  switch (i)
+    {
+      case 0: return Point<1>(-1.);
+      case 1: return Point<1>(1.);
+      case 2: return Point<1>(2.*p(0));
+      case 2: return Point<1>(3.*p(0)*p(0));
+           
+    }
+  return Point<1>();
+};
+
+
+template <>
+inline
+Tensor<2,1>
+FEDG_P3<1>::shape_grad_grad (const unsigned int i,
+                            const Point<1> &) const
+{
+  Assert(false, ExcNotImplemented());
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+  return Tensor<2,1>();
+};
+
+
+template <>
+void FEDG_P3<1>::get_unit_support_points (vector<Point<1> >  &support_points) const
+{
+  FiniteElement<1>::get_unit_support_points (support_points);
+};
+
+
+template <>
+void FEDG_P3<1>::get_support_points (const DoFHandler<1>::cell_iterator &cell,
+                                    vector<Point<1> >  &support_points) const
+{
+  FiniteElement<1>::get_support_points (cell, support_points);
+};
+
+
+template <>
+void FEDG_P3<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
+                                         vector<Point<1> >  &) const
+{
+  Assert (false, ExcInternalError());
+};
+
+
+template <>
+void FEDG_P3<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
+                                        FullMatrix<double> &local_mass_matrix) const
+{
+  Assert(false, ExcNotImplemented());
+  Assert (local_mass_matrix.n() == dofs_per_cell,
+         ExcWrongFieldDimension(local_mass_matrix.n(),dofs_per_cell));
+  Assert (local_mass_matrix.m() == dofs_per_cell,
+         ExcWrongFieldDimension(local_mass_matrix.m(),dofs_per_cell));
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+FEDG_P3<2>::FEDG_P3 () :
+               FEQ1Mapping<2> (0, 0, 10, 0, 1,
+                               vector<bool> (1, true))
+{
+//  initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P3<2>::initialize_matrices ()
+{
+  Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P3<2>::shape_value (const unsigned int i,
+                         const Point<2>& p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+  switch (i)
+    {
+      case 0: return 1;
+      case 1: return p(0);
+      case 2: return p(1);
+      case 3: return p(0)*p(0);
+      case 4: return p(0)*p(1);
+      case 5: return p(1)*p(1);
+      case 6: return p(0)*p(0)*p(0);
+      case 7: return p(0)*p(0)*p(1);
+      case 8: return p(0)*p(1)*p(1);
+      case 9: return p(1)*p(1)*p(1);
+    }
+  return 0.;
+};
+
+
+template <>
+inline
+Tensor<1,2>
+FEDG_P3<2>::shape_grad (const unsigned int i,
+                        const Point<2>& p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
+  switch (i)
+    {
+      case 0: return Point<2> (0,0);
+      case 1: return Point<2> (1,0);
+      case 2: return Point<2> (0,1);
+      case 3: return Point<2> (2*p(0),0);
+      case 4: return Point<2> (p(1),p(0));
+      case 5: return Point<2> (0,2*p(1));
+      case 6: return Point<2> (3*p(0)*p(0), 0);
+      case 7: return Point<2> (2*p(0)*p(1), p(0)*p(0));
+      case 8: return Point<2> (p(1)*p(1),   2*p(0)*p(1));
+      case 9: return Point<2> (0,           3*p(1)*p(1));
+    }
+  return Point<2> ();
+};
+
+
+template <>
+inline
+Tensor<2,2>
+FEDG_P3<2>::shape_grad_grad (const unsigned int i,
+                             const Point<2> &) const
+{
+  Assert(false, ExcNotImplemented());
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+  return Tensor<2,2>();
+};
+
+
+template <>
+void FEDG_P3<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &,
+                                        FullMatrix<double> &) const
+{
+  Assert(false, ExcNotImplemented ());
+};
+
+
+template <>
+void FEDG_P3<2>::get_unit_support_points (vector<Point<2> > &unit_points) const
+{
+  Assert(false, ExcNotImplemented ());
+  Assert (unit_points.size() == dofs_per_cell,
+         ExcWrongFieldDimension (unit_points.size(), dofs_per_cell));
+
+  unit_points[0] = Point<2> (.5,.5);
+  unit_points[1] = Point<2> (1,0);
+  unit_points[2] = Point<2> (0,1);
+  unit_points[3] = Point<2> (1,0);
+  unit_points[4] = Point<2> (0,1);
+  unit_points[5] = Point<2> (1,1);
+};
+
+
+#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+FEDG_P3<3>::FEDG_P3 () :
+               FEQ1Mapping<3> (0, 0, 0, 4, 1,
+                               vector<bool> (1, true))
+{
+  Assert(false, ExcNotImplemented ());
+//  initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P3<3>::initialize_matrices ()
+{
+  Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P3<3>::shape_value (const unsigned int i,
+                        const Point<3>& p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+  switch (i)
+    {
+      case 0: return 1.;
+      case 1: return p(0);
+      case 2: return p(1);
+      case 3: return p(2);
+    }
+  return 0.;
+};
+
+
+template <>
+inline
+Tensor<1,3>
+FEDG_P3<3>::shape_grad (const unsigned int i,
+                        const Point<3>& p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+                                  // originally, the return type of the
+                                  // function was Point<dim>, so we
+                                  // still construct it as that. it should
+                                  // make no difference in practice,
+                                  // however
+  switch (i)
+    {
+      case 0: return Point<3>(0,0,0);
+      case 1: return Point<3>(1,0,0);
+      case 2: return Point<3>(0,1,0);
+      case 3: return Point<3>(0,0,1);
+    }
+  return Point<3> ();
+};
+
+
+template <>
+inline
+Tensor<2,3>
+FEDG_P3<3>::shape_grad_grad (const unsigned int i,
+                             const Point<3> &p) const
+{
+  Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+
+  Tensor<2,3> return_value;
+  return return_value;
+};
+
+
+template <>
+void FEDG_P3<3>::get_local_mass_matrix (const DoFHandler<3>::cell_iterator &,
+                                        FullMatrix<double> &local_mass_matrix) const
+{
+  Assert (local_mass_matrix.n() == dofs_per_cell,
+         ExcWrongFieldDimension(local_mass_matrix.n(),dofs_per_cell));
+  Assert (local_mass_matrix.m() == dofs_per_cell,
+         ExcWrongFieldDimension(local_mass_matrix.m(),dofs_per_cell));
+
+  AssertThrow (false, ExcComputationNotUseful(3));
+};
+
+
+template <>
+void FEDG_P3<3>::get_unit_support_points (vector<Point<3> > &unit_points) const {
+  Assert (unit_points.size() == dofs_per_cell,
+         ExcWrongFieldDimension (unit_points.size(), dofs_per_cell));
+
+  unit_points[0] = Point<3> (.5,.5,.5);
+  unit_points[1] = Point<3> (1,0,0);
+  unit_points[2] = Point<3> (0,1,0);
+  unit_points[3] = Point<3> (0,0,1);
+};
+
+
+#endif
+
+
+template <int dim>
+void
+FEDG_P3<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
+                                  vector<Point<dim> >  &support_points) const
+{
+  Assert (support_points.size() == dofs_per_cell,
+         ExcWrongFieldDimension (support_points.size(), dofs_per_cell));
+  
+  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
+    support_points[vertex] = cell->vertex(vertex);
+};
+
+
+template <int dim>
+void
+FEDG_P3<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
+                                      vector<Point<dim> >  &support_points) const
+{
+  Assert ((support_points.size() == dofs_per_face) &&
+         (support_points.size() == GeometryInfo<dim>::vertices_per_face),
+         ExcWrongFieldDimension (support_points.size(),
+                                 GeometryInfo<dim>::vertices_per_face));
+
+  for (unsigned int vertex=0; vertex<dofs_per_face; ++vertex)
+    support_points[vertex] = face->vertex(vertex);
+};
+
+
+// explicit instantiations
+
+template class FEDG_P3<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.