FE_Nedelec can now be constructed for up to degree 12 in debug mode

Set a degree-dependent threshold for computing face embedding matrices.
Fixes #2176.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 5f0e39a585373a2dc0aa95b1ed88cc73b86cb1f9..65c4fed64e8edc1800f093274e35422c4822030f 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -57,6 +57,12 @@ inconvenience this causes.
 
 
 <ol> 
+ <li> Fixed: FE_Nedelec elements up to polynomial order 12 can now be
+ constructed.
+ <br>
+ (Jean-Paul Pelteret, 2016/02/12)
+ </li>
+
  <li> Fixed: The GridTools::build_triangulation_from_patches() function now 
  also copies the locations of vertices from the cells of the source 
  triangulation to the triangulation that is built from the list of patch cells.

diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index d9011e17d0965fd4d6c935a6bdef17a46d56b5e9..2954d25e5d3b6745a83f0bda31fba0406ff6d2d7 100644 (file)
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -40,6 +40,23 @@ DEAL_II_NAMESPACE_OPEN
 
 //#define DEBUG_NEDELEC
 
+namespace internal
+{
+  namespace
+  {
+    double
+    get_embedding_computation_tolerance(const unsigned int p)
+    {
+      // This heuristic was computed by monitoring the worst residual
+      // resulting from the least squares computation when computing
+      // the face embedding matrices in the FE_Nedelec constructor.
+      // The residual growth is exponential, but is bounded by this
+      // function up to degree 12.
+      return 1.e-15*std::exp(std::pow(p,1.075));
+    }
+  }
+}
+
 
 template <int dim>
 FE_Nedelec<dim>::FE_Nedelec (const unsigned int p) :
@@ -86,7 +103,8 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int p) :
     face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
 
   FETools::compute_face_embedding_matrices<dim,double>
-  (*this, face_embeddings, 0, 0);
+  (*this, face_embeddings, 0, 0,
+   internal::get_embedding_computation_tolerance(p));
 
   switch (dim)
     {
@@ -3002,7 +3020,8 @@ FE_Nedelec<dim>
 #endif
       this_nonconst.reinit_restriction_and_prolongation_matrices ();
       // Fill prolongation matrices with embedding operators
-      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true);
+      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
+                                           internal::get_embedding_computation_tolerance(this->degree));
 #ifdef DEBUG_NEDELEC
       deallog << "Restriction" << std::endl;
 #endif
@@ -3052,7 +3071,8 @@ FE_Nedelec<dim>
 #endif
       this_nonconst.reinit_restriction_and_prolongation_matrices ();
       // Fill prolongation matrices with embedding operators
-      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true);
+      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
+                                           internal::get_embedding_computation_tolerance(this->degree));
 #ifdef DEBUG_NEDELEC
       deallog << "Restriction" << std::endl;
 #endif