--- /dev/null
+//---------------------------- fe_lib.dgp.h ---------------------------
+// $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.
+//
+//---------------------------- fe_lib.dgp.h ---------------------------
+#ifndef __deal2__fe_lib_dgp_h
+#define __deal2__fe_lib_dgp_h
+
+
+/*---------------------------- fe_lib.h ---------------------------*/
+
+
+#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.
+ *
+ * The number of degrees of freedom equal the number of the respective vertex
+ * of the cell
+ *
+ * @author Wolfgang Bangerth, 1998, 1999
+ */
+template <int dim>
+class FEDG_P1 : public FEQ1Mapping<dim>
+{
+ public:
+ /**
+ * Constructor
+ */
+ FEDG_P1 ();
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the tensor of second derivatives
+ * of the #i#th shape function at
+ * point #p# on the unit cell.
+ *
+ * For linear elements, all second
+ * derivatives on the unit cell are zero.
+ */
+ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * Please note that as allowed in the
+ * documentation of the base class,
+ * this function does not implement
+ * the setting up of the local mass
+ * matrix in three space dimensions
+ * because of too high computational
+ * costs. The specified exception
+ * is thrown instead.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ FullMatrix<double> &local_mass_matrix) const;
+
+ private:
+
+ /**
+ * This function is simply singled out of
+ * the constructor; it sets up the
+ * #restriction# and #prolongation#
+ * matrices. Since we have two constructors
+ * which need this functionality, we
+ * provide a single function for this.
+ */
+ void initialize_matrices ();
+};
+
+
+/**
+ * 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}
+ *
+ * \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
+ *
+ * *-------* *-------*
+ * /| | / /|
+ * / | 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
+ */
+template <int dim>
+class FEDG_P2 : public FEQ1Mapping<dim>
+{
+ public:
+ /**
+ * Constructor
+ */
+ FEDG_P2 ();
+
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the tensor of second derivatives
+ * of the #i#th shape function at
+ * point #p# on the unit cell.
+ */
+ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * Please note that as allowed in the
+ * documentation of the base class,
+ * this function does not implement
+ * the setting up of the local mass
+ * matrix in three space dimensions
+ * because of too high computational
+ * costs. The specified exception
+ * is thrown instead.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ FullMatrix<double> &local_mass_matrix) const;
+
+ private:
+
+ /**
+ * This function is simply singled out of
+ * the constructor; it sets up the
+ * #restriction# and #prolongation#
+ * matrices. Since we have two constructors
+ * which need this functionality, we
+ * provide a single function for this.
+ */
+ void initialize_matrices ();
+};
+
+
+/**
+ * Subparametric Q3 finite element in #dim# space dimensions.
+ * A Q1 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--1
+ * \end{verbatim}
+ *
+ * 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
+ */
+template <int dim>
+class FEDG_P3 : public FEQ1Mapping<dim>
+{
+ public:
+ /**
+ * Constructor
+ */
+ FEDG_P3 ();
+
+ public:
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the tensor of second derivatives
+ * of the #i#th shape function at
+ * point #p# on the unit cell.
+ */
+ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ *
+ * Please note that as allowed in the
+ * documentation of the base class,
+ * this function does not implement
+ * the setting up of the local mass
+ * matrix in three space dimensions
+ * because of too high computational
+ * costs. The specified exception
+ * is thrown instead.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ FullMatrix<double> &local_mass_matrix) const;
+
+ private:
+
+ /**
+ * This function is simply singled out of
+ * the constructor; it sets up the
+ * #restriction# and #prolongation#
+ * matrices. Since we have two constructors
+ * which need this functionality, we
+ * provide a single function for this.
+ */
+ void initialize_matrices ();
+};
+
+
+/**
+ * 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}
+ *
+ * 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.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class FEDG_P4 : public FEQ1Mapping<dim>
+{
+ public:
+ /**
+ * Constructor
+ */
+ FEDG_P4 ();
+
+ public:
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual double shape_value(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ */
+ virtual Tensor<1,dim> shape_grad(const unsigned int i,
+ const Point<dim>& p) const;
+
+ /**
+ * Return the tensor of second derivatives
+ * of the #i#th shape function at
+ * point #p# on the unit cell.
+ */
+ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_unit_support_points (vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on this function.
+ */
+ virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
+ vector<Point<dim> > &support_points) const;
+
+ /**
+ * Refer to the base class for detailed
+ * information on what this function does.
+ *
+ * Please note that as allowed in the
+ * documentation of the base class,
+ * this function does not implement
+ * the setting up of the local mass
+ * matrix in three space dimensions
+ * because of too high computational
+ * costs. The specified exception
+ * is thrown instead.
+ */
+ virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+ FullMatrix<double> &local_mass_matrix) const;
+
+ private:
+
+ /**
+ * This function is simply singled out of
+ * the constructor; it sets up the
+ * #restriction# and #prolongation#
+ * matrices. Since we have two constructors
+ * which need this functionality, we
+ * provide a single function for this.
+ */
+ void initialize_matrices ();
+};
+
+
+/*---------------------------- fe_lib.h ---------------------------*/
+
+#endif
+/*---------------------------- fe_lib.h ---------------------------*/
+
+
--- /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_P1<deal_II_dimension>::initialize_matrices ();
+
+
+#if deal_II_dimension == 1
+
+template <>
+FEDG_P1<1>::FEDG_P1 () :
+ FEQ1Mapping<1> (0, 2, 0, 0, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P1<1>::initialize_matrices ()
+{
+ // for restriction and prolongation matrices:
+ // note that we do not add up all the
+ // contributions since then we would get
+ // two summands per vertex in 1d (four
+ // in 2d, etc), but only one per line dof.
+ // We could accomplish for that by dividing
+ // the vertex dof values by 2 (4, etc), but
+ // would get into trouble at the boundary
+ // of the domain since there only one
+ // cell contributes to a vertex. Rather,
+ // we do not add up the contributions but
+ // set them right into the matrices!
+ restriction[0](0,0) = 1.0;
+ restriction[1](1,1) = 1.0;
+
+ prolongation[0](0,0) = 1.0;
+ prolongation[0](1,0) = 1./2.;
+ prolongation[0](1,1) = 1./2.;
+
+ prolongation[1](0,0) = 1./2.;
+ prolongation[1](0,1) = 1./2.;
+ prolongation[1](1,1) = 1.0;
+};
+
+
+template <>
+double
+FEDG_P1<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.-p(0);
+ case 1: return p(0);
+ }
+ return 0.;
+}
+
+
+template <>
+inline
+Tensor<1,1>
+FEDG_P1<1>::shape_grad(const unsigned int i,
+ const Point<1>&) 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.);
+ }
+ return Point<1>();
+};
+
+
+template <>
+inline
+Tensor<2,1>
+FEDG_P1<1>::shape_grad_grad (const unsigned int i,
+ const Point<1> &) const
+{
+ Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+ // second derivatives on the unit cell
+ // are always zero
+ return Tensor<2,1>();
+};
+
+
+template <>
+void FEDG_P1<1>::get_unit_support_points (vector<Point<1> > &support_points) const
+{
+ FiniteElement<1>::get_unit_support_points (support_points);
+};
+
+
+template <>
+void FEDG_P1<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_P1<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
+ vector<Point<1> > &) const
+{
+ Assert (false, ExcInternalError());
+};
+
+
+template <>
+void FEDG_P1<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
+ 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));
+
+ const double h = cell->vertex(1)(0) - cell->vertex(0)(0);
+ Assert (h>0, ExcJacobiDeterminantHasWrongSign());
+
+ local_mass_matrix(0,0) = local_mass_matrix(1,1) = 1./3.*h;
+ local_mass_matrix(0,1) = local_mass_matrix(1,0) = 1./6.*h;
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+FEDG_P1<2>::FEDG_P1 () :
+ FEQ1Mapping<2> (0, 0, 3, 0, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P1<2>::initialize_matrices ()
+{
+ Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P1<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);
+ }
+ return 0.;
+};
+
+
+template <>
+inline
+Tensor<1,2>
+FEDG_P1<2>::shape_grad (const unsigned int i,
+ const Point<2>&) 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);
+ }
+ return Point<2> ();
+};
+
+
+template <>
+inline
+Tensor<2,2>
+FEDG_P1<2>::shape_grad_grad (const unsigned int i,
+ const Point<2> &) const
+{
+ Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+ return Tensor<2,2>();
+};
+
+
+template <>
+void FEDG_P1<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &,
+ 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));
+};
+
+
+template <>
+void FEDG_P1<2>::get_unit_support_points (vector<Point<2> > &unit_points) const
+{
+ 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);
+};
+
+
+#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+FEDG_P1<3>::FEDG_P1 () :
+ FEQ1Mapping<3> (0, 0, 0, 4, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P1<3>::initialize_matrices ()
+{
+ Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P1<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_P1<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_P1<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_P1<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_P1<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_P1<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_P1<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_P1<deal_II_dimension>;
--- /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_P2<deal_II_dimension>::initialize_matrices ();
+
+
+#if deal_II_dimension == 1
+
+template <>
+FEDG_P2<1>::FEDG_P2 () :
+ FEQ1Mapping<1> (0, 3, 0, 0, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P2<1>::initialize_matrices ()
+{
+ Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+double
+FEDG_P2<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);
+ }
+ return 0.;
+}
+
+
+template <>
+inline
+Tensor<1,1>
+FEDG_P2<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));
+
+ }
+ return Point<1>();
+};
+
+
+template <>
+inline
+Tensor<2,1>
+FEDG_P2<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_P2<1>::get_unit_support_points (vector<Point<1> > &support_points) const
+{
+ FiniteElement<1>::get_unit_support_points (support_points);
+};
+
+
+template <>
+void FEDG_P2<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_P2<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
+ vector<Point<1> > &) const
+{
+ Assert (false, ExcInternalError());
+};
+
+
+template <>
+void FEDG_P2<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_P2<2>::FEDG_P2 () :
+ FEQ1Mapping<2> (0, 0, 6, 0, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P2<2>::initialize_matrices ()
+{
+ Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P2<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);
+ }
+ return 0.;
+};
+
+
+template <>
+inline
+Tensor<1,2>
+FEDG_P2<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));
+ }
+ return Point<2> ();
+};
+
+
+template <>
+inline
+Tensor<2,2>
+FEDG_P2<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_P2<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &,
+ FullMatrix<double> &) const
+{
+ Assert(false, ExcNotImplemented ());
+};
+
+
+template <>
+void FEDG_P2<2>::get_unit_support_points (vector<Point<2> > &unit_points) const
+{
+ 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_P2<3>::FEDG_P2 () :
+ FEQ1Mapping<3> (0, 0, 0, 4, 1,
+ vector<bool> (1, true))
+{
+ initialize_matrices ();
+};
+
+
+template <>
+void FEDG_P2<3>::initialize_matrices ()
+{
+ Assert(false, ExcNotImplemented());
+};
+
+
+template <>
+inline
+double
+FEDG_P2<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_P2<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_P2<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_P2<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_P2<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_P2<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_P2<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_P2<deal_II_dimension>;