From: guido Date: Fri, 9 Jun 2000 20:51:38 +0000 (+0000) Subject: first step towards FEDG_Pk X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1726cd668857d20d9c9295213d3c0db9a4b1d47c;p=dealii-svn.git first step towards FEDG_Pk git-svn-id: https://svn.dealii.org/trunk@2999 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_lib.dgp.h b/deal.II/deal.II/include/fe/fe_lib.dgp.h new file mode 100644 index 0000000000..2e407a24fa --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_lib.dgp.h @@ -0,0 +1,507 @@ +//---------------------------- 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 + +/** + * 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 +class FEDG_P1 : public FEQ1Mapping +{ + 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& 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& 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 &p) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_unit_support_points (vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_support_points (const DoFHandler::cell_iterator &cell, + vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_face_support_points (const DoFHandler::face_iterator &face, + vector > &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::cell_iterator &cell, + FullMatrix &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 +class FEDG_P2 : public FEQ1Mapping +{ + 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& 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& 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 &p) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_unit_support_points (vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_support_points (const DoFHandler::cell_iterator &cell, + vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_face_support_points (const DoFHandler::face_iterator &face, + vector > &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::cell_iterator &cell, + FullMatrix &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 +class FEDG_P3 : public FEQ1Mapping +{ + 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& 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& 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 &p) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_unit_support_points (vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_support_points (const DoFHandler::cell_iterator &cell, + vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_face_support_points (const DoFHandler::face_iterator &face, + vector > &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::cell_iterator &cell, + FullMatrix &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 +class FEDG_P4 : public FEQ1Mapping +{ + 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& 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& 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 &p) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_unit_support_points (vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_support_points (const DoFHandler::cell_iterator &cell, + vector > &support_points) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_face_support_points (const DoFHandler::face_iterator &face, + vector > &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::cell_iterator &cell, + FullMatrix &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 ---------------------------*/ + + diff --git a/deal.II/deal.II/source/fe/fe_lib.dgp1.cc b/deal.II/deal.II/source/fe/fe_lib.dgp1.cc new file mode 100644 index 0000000000..aa65462607 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_lib.dgp1.cc @@ -0,0 +1,381 @@ +//---------------------------- $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 +#include +#include +#include +#include + + +// declare explicit specializations before use: +template <> void FEDG_P1::initialize_matrices (); + + +#if deal_II_dimension == 1 + +template <> +FEDG_P1<1>::FEDG_P1 () : + FEQ1Mapping<1> (0, 2, 0, 0, 1, + vector (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 +inline +Tensor<1,1> +FEDG_P1<1>::shape_grad(const unsigned int i, + const Point<1>&) const +{ + Assert((i, 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(); +}; + + +template <> +void FEDG_P1<1>::get_unit_support_points (vector > &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 > &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 > &) const +{ + Assert (false, ExcInternalError()); +}; + + +template <> +void FEDG_P1<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, + FullMatrix &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 (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 +inline +Tensor<1,2> +FEDG_P1<2>::shape_grad (const unsigned int i, + const Point<2>&) const +{ + Assert((i, 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(); +}; + + +template <> +void FEDG_P1<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &, + FullMatrix &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 > &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 (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 +inline +Tensor<1,3> +FEDG_P1<3>::shape_grad (const unsigned int i, + const Point<3>& p) const +{ + Assert((i, 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 return_value; + return return_value; +}; + + +template <> +void FEDG_P1<3>::get_local_mass_matrix (const DoFHandler<3>::cell_iterator &, + FullMatrix &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 > &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 +void +FEDG_P1::get_support_points (const typename DoFHandler::cell_iterator &cell, + vector > &support_points) const +{ + Assert (support_points.size() == dofs_per_cell, + ExcWrongFieldDimension (support_points.size(), dofs_per_cell)); + + for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) + support_points[vertex] = cell->vertex(vertex); +}; + + +template +void +FEDG_P1::get_face_support_points (const typename DoFHandler::face_iterator &face, + vector > &support_points) const +{ + Assert ((support_points.size() == dofs_per_face) && + (support_points.size() == GeometryInfo::vertices_per_face), + ExcWrongFieldDimension (support_points.size(), + GeometryInfo::vertices_per_face)); + + for (unsigned int vertex=0; vertexvertex(vertex); +}; + + +// explicit instantiations + +template class FEDG_P1; diff --git a/deal.II/deal.II/source/fe/fe_lib.dgp2.cc b/deal.II/deal.II/source/fe/fe_lib.dgp2.cc new file mode 100644 index 0000000000..a3c3cc9d3f --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_lib.dgp2.cc @@ -0,0 +1,363 @@ +//---------------------------- $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 +#include +#include +#include +#include + + +// declare explicit specializations before use: +template <> void FEDG_P2::initialize_matrices (); + + +#if deal_II_dimension == 1 + +template <> +FEDG_P2<1>::FEDG_P2 () : + FEQ1Mapping<1> (0, 3, 0, 0, 1, + vector (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 +inline +Tensor<1,1> +FEDG_P2<1>::shape_grad(const unsigned int i, + const Point<1>&p) const +{ + Assert((i, 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(); +}; + + +template <> +void FEDG_P2<1>::get_unit_support_points (vector > &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 > &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 > &) const +{ + Assert (false, ExcInternalError()); +}; + + +template <> +void FEDG_P2<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, + FullMatrix &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 (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 +inline +Tensor<1,2> +FEDG_P2<2>::shape_grad (const unsigned int i, + const Point<2>& p) const +{ + Assert((i, 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(); +}; + + +template <> +void FEDG_P2<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &, + FullMatrix &) const +{ + Assert(false, ExcNotImplemented ()); +}; + + +template <> +void FEDG_P2<2>::get_unit_support_points (vector > &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 (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 +inline +Tensor<1,3> +FEDG_P2<3>::shape_grad (const unsigned int i, + const Point<3>& p) const +{ + Assert((i, 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 return_value; + return return_value; +}; + + +template <> +void FEDG_P2<3>::get_local_mass_matrix (const DoFHandler<3>::cell_iterator &, + FullMatrix &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 > &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 +void +FEDG_P2::get_support_points (const typename DoFHandler::cell_iterator &cell, + vector > &support_points) const +{ + Assert (support_points.size() == dofs_per_cell, + ExcWrongFieldDimension (support_points.size(), dofs_per_cell)); + + for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) + support_points[vertex] = cell->vertex(vertex); +}; + + +template +void +FEDG_P2::get_face_support_points (const typename DoFHandler::face_iterator &face, + vector > &support_points) const +{ + Assert ((support_points.size() == dofs_per_face) && + (support_points.size() == GeometryInfo::vertices_per_face), + ExcWrongFieldDimension (support_points.size(), + GeometryInfo::vertices_per_face)); + + for (unsigned int vertex=0; vertexvertex(vertex); +}; + + +// explicit instantiations + +template class FEDG_P2;