]> https://gitweb.dealii.org/ - dealii.git/commitdiff
start getting FE_BDM to work in 3D 1245/head
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 2 Aug 2015 20:33:35 +0000 (22:33 +0200)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 5 Aug 2015 20:01:36 +0000 (22:01 +0200)
16 files changed:
include/deal.II/base/polynomial_space.h
include/deal.II/base/polynomials_bdm.h
include/deal.II/fe/fe_bdm.h
include/deal.II/fe/fe_tools.h
source/base/polynomial_space.cc
source/base/polynomials_bdm.cc
source/fe/fe_bdm.cc
source/fe/fe_tools.cc
source/fe/fe_tools.inst.in
tests/fe/bdm_1.cc
tests/fe/bdm_2.cc
tests/fe/bdm_3.cc
tests/fe/bdm_5.cc
tests/fe/bdm_8.cc
tests/fe/bdm_9.cc
tests/fe/fe_data_test.cc

index 486d8a6f53cd517b1059d55114bc3ee2f3f2aab7..3d7baf7afb58019ca0dc4b69321d885651190b8e 100644 (file)
@@ -188,6 +188,9 @@ public:
   /**
    * Static function used in the constructor to compute the number of
    * polynomials.
+   *
+   * @warning The argument `n` is not the maximal degree, but the
+   * number of onedimensional polynomials, thus the degree plus one.
    */
   static unsigned int compute_n_pols (const unsigned int n);
 
index 45a0fce2e0d0f349c9e657e4d3ade6c0bd7c9b5b..393573d05d1302e851d97e6d225fbbf6c5f07f29 100644 (file)
@@ -49,21 +49,34 @@ DEAL_II_NAMESPACE_OPEN
  *
  *   More specifically, for $k=1$, this space has shape functions
  *   @f{align*}
- *     \phi_0 = \begin{array}{cc} 1 \\ 0 \end{array},
- *     \phi_1 = \begin{array}{cc} -\sqrt{3}+2\sqrt{3}x \\ 0 \end{array},
- *     \phi_2 = \begin{array}{cc} -\sqrt{3}+2\sqrt{3}y \\ 0 \end{array},
- *     \phi_3 = \begin{array}{cc} 0 \\ 1 \end{array},
- *     \phi_4 = \begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}x \end{array},
- *     \phi_5 = \begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}y \end{array},
- *     \phi_6 = \begin{array}{cc} x^2 \\ -2xy \end{array},
- *     \phi_7 = \begin{array}{cc} 2xy \\ -y^2 \end{array},
+ *     \phi_0 = \left[\begin{array}{cc} 1 \\ 0 \end{array}\right],
+ *     \phi_1 = \left[\begin{array}{cc} -\sqrt{3}+2\sqrt{3}x \\ 0 \end{array}\right],
+ *     \phi_2 = \left[\begin{array}{cc} -\sqrt{3}+2\sqrt{3}y \\ 0 \end{array}\right],
+ *     \phi_3 = \left[\begin{array}{cc} 0 \\ 1 \end{array}\right],
+ *     \phi_4 = \left[\begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}x \end{array}\right],
+ *     \phi_5 = \left[\begin{array}{cc} 0 \\ -\sqrt{3}+2\sqrt{3}y \end{array}\right],
+ *     \phi_6 = \left[\begin{array}{cc} x^2 \\ -2xy \end{array}\right],
+ *     \phi_7 = \left[\begin{array}{cc} 2xy \\ -y^2 \end{array}\right],
  *   @f}
  *
+ *   Thus, the dimension of the shape function space is dimension
+ *   times the number of polynomials of degree $k$ plus two:
+ *   @f[
+ *     n = 2\frac{(k+1)(k+2)}2 + 2
+ *   @f]
+
  *   <dt>In 3D:
  *   <dd> For any <i>i=0,...,k</i> the
  * curls of <i>(0,0,xy<sup>i+1</sup>z<sup>k-i</sup>)</i>,
  * <i>(x<sup>k-i</sup>yz<sup>i+1</sup>,0,0)</i> and
  * <i>(0,x<sup>i+1</sup>y<sup>k-i</sup>z,0)</i>
