]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add Rannacher-Turek element. 1349/head
authorPatrick Esser <Patrick.Esser@gmx.net>
Tue, 30 Jun 2015 13:58:34 +0000 (15:58 +0200)
committerPatrick Esser <Patrick.Esser@gmx.net>
Mon, 17 Aug 2015 14:48:52 +0000 (16:48 +0200)
Lowest order Rannacher Turek polynomials and FE in 2D and tests.

15 files changed:
doc/news/changes.h
include/deal.II/base/polynomials_rannacher_turek.h [new file with mode: 0644]
include/deal.II/fe/fe_rannacher_turek.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/polynomials_rannacher_turek.cc [new file with mode: 0644]
source/base/polynomials_rannacher_turek.inst.in [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/fe_poly.cc
source/fe/fe_poly.inst.in
source/fe/fe_rannacher_turek.cc [new file with mode: 0644]
source/fe/fe_rannacher_turek.inst.in [new file with mode: 0644]
tests/fe/fe_rannacher_turek_01.cc [new file with mode: 0644]
tests/fe/fe_rannacher_turek_01.output [new file with mode: 0644]
tests/fe/fe_rannacher_turek_02.cc [new file with mode: 0644]
tests/fe/fe_rannacher_turek_02.output [new file with mode: 0644]

index b4f44d535b74b60f55c8c9c70ad93eadb61004da..ebcc17138e95e1dad1eff04ff2e9d699cc0e3930 100644 (file)
@@ -93,6 +93,12 @@ inconvenience this causes.
   (Jason Sheldon, Wolfgang Bangerth, 2015/08/13)
   </li>
 
+  <li> New: FE_RannacherTurek describes a discontinuous FiniteElement
+  with vanishing mean values of jumps across faces.
+  <br>
+  (Patrick Esser, 2015/08/17)
+  </li>
+
   <li> New: FE_Q_Bubbles describes a FiniteElement based on FE_Q 
   enriched by bubble functions.
   <br>
diff --git a/include/deal.II/base/polynomials_rannacher_turek.h b/include/deal.II/base/polynomials_rannacher_turek.h
new file mode 100644 (file)
index 0000000..64302be
--- /dev/null
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii__polynomials_rannacher_turek_h
+#define dealii__polynomials_rannacher_turek_h
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/**
+ * Basis for polynomial space on the unit square used for lowest order
+ * Rannacher Turek element.
+ *
+ * The i-th basis function is the dual basis element corresponding to
+ * the dof which evaluates the function's mean value across the i-th
+ * face. The numbering can be found in GeometryInfo.
+ *
+ * @ingroup Polynomials
+ * @author Patrick Esser
+ * @date 2015
+ **/
+template <int dim>
+class PolynomialsRannacherTurek
+{
+public:
+  /**
+   * Dimension we are working in.
+   */
+  static const unsigned int dimension = dim;
+
+  /**
+   * Constructor, checking that the basis is implemented in this
+   * dimension.
+   */
+  PolynomialsRannacherTurek();
+
+  /** Value of basis function @p i at @p p.
+  */
+  double compute_value(const unsigned int i,
+                       const Point<dim> &p) const;
+  /** Gradient of basis function @p i at @p p.
+  */
+  Tensor<1, dim> compute_grad(const unsigned int i,
+                              const Point<dim> &p) const;
+  /** Gradient of gradient of basis function @p i at @p p.
+  */
+  Tensor<2, dim> compute_grad_grad(const unsigned int i,
+                                   const Point<dim> &p) const;
+
+  /**
+   * Compute values and derivatives of all basis functions at @p
+   * unit_point.
+   *
+   * Size of the vectors must be either equal to the number of
+   * polynomials or zero. A size of zero means that we are not
+   * computing the vector entries.
+   */
+  void compute(const Point<dim> &unit_point,
+               std::vector<double> &values,
+               std::vector<Tensor<1, dim> > &grads,
+               std::vector<Tensor<2, dim> > &grad_grads) const;
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h
new file mode 100644 (file)
index 0000000..52dc09e
--- /dev/null
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii__fe_rannacher_turek_h
+#define dealii__fe_rannacher_turek_h
+
+#include <deal.II/base/polynomials_rannacher_turek.h>
+#include <deal.II/fe/fe_poly.h>
+#include <deal.II/fe/fe_base.h>
+#include <string>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/**
+ * Implementation of Rannacher-Turek elements. Functions generated by this
+ * element will be discontinuous, but their jump along faces is mean value
+ * free.
+ *
+ * Implemented only in dimension 2, lowest order, without hanging nodes and
+ * restriction/prolongation.
+ *
+ * <h3>Interpolation</h3>
+ *
+ * <h4>Node values</h4>
+ * The @ref GlossNodes "node values" are moments on faces.
+ *
+ * <h4>Generalized support points</h4>
+ * To calculate the node values, we are using a QGauss rule on each face.
+ * By default, we are using a two point rule to integrate Rannacher-Turek
+ * functions exactly. But in order to be able to interpolate other
+ * functions with sufficient accuracy, the number of quadrature points
+ * used on a face can be adjusted in the constructor.
+ *
+ * @ingroup fe
+ * @author Patrick Esser
+ * @date 2015
+ */
+template <int dim>
+class FE_RannacherTurek : public FE_Poly<PolynomialsRannacherTurek<dim>, dim>
+{
+public:
+  /**
+   * Constructor for Rannacher-Turek element of degree @p degree, using
+   * @p n_face_support_points quadrature points on each face for
+   * interpolation. Notice that the element of degree 0 contains
+   * polynomials of degree 2.
+   *
+   * Only implemented for degree 0 in 2D.
+   */
+  FE_RannacherTurek(const unsigned int degree = 0,
+                    const unsigned int n_face_support_points = 2);
+
+  virtual std::string get_name() const;
+  virtual FiniteElement<dim> *clone() const;
+
+  virtual void interpolate(
+    std::vector<double> &local_dofs,
+    const std::vector<double> &values) const;
+private:
+  /**
+   * Degree of this element.
+   */
+  const unsigned int degree;
+  /**
+   * The number of quadrature points used on each face to evaluate node
+   * functionals during interpolation.
+   */
+  const unsigned int n_face_support_points;
+  /**
+   * The weights used on the faces to evaluate node functionals.
+   */
+  std::vector<double> weights;
+
+  /**
+   * Compute generalized support points and their weights.
+   */
+  void initialize_support_points();
+  /**
+   * Return information about degrees of freedom per object as needed
+   * during construction.
+   */
+  std::vector<unsigned int> get_dpo_vector();
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 9102a1e0ecc83ff8a6425f678e3ce20d498c51f3..1f45f6ed84b7324292994bd88733d3b327ee8dba 100644 (file)
@@ -52,6 +52,7 @@ SET(_src
   polynomial_space.cc
   polynomials_p.cc
   polynomials_piecewise.cc
+  polynomials_rannacher_turek.cc
   polynomials_raviart_thomas.cc
   quadrature.cc
   quadrature_lib.cc
@@ -71,6 +72,7 @@ SET(_src
 
 SET(_inst
   data_out_base.inst.in
+  polynomials_rannacher_turek.inst.in
   time_stepping.inst.in
   )
 
diff --git a/source/base/polynomials_rannacher_turek.cc b/source/base/polynomials_rannacher_turek.cc
new file mode 100644 (file)
index 0000000..9167154
--- /dev/null
@@ -0,0 +1,172 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/polynomials_rannacher_turek.h>
+#include <deal.II/base/geometry_info.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim>
+PolynomialsRannacherTurek<dim>::PolynomialsRannacherTurek()
+{
+  Assert(dim == 2, ExcNotImplemented());
+}
+
+
+
+template <int dim>
+double PolynomialsRannacherTurek<dim>::compute_value(const unsigned int i,
+                                                     const Point<dim> &p) const
+{
+  Assert(dim == 2, ExcNotImplemented());
+  if (i == 0)
+    {
+      return (0.75 - 2.5*p(0) + 1.5*p(1) + 1.5*(p(0)*p(0) - p(1)*p(1)));
+    }
+  else if (i == 1)
+    {
+      return (-0.25 - 0.5*p(0) + 1.5*p(1) + 1.5*(p(0)*p(0) - p(1)*p(1)));
+    }
+  else if (i == 2)
+    {
+      return (0.75 + 1.5*p(0) - 2.5*p(1) - 1.5*(p(0)*p(0) - p(1)*p(1)));
+    }
+  else if (i == 3)
+    {
+      return (-0.25 + 1.5*p(0) - 0.5*p(1) - 1.5*(p(0)*p(0) - p(1)*p(1)));
+    }
+
+  Assert(false, ExcNotImplemented());
+  return 0;
+}
+
+
+
+template <int dim>
+Tensor<1, dim> PolynomialsRannacherTurek<dim>::compute_grad(
+  const unsigned int i,
+  const Point<dim> &p) const
+{
+  Assert(dim == 2, ExcNotImplemented());
+  Tensor<1, dim> grad;
+  if (i == 0)
+    {
+      grad[0] = -2.5 + 3*p(0);
+      grad[1] = 1.5 - 3*p(1);
+    }
+  else if (i == 1)
+    {
+      grad[0] = -0.5 + 3.0*p(0);
+      grad[1] = 1.5 - 3.0*p(1);
+    }
+  else if (i == 2)
+    {
+      grad[0] = 1.5 - 3.0*p(0);
+      grad[1] = -2.5 + 3.0*p(1);
+    }
+  else if (i == 3)
+    {
+      grad[0] = 1.5 - 3.0*p(0);
+      grad[1] = -0.5 + 3.0*p(1);
+    }
+  else
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+  return grad;
+}
+
+
+
+template <int dim> Tensor<2, dim>
+PolynomialsRannacherTurek<dim>::compute_grad_grad(const unsigned int i,
+                                                  const Point<dim> & /*p*/)
+const
+{
+  Assert(dim == 2, ExcNotImplemented());
+  Tensor<2, dim> grad_grad;
+  if (i == 0)
+    {
+      grad_grad[0][0] = 3;
+      grad_grad[0][1] = 0;
+      grad_grad[1][0] = 0;
+      grad_grad[1][1] = -3;
+    }
+  else if (i == 1)
+    {
+      grad_grad[0][0] = 3;
+      grad_grad[0][1] = 0;
+      grad_grad[1][0] = 0;
+      grad_grad[1][1] = -3;
+    }
+  else if (i == 2)
+    {
+      grad_grad[0][0] = -3;
+      grad_grad[0][1] = 0;
+      grad_grad[1][0] = 0;
+      grad_grad[1][1] = 3;
+    }
+  else if (i == 3)
+    {
+      grad_grad[0][0] = -3;
+      grad_grad[0][1] = 0;
+      grad_grad[1][0] = 0;
+      grad_grad[1][1] = 3;
+    }
+  return grad_grad;
+}
+
+
+
+template <int dim>
+void PolynomialsRannacherTurek<dim>::compute(
+  const Point<dim> &unit_point,
+  std::vector<double> &values,
+  std::vector<Tensor<1, dim> > &grads,
+  std::vector<Tensor<2, dim> > &grad_grads) const
+{
+  const unsigned int n_pols = dealii::GeometryInfo<dim>::faces_per_cell;
+  Assert(values.size() == n_pols || values.size() == 0,
+         ExcDimensionMismatch(values.size(), n_pols));
+  Assert(grads.size() == n_pols || grads.size() == 0,
+         ExcDimensionMismatch(grads.size(), n_pols));
+  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
+         ExcDimensionMismatch(grad_grads.size(), n_pols));
+
+  for (unsigned int i = 0; i < n_pols; ++i)
+    {
+      if (values.size() != 0)
+        {
+          values[i] = compute_value(i, unit_point);
+        }
+      if (grads.size() != 0)
+        {
+          grads[i] = compute_grad(i, unit_point);
+        }
+      if (grad_grads.size() != 0)
+        {
+          grad_grads[i] = compute_grad_grad(i, unit_point);
+        }
+    }
+}
+
+
+// explicit instantiations
+#include "polynomials_rannacher_turek.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/base/polynomials_rannacher_turek.inst.in b/source/base/polynomials_rannacher_turek.inst.in
new file mode 100644 (file)
index 0000000..32c95b6
--- /dev/null
@@ -0,0 +1,22 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template class PolynomialsRannacherTurek<deal_II_dimension>;
+  }
+
index 73fee37d1e3a7b0065367cafbb184bdeae5e23c0..228cd4e686c1170d200824500a49dc77dbdfe916 100644 (file)
@@ -39,6 +39,7 @@ SET(_src
   fe_q_dg0.cc
   fe_q_hierarchical.cc
   fe_q_iso_q1.cc
+  fe_rannacher_turek.cc
   fe_raviart_thomas.cc
   fe_raviart_thomas_nodal.cc
   fe_system.cc
@@ -78,6 +79,7 @@ SET(_inst
   fe_q_hierarchical.inst.in
   fe_q.inst.in
   fe_q_iso_q1.inst.in
+  fe_rannacher_turek.inst.in
   fe_raviart_thomas.inst.in
   fe_raviart_thomas_nodal.inst.in
   fe_system.inst.in
index ebd25e79063c0d4411c2710ae769b69402d73632..9fe3a4650d714b1a181050a9822440d7e5631de3 100644 (file)
@@ -20,6 +20,7 @@
 #include <deal.II/base/tensor_product_polynomials_bubbles.h>
 #include <deal.II/base/polynomials_p.h>
 #include <deal.II/base/polynomials_piecewise.h>
+#include <deal.II/base/polynomials_rannacher_turek.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_poly.templates.h>
index dcd33ed098ef6011c91da0a64ee403a57914e940..d9dd85373abf9efd91f22c35fe888f3fec17b2a3 100644 (file)
@@ -24,6 +24,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     template class FE_Poly<TensorProductPolynomials<deal_II_dimension,Polynomials::PiecewisePolynomial<double> >, deal_II_dimension, deal_II_space_dimension>;
     template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
     template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
+    template class FE_Poly<PolynomialsRannacherTurek<deal_II_dimension>, deal_II_dimension, deal_II_space_dimension>;
 #endif
   }
 
diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc
new file mode 100644 (file)
index 0000000..4f901e5
--- /dev/null
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/fe/fe_rannacher_turek.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <algorithm>
+
+#include <sstream>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim>
+FE_RannacherTurek<dim>::FE_RannacherTurek(const unsigned int degree,
+                                          const unsigned int n_face_support_points) :
+  FE_Poly<PolynomialsRannacherTurek<dim>, dim>(
+    PolynomialsRannacherTurek<dim>(),
+    FiniteElementData<dim>(this->get_dpo_vector(),
+                           1,
+                           2,
+                           FiniteElementData<dim>::L2),
+    std::vector<bool>(4, false), // restriction not implemented
+    std::vector<ComponentMask>(4, std::vector<bool>(1, true))),
+  degree(degree),
+  n_face_support_points(n_face_support_points)
+{
+  Assert(dim == 2, ExcNotImplemented());
+  Assert(degree == 0, ExcNotImplemented());
+  this->initialize_support_points();
+}
+
+
+
+template <int dim>
+std::vector<unsigned int> FE_RannacherTurek<dim>::get_dpo_vector()
+{
+  std::vector<unsigned int> dpo(dim + 1, 0);
+  dpo[dim - 1] = 1;
+
+  return dpo;
+}
+
+
+
+template <int dim>
+std::string FE_RannacherTurek<dim>::get_name() const
+{
+  std::ostringstream namebuf;
+  namebuf << "FE_RannacherTurek"
+          << "<" << dim << ">"
+          << "(" << this->degree << ", " << this-n_face_support_points << ")";
+  return namebuf.str();
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *FE_RannacherTurek<dim>::clone() const
+{
+  return new FE_RannacherTurek<dim>(this->degree, this->n_face_support_points);
+}
+
+
+
+template <int dim>
+void FE_RannacherTurek<dim>::initialize_support_points()
+{
+  Assert(dim == 2, ExcNotImplemented());
+  dealii::QGauss<dim-1> face_quadrature(this->n_face_support_points);
+  this->weights = face_quadrature.get_weights();
+  this->generalized_support_points.resize(4*face_quadrature.size());
+  for (unsigned int q = 0;
+       q < face_quadrature.size();
+       ++q)
+    {
+      this->generalized_support_points[0*face_quadrature.size() + q] =
+        dealii::Point<dim>(0, 1 - face_quadrature.point(q)(0));
+      this->generalized_support_points[1*face_quadrature.size() + q] =
+        dealii::Point<dim>(1, 1 - face_quadrature.point(q)(0));
+      this->generalized_support_points[2*face_quadrature.size() + q] =
+        dealii::Point<dim>(face_quadrature.point(q)(0), 0);
+      this->generalized_support_points[3*face_quadrature.size() + q] =
+        dealii::Point<dim>(face_quadrature.point(q)(0), 1);
+    }
+}
+
+
+
+template <int dim>
+void FE_RannacherTurek<dim>::interpolate(
+  std::vector<double> &local_dofs,
+  const std::vector<double> &values) const
+{
+  AssertDimension(values.size(), this->generalized_support_points.size());
+  AssertDimension(local_dofs.size(), this->dofs_per_cell);
+
+  const unsigned int q_points_per_face = this->weights.size();
+  std::fill(local_dofs.begin(), local_dofs.end(), 0.0);
+
+  std::vector<double>::const_iterator value = values.begin();
+  for (unsigned int face = 0;
+       face < dealii::GeometryInfo<dim>::faces_per_cell;
+       ++face)
+    {
+      for (unsigned int q = 0;
+           q < q_points_per_face;
+           ++q)
+        {
+          local_dofs[face] += (*value) * this->weights[q];
+          ++value;
+        }
+    }
+}
+
+
+
+// explicit instantiations
+#include "fe_rannacher_turek.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_rannacher_turek.inst.in b/source/fe/fe_rannacher_turek.inst.in
new file mode 100644 (file)
index 0000000..5c165e2
--- /dev/null
@@ -0,0 +1,22 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template class FE_RannacherTurek<deal_II_dimension>;
+  }
+
diff --git a/tests/fe/fe_rannacher_turek_01.cc b/tests/fe/fe_rannacher_turek_01.cc
new file mode 100644 (file)
index 0000000..55b0d66
--- /dev/null
@@ -0,0 +1,151 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Interfaces being tested
+#include <deal.II/fe/fe_rannacher_turek.h>
+#include <deal.II/base/polynomials_rannacher_turek.h>
+// Interfaces needed for testing
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include "../tests.h"
+
+// Known results for the Rannacher Turek element. All output should be zero.
+
+void test_known_values()
+{
+  PolynomialsRannacherTurek<2> pols;
+  Point<2> p(0.5, 0.5);
+  for (unsigned int i = 0; i < 4; ++i)
+    {
+      deallog << pols.compute_value(i, p) - 0.25 << std::endl;
+    }
+}
+
+void test_n_dofs()
+{
+  FE_RannacherTurek<2> fe_ratu;
+
+  Triangulation<2> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  tria.refine_global(1);
+
+  DoFHandler<2> dofh;
+  dofh.initialize(tria, fe_ratu);
+
+  deallog << dofh.n_dofs() - 12 << std::endl;
+}
+
+void test_nodal_matrix()
+{
+  FE_RannacherTurek<2> fe;
+
+  FullMatrix<double> N(4, 4);
+  // compute_node_matrix does not check for number of components but
+  // simply assumes that there are dim components. Thus we do it ourselves.
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  const std::vector<Point<2> > &points = fe.get_generalized_support_points();
+  std::vector<double> values(points.size());
+  std::vector<double> local_dofs(n_dofs);
+
+  for (unsigned int i = 0; i < n_dofs; ++i)
+    {
+      for (unsigned int k = 0; k < values.size(); ++k)
+        values[k] = fe.shape_value(i, points[k]);
+      fe.interpolate(local_dofs, values);
+
+      for (unsigned int j=0; j < n_dofs; ++j)
+        N(j,i) = local_dofs[j];
+    }
+
+  for (unsigned int i = 0; i < 4; ++i)
+    {
+      for (unsigned int j = 0; j < 4; ++j)
+        {
+          if (i == j)
+            {
+              deallog << N(i, j) - 1.0;
+            }
+          else
+            {
+              deallog << N(i, j) - 0.0;
+            }
+          if (j + 1 < 4)
+            deallog << " ";
+          else
+            deallog << std::endl;
+        }
+    }
+}
+
+void test_interpolation()
+{
+  Triangulation<2> tr;
+  GridGenerator::hyper_cube(tr, -1, 1);
+  tr.refine_global(2);
+
+  FE_RannacherTurek<2> fe;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+
+  DoFHandler<2> dofh;
+  dofh.initialize(tr, fe);
+
+  Vector<double> input_vector(dofh.n_dofs());
+  for (unsigned int i = 0; i < input_vector.size(); ++i)
+    {
+      input_vector[i] = double(i);
+    }
+
+  Quadrature<2> quadrature(fe.get_generalized_support_points());
+  FEValues<2> fev(fe, quadrature, update_values | update_JxW_values);
+
+  typedef DoFHandler<2>::cell_iterator cell_it;
+  cell_it cell = dofh.begin_active();
+  for (; cell != dofh.end(); ++cell)
+    {
+      fev.reinit(cell);
+
+      std::vector<double> values(quadrature.size());
+      fev.get_function_values(input_vector, values);
+
+      std::vector<double> interpolated_local_dofs(n_dofs);
+      fe.interpolate(interpolated_local_dofs, values);
+
+      Vector<double> local_dofs(n_dofs);
+      cell->get_dof_values(input_vector, local_dofs);
+
+      for (unsigned int j = 0; j < n_dofs; ++j)
+        {
+          deallog << local_dofs[j] - interpolated_local_dofs[j] << " ";
+        }
+      deallog << std::endl;
+    }
+}
+
+int main()
+{
+  initlog();
+
+  test_known_values();
+  test_n_dofs();
+  test_nodal_matrix();
+  test_interpolation();
+
+  return 0;
+}
diff --git a/tests/fe/fe_rannacher_turek_01.output b/tests/fe/fe_rannacher_turek_01.output
new file mode 100644 (file)
index 0000000..46c8ddb
--- /dev/null
@@ -0,0 +1,26 @@
+
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0
+DEAL::0.00000 -1.04083e-16 0.00000 5.55112e-17
+DEAL::-1.11022e-16 -2.22045e-16 1.11022e-16 1.66533e-16
+DEAL::0.00000 5.55112e-17 0.00000 -1.04083e-16
+DEAL::1.11022e-16 1.66533e-16 -1.11022e-16 -2.22045e-16
+DEAL::-5.55112e-17 -6.66134e-16 0.00000 8.88178e-16 
+DEAL::1.11022e-16 -8.88178e-16 0.00000 1.77636e-15 
+DEAL::0.00000 8.88178e-16 4.44089e-16 0.00000 
+DEAL::0.00000 0.00000 8.88178e-16 0.00000 
+DEAL::4.44089e-16 0.00000 0.00000 1.77636e-15 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::3.55271e-15 3.55271e-15 1.77636e-15 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::3.55271e-15 3.55271e-15 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 
diff --git a/tests/fe/fe_rannacher_turek_02.cc b/tests/fe/fe_rannacher_turek_02.cc
new file mode 100644 (file)
index 0000000..22a26f2
--- /dev/null
@@ -0,0 +1,66 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Interfaces being tested
+#include <deal.II/fe/fe_rannacher_turek.h>
+#include <deal.II/base/polynomials_rannacher_turek.h>
+// Interfaces needed for testing
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include "../tests.h"
+
+
+// Regression test for some values of Rannacher Turek element.
+
+void test_values()
+{
+  Triangulation<2> tria;
+  GridGenerator::hyper_cube(tria, 0.0, 1.0);
+
+  FE_RannacherTurek<2> fe;
+  DoFHandler<2> dofh;
+  dofh.initialize(tria, fe);
+
+  QGauss<2> quadrature(8);
+  FEValues<2> fev(fe, quadrature, update_values | update_gradients | update_JxW_values);
+  fev.reinit(dofh.begin_active());
+
+  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+    {
+      for (unsigned int q = 0; q < quadrature.size(); ++q)
+        {
+          deallog << fev.shape_value(i, q) << " ";
+          for (unsigned int d = 0; d < 2; ++d)
+            {
+              deallog << fev.shape_grad(i, q)[d] << " ";
+            }
+        }
+      deallog << std::endl;
+    }
+}
+
+int main()
+{
+  initlog();
+
+  test_values();
+
+  return 0;
+}
diff --git a/tests/fe/fe_rannacher_turek_02.output b/tests/fe/fe_rannacher_turek_02.output
new file mode 100644 (file)
index 0000000..5eebd5b
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::0.730145 -2.44043 1.44043 0.540529 -2.19500 1.44043 0.270527 -1.78830 1.44043 0.00852669 -1.27515 1.44043 -0.174908 -0.724848 1.44043 -0.255006 -0.211701 1.44043 -0.256138 0.195000 1.44043 -0.230145 0.440435 1.44043 0.837950 -2.44043 1.19500 0.648333 -2.19500 1.19500 0.378331 -1.78830 1.19500 0.116331 -1.27515 1.19500 -0.0671033 -0.724848 1.19500 -0.147201 -0.211701 1.19500 -0.148333 0.195000 1.19500 -0.122340 0.440435 1.19500 0.972385 -2.44043 0.788299 0.782768 -2.19500 0.788299 0.512766 -1.78830 0.788299 0.250766 -1.27515 0.788299 0.0673317 -0.724848 0.788299 -0.0127662 -0.211701 0.788299 -0.0138983 0.195000 0.788299 0.0120947 0.440435 0.788299 1.06334 -2.44043 0.275152 0.873719 -2.19500 0.275152 0.603717 -1.78830 0.275152 0.341717 -1.27515 0.275152 0.158283 -0.724848 0.275152 0.0781848 -0.211701 0.275152 0.0770527 0.195000 0.275152 0.103046 0.440435 0.275152 1.06334 -2.44043 -0.275152 0.873719 -2.19500 -0.275152 0.603717 -1.78830 -0.275152 0.341717 -1.27515 -0.275152 0.158283 -0.724848 -0.275152 0.0781848 -0.211701 -0.275152 0.0770527 0.195000 -0.275152 0.103046 0.440435 -0.275152 0.972385 -2.44043 -0.788299 0.782768 -2.19500 -0.788299 0.512766 -1.78830 -0.788299 0.250766 -1.27515 -0.788299 0.0673317 -0.724848 -0.788299 -0.0127662 -0.211701 -0.788299 -0.0138983 0.195000 -0.788299 0.0120947 0.440435 -0.788299 0.837950 -2.44043 -1.19500 0.648333 -2.19500 -1.19500 0.378331 -1.78830 -1.19500 0.116331 -1.27515 -1.19500 -0.0671033 -0.724848 -1.19500 -0.147201 -0.211701 -1.19500 -0.148333 0.195000 -1.19500 -0.122340 0.440435 -1.19500 0.730145 -2.44043 -1.44043 0.540529 -2.19500 -1.44043 0.270527 -1.78830 -1.44043 0.00852669 -1.27515 -1.44043 -0.174908 -0.724848 -1.44043 -0.255006 -0.211701 -1.44043 -0.256138 0.195000 -1.44043 -0.230145 0.440435 -1.44043 
+DEAL::-0.230145 -0.440435 1.44043 -0.256138 -0.195000 1.44043 -0.255006 0.211701 1.44043 -0.174908 0.724848 1.44043 0.00852669 1.27515 1.44043 0.270527 1.78830 1.44043 0.540529 2.19500 1.44043 0.730145 2.44043 1.44043 -0.122340 -0.440435 1.19500 -0.148333 -0.195000 1.19500 -0.147201 0.211701 1.19500 -0.0671033 0.724848 1.19500 0.116331 1.27515 1.19500 0.378331 1.78830 1.19500 0.648333 2.19500 1.19500 0.837950 2.44043 1.19500 0.0120947 -0.440435 0.788299 -0.0138983 -0.195000 0.788299 -0.0127662 0.211701 0.788299 0.0673317 0.724848 0.788299 0.250766 1.27515 0.788299 0.512766 1.78830 0.788299 0.782768 2.19500 0.788299 0.972385 2.44043 0.788299 0.103046 -0.440435 0.275152 0.0770527 -0.195000 0.275152 0.0781848 0.211701 0.275152 0.158283 0.724848 0.275152 0.341717 1.27515 0.275152 0.603717 1.78830 0.275152 0.873719 2.19500 0.275152 1.06334 2.44043 0.275152 0.103046 -0.440435 -0.275152 0.0770527 -0.195000 -0.275152 0.0781848 0.211701 -0.275152 0.158283 0.724848 -0.275152 0.341717 1.27515 -0.275152 0.603717 1.78830 -0.275152 0.873719 2.19500 -0.275152 1.06334 2.44043 -0.275152 0.0120947 -0.440435 -0.788299 -0.0138983 -0.195000 -0.788299 -0.0127662 0.211701 -0.788299 0.0673317 0.724848 -0.788299 0.250766 1.27515 -0.788299 0.512766 1.78830 -0.788299 0.782768 2.19500 -0.788299 0.972385 2.44043 -0.788299 -0.122340 -0.440435 -1.19500 -0.148333 -0.195000 -1.19500 -0.147201 0.211701 -1.19500 -0.0671033 0.724848 -1.19500 0.116331 1.27515 -1.19500 0.378331 1.78830 -1.19500 0.648333 2.19500 -1.19500 0.837950 2.44043 -1.19500 -0.230145 -0.440435 -1.44043 -0.256138 -0.195000 -1.44043 -0.255006 0.211701 -1.44043 -0.174908 0.724848 -1.44043 0.00852669 1.27515 -1.44043 0.270527 1.78830 -1.44043 0.540529 2.19500 -1.44043 0.730145 2.44043 -1.44043 
+DEAL::0.730145 1.44043 -2.44043 0.837950 1.19500 -2.44043 0.972385 0.788299 -2.44043 1.06334 0.275152 -2.44043 1.06334 -0.275152 -2.44043 0.972385 -0.788299 -2.44043 0.837950 -1.19500 -2.44043 0.730145 -1.44043 -2.44043 0.540529 1.44043 -2.19500 0.648333 1.19500 -2.19500 0.782768 0.788299 -2.19500 0.873719 0.275152 -2.19500 0.873719 -0.275152 -2.19500 0.782768 -0.788299 -2.19500 0.648333 -1.19500 -2.19500 0.540529 -1.44043 -2.19500 0.270527 1.44043 -1.78830 0.378331 1.19500 -1.78830 0.512766 0.788299 -1.78830 0.603717 0.275152 -1.78830 0.603717 -0.275152 -1.78830 0.512766 -0.788299 -1.78830 0.378331 -1.19500 -1.78830 0.270527 -1.44043 -1.78830 0.00852669 1.44043 -1.27515 0.116331 1.19500 -1.27515 0.250766 0.788299 -1.27515 0.341717 0.275152 -1.27515 0.341717 -0.275152 -1.27515 0.250766 -0.788299 -1.27515 0.116331 -1.19500 -1.27515 0.00852669 -1.44043 -1.27515 -0.174908 1.44043 -0.724848 -0.0671033 1.19500 -0.724848 0.0673317 0.788299 -0.724848 0.158283 0.275152 -0.724848 0.158283 -0.275152 -0.724848 0.0673317 -0.788299 -0.724848 -0.0671033 -1.19500 -0.724848 -0.174908 -1.44043 -0.724848 -0.255006 1.44043 -0.211701 -0.147201 1.19500 -0.211701 -0.0127662 0.788299 -0.211701 0.0781848 0.275152 -0.211701 0.0781848 -0.275152 -0.211701 -0.0127662 -0.788299 -0.211701 -0.147201 -1.19500 -0.211701 -0.255006 -1.44043 -0.211701 -0.256138 1.44043 0.195000 -0.148333 1.19500 0.195000 -0.0138983 0.788299 0.195000 0.0770527 0.275152 0.195000 0.0770527 -0.275152 0.195000 -0.0138983 -0.788299 0.195000 -0.148333 -1.19500 0.195000 -0.256138 -1.44043 0.195000 -0.230145 1.44043 0.440435 -0.122340 1.19500 0.440435 0.0120947 0.788299 0.440435 0.103046 0.275152 0.440435 0.103046 -0.275152 0.440435 0.0120947 -0.788299 0.440435 -0.122340 -1.19500 0.440435 -0.230145 -1.44043 0.440435 
+DEAL::-0.230145 1.44043 -0.440435 -0.122340 1.19500 -0.440435 0.0120947 0.788299 -0.440435 0.103046 0.275152 -0.440435 0.103046 -0.275152 -0.440435 0.0120947 -0.788299 -0.440435 -0.122340 -1.19500 -0.440435 -0.230145 -1.44043 -0.440435 -0.256138 1.44043 -0.195000 -0.148333 1.19500 -0.195000 -0.0138983 0.788299 -0.195000 0.0770527 0.275152 -0.195000 0.0770527 -0.275152 -0.195000 -0.0138983 -0.788299 -0.195000 -0.148333 -1.19500 -0.195000 -0.256138 -1.44043 -0.195000 -0.255006 1.44043 0.211701 -0.147201 1.19500 0.211701 -0.0127662 0.788299 0.211701 0.0781848 0.275152 0.211701 0.0781848 -0.275152 0.211701 -0.0127662 -0.788299 0.211701 -0.147201 -1.19500 0.211701 -0.255006 -1.44043 0.211701 -0.174908 1.44043 0.724848 -0.0671033 1.19500 0.724848 0.0673317 0.788299 0.724848 0.158283 0.275152 0.724848 0.158283 -0.275152 0.724848 0.0673317 -0.788299 0.724848 -0.0671033 -1.19500 0.724848 -0.174908 -1.44043 0.724848 0.00852669 1.44043 1.27515 0.116331 1.19500 1.27515 0.250766 0.788299 1.27515 0.341717 0.275152 1.27515 0.341717 -0.275152 1.27515 0.250766 -0.788299 1.27515 0.116331 -1.19500 1.27515 0.00852669 -1.44043 1.27515 0.270527 1.44043 1.78830 0.378331 1.19500 1.78830 0.512766 0.788299 1.78830 0.603717 0.275152 1.78830 0.603717 -0.275152 1.78830 0.512766 -0.788299 1.78830 0.378331 -1.19500 1.78830 0.270527 -1.44043 1.78830 0.540529 1.44043 2.19500 0.648333 1.19500 2.19500 0.782768 0.788299 2.19500 0.873719 0.275152 2.19500 0.873719 -0.275152 2.19500 0.782768 -0.788299 2.19500 0.648333 -1.19500 2.19500 0.540529 -1.44043 2.19500 0.730145 1.44043 2.44043 0.837950 1.19500 2.44043 0.972385 0.788299 2.44043 1.06334 0.275152 2.44043 1.06334 -0.275152 2.44043 0.972385 -0.788299 2.44043 0.837950 -1.19500 2.44043 0.730145 -1.44043 2.44043 

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.