]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
first step towards FEDG_Pk
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jun 2000 20:51:38 +0000 (20:51 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Jun 2000 20:51:38 +0000 (20:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@2999 0785d39b-7218-0410-832d-ea1e28bc413d

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

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 (file)
index 0000000..2e407a2
--- /dev/null
@@ -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 <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     ---------------------------*/
+
+
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 (file)
index 0000000..aa65462
--- /dev/null
@@ -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 <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>;
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 (file)
index 0000000..a3c3cc9
--- /dev/null
@@ -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 <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>;

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.