+ *
+ * The size of this function space is dimension times the number of
+ * polynomials of degree $k$ plus 3 times k+1:
+ *   @f[
+ *     n = 3\frac{(k+1)(k+2)(k+3)}6 + 3(k+1)
+ *   @f]
+ *
  * </dl>
  *
  * @todo Second derivatives in 3D are missing.
index 3e7061efc5e7564f812daebfb7568212889f4199..5ae1e3ee7fc9855f8a6d74409018655f30f32209 100644 (file)
@@ -34,7 +34,8 @@ DEAL_II_NAMESPACE_OPEN
  *
  * <h3>Degrees of freedom</h3>
  *
- * @todo This is for 2D only.
+ * @todo The 3D version exhibits some numerical instabilities, in
+ * particular for higher order
  *
  * @todo Restriction matrices are missing.
  *
@@ -102,13 +103,19 @@ private:
    * @ref GlossGeneralizedSupport "glossary entry on generalized support points"
    * for more information.
    */
-  void initialize_support_points (const unsigned int rt_degree);
+  void initialize_support_points (const unsigned int bdm_degree);
+  /**
+   * The values in the face support points of the polynomials needed as
+   * test functions. The outer vector is indexed by quadrature points, the
+   * inner by the test function. The test function space is PolynomialsP<dim-1>.
+   */
+  std::vector<std::vector<double> > test_values_face;
   /**
    * The values in the interior support points of the polynomials needed as
    * test functions. The outer vector is indexed by quadrature points, the
-   * inner by the test function.
+   * inner by the test function. The test function space is PolynomialsP<dim>.
    */
-  std::vector<std::vector<double> > test_values;
+  std::vector<std::vector<double> > test_values_cell;
 };
 
 DEAL_II_NAMESPACE_CLOSE
index 20fe8a291e2b631990522b4c37f1d429a0e94c16..203fbd5770a21d59b48c85478cc6a3de106ea14e 100644 (file)
@@ -296,11 +296,15 @@ namespace FETools
    *
    * @param isotropic_only Set to <code>true</code> if you only want to
    * compute matrices for isotropic refinement.
+   *
+   * @param threshold is the gap allowed in the least squares
+   * algorithm computing the embedding.
    */
   template <int dim, typename number, int spacedim>
   void compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
                                   std::vector<std::vector<FullMatrix<number> > > &matrices,
-                                  const bool isotropic_only = false);
+                                  const bool isotropic_only = false,
+                                  const double threshold = 1.e-12);
 
   /**
    * Compute the embedding matrices on faces needed for constraint matrices.
@@ -317,6 +321,9 @@ namespace FETools
    * @param face_fine The number of the face on the refined side of the face
    * for which this is computed.
    *
+   * @param threshold is the gap allowed in the least squares
+   * algorithm computing the embedding.
+   *
    * @warning This function will be used in computing constraint matrices. It
    * is not sufficiently tested yet.
    */
@@ -325,7 +332,8 @@ namespace FETools
   compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
                                   FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
                                   const unsigned int face_coarse,
