From 8754c480b22c0c66d3c8a29804de9c6bec1360f1 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Mon, 1 Feb 2016 11:26:01 +0100
Subject: [PATCH] remove invalid assert in compute_line_dof_identities()
 related to element domination; add a unit test

---
 source/hp/dof_handler.cc   |  13 ++---
 tests/fe/fe_nothing.cc     | 103 +++++++++++++++++++++++++++++++++++++
 tests/fe/fe_nothing.output |   9 ++++
 3 files changed, 116 insertions(+), 9 deletions(-)
 create mode 100644 tests/fe/fe_nothing.cc
 create mode 100644 tests/fe/fe_nothing.output

diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 5bc2f527d5..79b267ef3f 100644
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -2321,15 +2321,10 @@ namespace hp
 
                           if (i == (*finite_elements)[fe_index_1].dofs_per_line)
                             {
-                              // the dofs of these two finite elements are
-                              // identical. as a safety check, ensure that none
-                              // of the two FEs is trying to dominate the other,
-                              // which wouldn't make any sense in this case
-                              Assert ((*finite_elements)[fe_index_1].compare_for_face_domination
-                                      ((*finite_elements)[fe_index_2])
-                                      ==
-                                      FiniteElementDomination::either_element_can_dominate,
-                                      ExcInternalError());
+                              // The line dofs (i.e., the ones interior to a line) of these two finite elements are identical.
+                              // Note that there could be situations when one element still dominates another, e.g.:
+                              // FE_Q(2) x FE_Nothing(dominate) vs
+                              // FE_Q(2) x FE_Q(1)
 
                               --unique_sets_of_dofs;
 
diff --git a/tests/fe/fe_nothing.cc b/tests/fe/fe_nothing.cc
new file mode 100644
index 0000000000..43be9c2769
--- /dev/null
+++ b/tests/fe/fe_nothing.cc
@@ -0,0 +1,103 @@
+#include <deal.II/base/logstream.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include "../tests.h"
+
+#include <fstream>
+#include <iostream>
+
+const double eps = 1e-10;
+
+
+using namespace dealii;
+
+template<int dim>
+void test2cells(const unsigned int p1=2,
+                const unsigned int p2=1)
+{
+  deallog << "2cells: "<<dim<<"d, p1="<<p1<<", p2="<<p2<<std::endl;
+  Triangulation<dim>   triangulation;
+  {
+    Triangulation<dim>   triangulationL;
+    Triangulation<dim>   triangulationR;
+    GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain
+    GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain
+    Point<dim> shift_vector;
+    shift_vector[0] = 1.0;
+    GridTools::shift(shift_vector,triangulationR);
+    GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation);
+  }
+
+  hp::DoFHandler<dim> dof_handler(triangulation);
+
+  hp::FECollection<dim> fe_collection;
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p1),1,
+                                        FE_Q<dim>(p2),1));
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p1),1,
+                                        FE_Nothing<dim>(1,true),1));
+
+  // default is 0
+  dof_handler.begin_active()->set_active_fe_index(1);
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  ConstraintMatrix constraints;
+  constraints.clear();
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints.close ();
+
+  constraints.print(deallog.get_file_stream());
+
+  dof_handler.clear();
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(4);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+
+  try
+    {
+      test2cells<2>();
+      test2cells<3>();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/fe/fe_nothing.output b/tests/fe/fe_nothing.output
new file mode 100644
index 0000000000..d5ac6116e1
--- /dev/null
+++ b/tests/fe/fe_nothing.output
@@ -0,0 +1,9 @@
+
+DEAL::2cells: 2d, p1=2, p2=1
+    8 = 0
+    12 = 0
+DEAL::2cells: 3d, p1=2, p2=1
+    20 = 0
+    24 = 0
+    28 = 0
+    32 = 0
-- 
2.39.5