#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>
/**
- * 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>
/**
- * 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>
/**
- * 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>
--- /dev/null
+//---------------------------- $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>;