-                                  const unsigned int face_fine);
+                                  const unsigned int face_fine,
+                                  const double threshold = 1.e-12);
 
   /**
    * For all possible (isotropic and anisotropic) refinement cases compute the
index d430f222298fc5d465da5bc4f2dcff17b1affd27..b1721a988a6f8d1bfa2f8dadd2bead33f84eba97 100644 (file)
@@ -34,6 +34,14 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
 }
 
 
+template <>
+unsigned int
+PolynomialSpace<0>::compute_n_pols (const unsigned int n)
+{
+  return 0;
+}
+
+
 template <>
 void
 PolynomialSpace<1>::
index f3904766c9ee777a1dc36b4f14422e871574ea8b..8b32c2ffaf8714d815ffbdcb95b6ff43b0c6fefd 100644 (file)
@@ -14,7 +14,9 @@
 // ---------------------------------------------------------------------
 
 
+#include <deal.II/base/geometry_info.h>
 #include <deal.II/base/polynomials_bdm.h>
+#include <deal.II/base/polynomial_space.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <iostream>
 #include <iomanip>
@@ -171,7 +173,7 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &unit_point,
               // p(t) = t^(i+1)
               monomials[i+1].value(unit_point(d), monovali[d]);
               // q(t) = t^(k-i)
-              monomials[degree()-i].value(unit_point(d), monovalk[d]);
+              monomials[degree()-i-1].value(unit_point(d), monovalk[d]);
             }
           if (values.size() != 0)
             {
index 7fb4e2fa3b28039c665b7691b580e5724fe7d90a..93542ad4176b0d615f50570f9c3152d65dcfa8cc 100644 (file)
@@ -45,7 +45,6 @@ FE_BDM<dim>::FE_BDM (const unsigned int deg)
                                std::vector<bool>(dim,true)))
 {
   Assert (dim >= 2, ExcImpossibleInDim(dim));
-  Assert (dim<3, ExcNotImplemented());
   Assert (deg > 0, ExcMessage("Lowest order BDM element are degree 1, but you asked for degree 0"));
 
   const unsigned int n_dofs = this->dofs_per_cell;
@@ -80,13 +79,15 @@ FE_BDM<dim>::FE_BDM (const unsigned int deg)
   // will be the correct ones, not
   // the raw shape functions anymore.
 
+  // Embedding errors become pretty large, so we just replace the
+  // regular threshold in both "computing_..." functions by 1.
   this->reinit_restriction_and_prolongation_matrices(true, true);
-  FETools::compute_embedding_matrices (*this, this->prolongation, true);
+  FETools::compute_embedding_matrices (*this, this->prolongation, true, 1.);
 
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i=0; i<GeometryInfo<dim>::max_children_per_face; ++i)
     face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
-  FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0);
+  FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.);
   this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face,
                                      this->dofs_per_face);
   unsigned int target_row=0;
@@ -168,44 +169,59 @@ FE_BDM<dim>::interpolate(
 
   // First do interpolation on faces. There, the component evaluated
   // depends on the face direction and orientation.
-  unsigned int fbase = 0;
-  unsigned int f=0;
-  for (; f<GeometryInfo<dim>::faces_per_cell;
-       ++f, fbase+=this->dofs_per_face)
+
+  // The index of the first dof on this face or the cell
+  unsigned int dbase = 0;
+  // The index of the first generalized support point on this face or the cell
+  unsigned int pbase = 0;
+  for (unsigned int f = 0; f<GeometryInfo<dim>::faces_per_cell; ++f)
     {
-      for (unsigned int i=0; i<this->dofs_per_face; ++i)
+      // Old version with no moments in 2D. See comment below in
+      // initialize_support_points()
+      if (test_values_face.size() == 0)
+        {
+          for (unsigned int i=0; i<this->dofs_per_face; ++i)
+            local_dofs[dbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][pbase+i];
+          pbase += this->dofs_per_face;
+        }
+      else
         {
-          local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i];
+          for (unsigned int i=0; i<this->dofs_per_face; ++i)
+            {
+              double s = 0.;
+              for (unsigned int k=0; k<test_values_face.size(); ++k)
+                s += values[GeometryInfo<dim>::unit_normal_direction[f]][pbase+k] * test_values_face[k][i];
+              local_dofs[dbase+i] = s;
+            }
+          pbase += test_values_face.size();
         }
+      dbase += this->dofs_per_face;
     }
 
+  AssertDimension (dbase, this->dofs_per_face * GeometryInfo<dim>::faces_per_cell);
+  AssertDimension (pbase, this->generalized_support_points.size() - test_values_cell.size());
+
   // Done for BDM1
-  if (fbase == this->dofs_per_cell) return;
+  if (dbase == this->dofs_per_cell) return;
 
   // What's missing are the interior
   // degrees of freedom. In each
   // point, we take all components of
   // the solution.
-  Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
-
-  // Here, the number of the point
-  // and of the shape function
-  // coincides. This will change
-  // below, since we have more
-  // support points than test
-  // functions in the interior.
-  const unsigned int pbase = fbase;
-  for (unsigned int d=0; d<dim; ++d, fbase += test_values[0].size())
+  Assert ((this->dofs_per_cell - dbase) % dim == 0, ExcInternalError());
+
+  for (unsigned int d=0; d<dim; ++d, dbase += test_values_cell[0].size())
     {
-      for (unsigned int i=0; i<test_values[0].size(); ++i)
+      for (unsigned int i=0; i<test_values_cell[0].size(); ++i)
         {
-          local_dofs[fbase+i] = 0.;
-          for (unsigned int k=0; k<test_values.size(); ++k)
-            local_dofs[fbase+i] += values[d][pbase+k] * test_values[k][i];
+          double s = 0.;
+          for (unsigned int k=0; k<test_values_cell.size(); ++k)
+            s += values[d][pbase+k] * test_values_cell[k][i];
+          local_dofs[dbase+i] = s;
         }
     }
 
-  Assert (fbase == this->dofs_per_cell, ExcInternalError());
+  Assert (dbase == this->dofs_per_cell, ExcInternalError());
 }
 
 
@@ -215,20 +231,17 @@ template <int dim>
 std::vector<unsigned int>
 FE_BDM<dim>::get_dpo_vector (const unsigned int deg)
 {
-  // the element is face-based and we have
-  // (deg+1)^(dim-1) DoFs per face
-  unsigned int dofs_per_face = 1;
-  for (unsigned int d=1; d<dim; ++d)
-    dofs_per_face *= deg+1;
-
-  // and then there are interior dofs
-  unsigned int
-  interior_dofs = dim==1?deg:dim*deg*(deg-1)/2;
-  if (dim>2)
-    {
-      interior_dofs *= deg-2;
-      interior_dofs /= 3;
-    }
+  // the element is face-based and we have as many degrees of freedom
+  // on the faces as there are polynomials of degree up to
+  // deg. Observe the odd convention of
+  // PolynomialSpace::compute_n_pols()!
+  unsigned int dofs_per_face = PolynomialSpace<dim-1>::compute_n_pols(deg+1);
+
+  // and then there are interior dofs, namely the number of
+  // polynomials up to degree deg-2 in dim dimensions.
+  unsigned int interior_dofs = 0;
+  if (deg>1)
+    interior_dofs = dim * PolynomialSpace<dim>::compute_n_pols(deg-1);
 
   std::vector<unsigned int> dpo(dim+1);
   dpo[dim-1] = dofs_per_face;
@@ -249,9 +262,12 @@ FE_BDM<dim>::get_ria_vector (const unsigned int deg)
       return std::vector<bool>();
     }
 
-  Assert(dim==2, ExcNotImplemented());
   const unsigned int dofs_per_cell = PolynomialsBDM<dim>::compute_n_pols(deg);
-  const unsigned int dofs_per_face = deg+1;
+  const unsigned int dofs_per_face = PolynomialSpace<dim-1>::compute_n_pols(deg);
+
+  Assert(GeometryInfo<dim>::faces_per_cell*dofs_per_face < dofs_per_cell,
+         ExcInternalError());
+
   // all dofs need to be
   // non-additive, since they have
   // continuity requirements.
@@ -266,94 +282,101 @@ FE_BDM<dim>::get_ria_vector (const unsigned int deg)
 }
 
 
+namespace
+{
+  // This function sets up the values of the polynomials we want to
+  // take moments with in the quadrature points. In fact, we multiply
+  // thos by the weights, such that the sum of function values and
+  // test_values over quadrature points yields the interpolated degree
+  // of freedom.
+  template <int dim>
+  void
+  initialize_test_values (std::vector<std::vector<double> > &test_values,
+                          const Quadrature<dim> &quadrature,
+                          const unsigned int deg)
+  {
+    PolynomialsP<dim> poly(deg);
+    std::vector<Tensor<1,dim> > dummy1;
+    std::vector<Tensor<2,dim> > dummy2;
+
+    test_values.resize(quadrature.size());
+
+    for (unsigned int k=0; k<quadrature.size(); ++k)
+      {
+        test_values[k].resize(poly.n());
+        poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2);
+        for (unsigned int i=0; i < poly.n(); ++i)
+          {
+            test_values[k][i] *= quadrature.weight(k);
+          }
+      }
+  }
+
+  // This specialization only serves to avoid error messages. Nothing
+  // useful can be computed in dimension zero and thus the vector
+  // length stays zero.
+  template <>
+  void
+  initialize_test_values (std::vector<std::vector<double> > &,
+                          const Quadrature<0> &,
+                          const unsigned int)
+  {}
+}
+
+
 template <int dim>
 void
 FE_BDM<dim>::initialize_support_points (const unsigned int deg)
 {
-  // interior point in 1d
-  unsigned int npoints = deg;
-  // interior point in 2d
-  if (dim >= 2)
-    {
-      npoints *= deg;
-//      npoints /= 2;
-    }
-  // interior point in 2d
-  if (dim >= 3)
-    {
-      npoints *= deg;
-//      npoints /= 3;
-    }
-  npoints += GeometryInfo<dim>::faces_per_cell * this->dofs_per_face;
+  // Our support points are quadrature points on faces and inside the
+  // cell. First on the faces, we have to test polynomials of degree
+  // up to deg, which means we need dg+1 points in each direction. The
+  // fact that we do not have tensor product polynomials will be
+  // considered later.
+  QGauss<dim-1> face_points (deg+1);
+
+  // Copy the quadrature formula to the face points.
+  this->generalized_face_support_points.resize (face_points.size());
+  for (unsigned int k=0; k<face_points.size(); ++k)
+    this->generalized_face_support_points[k] = face_points.point(k);
+
+  // In the interior, we only test with polynomials of degree up to
+  // deg-2, thus we use deg-1 points. Note that deg>=1 and the lowest
+  // order element has no points in the cell.
+  QGauss<dim> cell_points(deg-1);
+
+  // Compute the size of the whole support point set
+  const unsigned int npoints
+    = cell_points.size() + GeometryInfo<dim>::faces_per_cell * face_points.size();
 
   this->generalized_support_points.resize (npoints);
-  this->generalized_face_support_points.resize (this->dofs_per_face);
-
-  // Number of the point being entered
-  unsigned int current = 0;
-
-  // On the faces, we choose as many
-  // Gauss points as necessary to
-  // determine the normal component
-  // uniquely. This is the deg of
-  // the BDM element plus
-  // one.
-  if (dim>1)
-    {
-      QGauss<dim-1> face_points (deg+1);
-      Assert (face_points.size() == this->dofs_per_face,
-              ExcInternalError());
-      for (unsigned int k=0; k<this->dofs_per_face; ++k)
-        this->generalized_face_support_points[k] = face_points.point(k);
-      Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
-      for (unsigned int k=0;
-           k<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
-           ++k)
-        this->generalized_support_points[k] = faces.point(k+QProjector<dim>
-                                                          ::DataSetDescriptor::face(0,
-                                                              true,
-                                                              false,
-                                                              false,
-                                                              this->dofs_per_face));
-
-      current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
-    }
+
+  Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
+  for (unsigned int k=0; k < face_points.size()*GeometryInfo<dim>::faces_per_cell; ++k)
+    this->generalized_support_points[k]
+      = faces.point(k+QProjector<dim>
+                    ::DataSetDescriptor::face(0, true, false, false,
+                                              this->dofs_per_face));
+
+  // Currently, for backward compatibility, we do not use moments, but
+  // point values on faces in 2D. In 3D, this is impossible, since the
+  // moments are only taken with respect to PolynomialsP.
+  if (dim>2)
+    initialize_test_values(test_values_face, face_points, deg);
 
   if (deg<=1) return;
-  // Although the polynomial space is
-  // only P_{k-2}, we use the tensor
-  // product points for Q_{k-2}
-  QGauss<dim> quadrature(deg);
 
   // Remember where interior points start
-  const unsigned int ibase=current;
-//  for (unsigned int k=0;k<deg-1;++k)
-  for (unsigned int j=0; j<deg; ++j)
-    for (unsigned int i=0; i<deg; ++i)
-      {
-        this->generalized_support_points[current] = quadrature.point(current-ibase);
-        ++current;
-      }
-  Assert(current == npoints, ExcInternalError());
-
+  const unsigned int ibase = face_points.size()*GeometryInfo<dim>::faces_per_cell;
+  for (unsigned int k=0; k<cell_points.size(); ++k)
+    {
+      this->generalized_support_points[ibase+k] = cell_points.point(k);
+    }
   // Finally, compute the values of
   // the test functions in the
   // interior quadrature points
-  PolynomialsP<dim> poly(deg-2);
 
-  test_values.resize(quadrature.size());
-  std::vector<Tensor<1,dim> > dummy1;
-  std::vector<Tensor<2,dim> > dummy2;
-
-  for (unsigned int k=0; k<quadrature.size(); ++k)
-    {
-      test_values[k].resize(poly.n());
-      poly.compute(quadrature.point(k), test_values[k], dummy1, dummy2);
-      for (unsigned int i=0; i < poly.n(); ++i)
-        {
-          test_values[k][i] *= quadrature.weight(k);
-        }
-    }
+  initialize_test_values(test_values_cell, cell_points, deg-2);
 }
 
 
@@ -362,3 +385,4 @@ FE_BDM<dim>::initialize_support_points (const unsigned int deg)
 #include "fe_bdm.inst"
 
 DEAL_II_NAMESPACE_CLOSE
+
index 2510e6216b1f294d9db22e98b16d647323abc000..05004834dcf1f92861413cbb2ea4415d6a14c301 100644 (file)
@@ -34,6 +34,8 @@
 #include <deal.II/fe/fe_dgp_monomial.h>
 #include <deal.II/fe/fe_dgp_nonparametric.h>
 #include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_abf.h>
+#include <deal.II/fe/fe_bdm.h>
 #include <deal.II/fe/fe_raviart_thomas.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_system.h>
@@ -123,7 +125,9 @@ namespace
     result["FE_Q_Hierarchical"]
       = FEFactoryPointer(new FETools::FEFactory<FE_Q_Hierarchical<dim> >);
     result["FE_ABF"]
-      = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >);
+      = FEFactoryPointer(new FETools::FEFactory<FE_ABF<dim> >);
+    result["FE_BDM"]
+      = FEFactoryPointer(new FETools::FEFactory<FE_BDM<dim> >);
     result["FE_RaviartThomas"]
       = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >);
     result["FE_RaviartThomasNodal"]
@@ -708,7 +712,8 @@ namespace FETools
       const FiniteElement<dim, spacedim> &fe,
       const FEValues<dim, spacedim> &coarse,
       const Householder<double> &H,
-      FullMatrix<number> &this_matrix)
+      FullMatrix<number> &this_matrix,
+      const double threshold)
     {
       const unsigned int n  = fe.dofs_per_cell;
       const unsigned int nd = fe.n_components ();
@@ -741,8 +746,10 @@ namespace FETools
       // solve the least squares
       // problem.
       const double result = H.least_squares (v_fine, v_coarse);
+      Assert (result <= threshold, ExcLeastSquaresError (result));
+      // Avoid warnings in release mode
       (void)result;
-      Assert (result < 1.e-12, ExcLeastSquaresError (result));
+      (void)threshold;
 
       // Copy into the result
       // matrix. Since the matrix
@@ -760,7 +767,8 @@ namespace FETools
     compute_embedding_matrices_for_refinement_case (
       const FiniteElement<dim, spacedim> &fe,
       std::vector<FullMatrix<number> > &matrices,
-      const unsigned int ref_case)
+      const unsigned int ref_case,
+      const double threshold)
     {
       const unsigned int n  = fe.dofs_per_cell;
       const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
@@ -851,7 +859,7 @@ namespace FETools
                 {
                   task_group +=
                     Threads::new_task (&compute_embedding_for_shape_function<dim, number, spacedim>,
-                                       i, fe, coarse, H, this_matrix);
+                                       i, fe, coarse, H, this_matrix, threshold);
                 }
               task_group.join_all();
             }
@@ -860,7 +868,7 @@ namespace FETools
               for (unsigned int i = 0; i < n; ++i)
                 {
                   compute_embedding_for_shape_function<dim, number, spacedim>
-                  (i, fe, coarse, H, this_matrix);
+                  (i, fe, coarse, H, this_matrix, threshold);
                 }
             }
 
@@ -882,7 +890,8 @@ namespace FETools
   void
   compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
                              std::vector<std::vector<FullMatrix<number> > > &matrices,
-                             const bool isotropic_only)
+                             const bool isotropic_only,
+                             const double threshold)
   {
     Threads::TaskGroup<void> task_group;
 
@@ -893,7 +902,7 @@ namespace FETools
 
     for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
       task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>,
-                                       fe, matrices[ref_case-1], ref_case);
+                                       fe, matrices[ref_case-1], ref_case, threshold);
 
     task_group.join_all ();
   }
@@ -905,7 +914,8 @@ namespace FETools
   compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
                                   FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
                                   const unsigned int face_coarse,
-                                  const unsigned int face_fine)
+                                  const unsigned int face_fine,
+                                  const double threshold)
   {
     Assert(face_coarse==0, ExcNotImplemented());
     Assert(face_fine==0, ExcNotImplemented());
@@ -1077,8 +1087,10 @@ namespace FETools
             // solve the least squares
             // problem.
             const double result = H.least_squares(v_fine, v_coarse);
+            Assert (result <= threshold, ExcLeastSquaresError(result));
+            // Avoid compiler warnings in Release mode
             (void)result;
-            Assert (result < 1.e-12, ExcLeastSquaresError(result));
+            (void)threshold;
 
             // Copy into the result
             // matrix. Since the matrix
index 42291f8ecf4b9d58019c73993e52f49b200a3429..4ec919f04b6481865d0c909b9d06b4e9769cf71f 100644 (file)
@@ -32,7 +32,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       template
       void compute_embedding_matrices<deal_II_dimension, double, deal_II_space_dimension>
       (const FiniteElement<deal_II_dimension,deal_II_space_dimension> &,
-       std::vector<std::vector<FullMatrix<double> > > &,bool);
+       std::vector<std::vector<FullMatrix<double> > > &, const bool, const double);
 #endif
       \}
   }
@@ -111,7 +111,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       template
       void compute_face_embedding_matrices<deal_II_dimension,double>
       (const FiniteElement<deal_II_dimension> &, FullMatrix<double> ( &)[GeometryInfo<deal_II_dimension>::max_children_per_face],
-       unsigned int, unsigned int);
+       unsigned int, unsigned int, const double);
 
       template
       void
index 5c61614f92e68c0c8d83bf72055f404aa3c0c908..c7e8448bd2dc407d4221bae847ca2c58691e9fba 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2003 - 2014 by the deal.II authors
+// Copyright (C) 2003 - 2015 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -15,7 +15,7 @@
 
 
 
-// Show the shape functions of the Raviart-Thomas element on the unit cell
+// Show the shape functions of the BDM element on the unit cell
 // Plots are gnuplot compatible if lines with desired prefix are selected.
 
 #include "../tests.h"
@@ -91,9 +91,11 @@ main()
   deallog.threshold_double(1.e-10);
 
   for (unsigned int degree=1; degree<4; ++degree)
-    plot_shape_functions<2>(degree);
-//  plot_shape_functions<3>(degree);
-
+    {
+      plot_shape_functions<2>(degree);
+      plot_shape_functions<3>(degree);
+    }
+  
   return 0;
 }
 
index b50a69d664ecdcc6da866758ef787092d33188f7..d26b877186d15b2f5dc63446d5c795c6b9d89a5e 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-// Show the shape functions of the Raviart-Thomas element on a grid
+// Show the shape functions of the BDM element on a grid
 // with only one cell. This cell is rotated, stretched, scaled, etc,
 // and on each of these cells each time we evaluate the shape
 // functions.
@@ -172,7 +172,10 @@ main()
   deallog.threshold_double(1.e-10);
 
   for (unsigned int degree=1; degree<4; ++degree)
-    plot_shape_functions<2>(degree);
-
+    {
+      plot_shape_functions<2>(degree);
+      plot_shape_functions<3>(degree);
+    }
+  
   return 0;
 }
index f8ef4e9670c110288e5dba3c5a53c6781a254422..b8a9266d6a5cde473e6c27233dffc8c35c9107a5 100644 (file)
@@ -15,7 +15,7 @@
 
 
 
-// Just output the constraint matrices of the RT element
+// Just output the constraint matrices of the BDM element
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
@@ -60,8 +60,11 @@ main()
   deallog.threshold_double(1.e-10);
 
   for (unsigned int degree=1; degree<4; ++degree)
-    test<2>(degree);
-
+    {
+      test<2>(degree);
+      test<3>(degree);
+    }
+  
   return 0;
 }
 
index 64d5614302a3eca8471fbe006bbd08f2ec0de1f2..4a0121df09956b1023f2fa9475b157e4e37e189e 100644 (file)
@@ -66,7 +66,7 @@ main()
   for (unsigned int degree=1; degree<4; ++degree)
     {
       test<2>(degree);
-//      test<3>(degree);
+      test<3>(degree);
     }
 
   return 0;
index b0d72bbb7da2ae77653ddcdf62c1068c3950709f..48c046f239d7d4e3f9eb4371400b96b877e5eb79 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-// build a mass matrix for the RT element and try to invert it. we had trouble
+// build a mass matrix for the BDM element and try to invert it. we had trouble
 // with this at one time
 
 #include "../tests.h"
@@ -106,8 +106,11 @@ main()
   deallog.threshold_double(1.e-10);
 
   for (unsigned int i=1; i<4; ++i)
-    test<2>(i);
-
+    {
+      test<2>(i);
+      test<3>(i);
+    }
+  
   return 0;
 }
 
index 24397747bce3205564d6567851348893c4132eeb..050b40e782eeacbd65931dd31a4d99a350e92d18 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-// build a mass matrix for the RT element and try to invert it. like the rt_8
+// build a mass matrix for the BDM element and try to invert it. like the rt_8
 // test, except that we use a library function to build the mass matrix
 
 #include "../tests.h"
@@ -97,8 +97,11 @@ main()
   deallog.threshold_double(1.e-10);
 
   for (unsigned int i=1; i<4; ++i)
-    test<2>(i);
-
+    {
+      test<2>(i);
+      test<3>(i);
+    }
+  
   return 0;
 }
 
index 5dd65a10d980426e1c0c9f1c320ceef7d9262b11..77372b6289857cd10190234eb8ed71796f178919 100644 (file)
@@ -73,6 +73,13 @@ void test_2d_3d (std::vector<FiniteElement<dim> *> &fe_datas)
                                       FE_DGQ<dim> (1), 1));
   deallog << (*fe_datas.rbegin())->get_name() << std::endl;
 
+  FE_BDM<dim> *bdm1 = new FE_BDM<dim>(1);
+  fe_datas.push_back(bdm1);
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+  FE_BDM<dim> *bdm2 = new FE_BDM<dim>(2);
+  fe_datas.push_back(bdm2);
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+
                                   // Hcurl elements
   FE_Nedelec<dim> *ned0 = new FE_Nedelec<dim>(0);
   fe_datas.push_back(ned0);

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.