From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Fri, 12 Feb 2016 22:20:11 +0000 (+0100)
Subject: FE_Nedelec can now be constructed for up to degree 12 in debug mode
X-Git-Tag: v8.5.0-rc1~1318^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f8d6a470b9b44645dd419050a67c0c913cfd0b1f;p=dealii.git

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 5f0e39a585..65c4fed64e 100644
--- 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 d9011e17d0..2954d25e5d 100644
--- 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