From fd26695f99680c18cd22a99f325742571d69fd2e Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Tue, 1 Sep 2015 20:52:44 +0200
Subject: [PATCH] test make_hanging_node_constraints() with p-ref and
 neither_element_dominates

---
 .../hp/hp_constraints_neither_dominate_01.cc  | 260 ++++++++++++++++++
 .../hp_constraints_neither_dominate_01.output | 101 +++++++
 2 files changed, 361 insertions(+)
 create mode 100644 tests/hp/hp_constraints_neither_dominate_01.cc
 create mode 100644 tests/hp/hp_constraints_neither_dominate_01.output

diff --git a/tests/hp/hp_constraints_neither_dominate_01.cc b/tests/hp/hp_constraints_neither_dominate_01.cc
new file mode 100644
index 0000000000..36824cecce
--- /dev/null
+++ b/tests/hp/hp_constraints_neither_dominate_01.cc
@@ -0,0 +1,260 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test hp_constraints when neither_element_dominates.
+// 2D and 3D p-refinement only for 2 elements.
+// Check continuity across the face by evaluating the
+// field from both sides.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_parser.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/hp/fe_collection.h>
+
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.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/fe/fe_values.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+
+#include <deal.II/lac/slepc_solver.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+#include <fstream>
+#include <iostream>
+
+// #define DEBUG_OUTPUT_VTK
+
+const double eps = 1e-10;
+
+using namespace dealii;
+
+template<int dim>
+void test2cells(const unsigned int p1=2,
+                const unsigned int p2=1)
+{
+  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_Nothing<dim>(1,true),1));
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p2),1,
+                                        FE_Q<dim>(1),1));
+  // push back to be able to resolve hp constrains no matter what:
+  fe_collection.push_back(FESystem<dim>(FE_Q<dim>(p2),1,
+                                        FE_Nothing<dim>(1,true),1));
+
+  hp::QCollection< dim-1 > q_face_collection;
+  q_face_collection.push_back(QIterated<dim-1>(QTrapez<1>(),2));
+  q_face_collection.push_back(QIterated<dim-1>(QTrapez<1>(),2));
+  q_face_collection.push_back(QIterated<dim-1>(QTrapez<1>(),2));
+
+  deallog << "2cells: "<<fe_collection[0].get_name()<< " vs "<<fe_collection[1].get_name()<<std::endl;
+
+  dof_handler.begin_active()->set_active_fe_index(1);//Q(p2)xQ(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());
+
+#ifdef DEBUG_OUTPUT_VTK
+  // output to check if all is good:
+  std::vector<Vector<double>> shape_functions;
+  std::vector<std::string> names;
+  for (unsigned int s=0; s < dof_handler.n_dofs(); s++)
+    {
+      Vector<double> shape_function;
+      shape_function.reinit(dof_handler.n_dofs());
+      shape_function[s] = 1.0;
+
+      // if the dof is constrained, first output unconstrained vector
+      if (constraints.is_constrained(s))
+        {
+          names.push_back(std::string("UN_") +
+                          dealii::Utilities::int_to_string(s,2));
+          shape_functions.push_back(shape_function);
+        }
+
+      names.push_back(std::string("N_") +
+                      dealii::Utilities::int_to_string(s,2));
+
+      // make continuous/constrain:
+      constraints.distribute(shape_function);
+      shape_functions.push_back(shape_function);
+    }
+
+  DataOut<dim,hp::DoFHandler<dim>> data_out;
+  data_out.attach_dof_handler (dof_handler);
+
+  // get material ids:
+  Vector<float> fe_index(triangulation.n_active_cells());
+  typename hp::DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (unsigned int index=0; cell!=endc; ++cell,++index)
+    {
+      fe_index[index] = cell->active_fe_index();
+    }
+  data_out.add_data_vector(fe_index, "fe_index");
+
+  for (unsigned int i = 0; i < shape_functions.size(); i++)
+    data_out.add_data_vector (shape_functions[i], names[i]);
+
+  data_out.build_patches(15);
+
+  std::string filename = "2cell-shape_functions_p1="+std::to_string(p1)
+    +"_p2="+std::to_string(p2)
+    +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+  std::ofstream output (filename.c_str ());
+  data_out.write_vtu (output);
+#endif
+
+  // fill some vector
+  Vector<double> solution(dof_handler.n_dofs());
+  for (unsigned int dof=0; dof < dof_handler.n_dofs(); dof++)
+    solution[dof] = 21.0*(dof%2) + 0.5 + dof%3;
+
+  constraints.distribute(solution);
+
+
+  // evaluate field at the interface:
+  hp::FEFaceValues< dim > fe_face_values_hp(fe_collection,
+                                            q_face_collection,
+                                            update_values | update_quadrature_points);
+
+  std::vector<unsigned int> local_face_dof_indices;
+  std::vector<Vector<double> > values;
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+      {
+        const unsigned int fe_index = cell->active_fe_index();
+        local_face_dof_indices.resize(fe_collection[fe_index].dofs_per_face);
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (!cell->at_boundary(f)) // that's enough for our simple mesh
+            {
+              deallog << "cell="<<cell<<" face="<<f<<std::endl;
+              fe_face_values_hp.reinit(cell,f);
+              const FEFaceValues<dim> &fe_face_values = fe_face_values_hp.get_present_fe_values();
+
+              const unsigned int n_q_points = fe_face_values.n_quadrature_points;
+              values.resize (n_q_points,
+                             Vector<double>(2));
+
+              fe_face_values.get_function_values (solution,values);
+
+              const std::vector<dealii::Point<dim> > &q_points = fe_face_values.get_quadrature_points();
+
+              for (unsigned int q = 0; q < n_q_points; q++)
+                deallog << "u["<<q_points[q]<<"]={"<<values[q][0]<<","<<values[q][1]<<"}"<<std::endl;
+            }
+      }
+
+  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>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
+    test2cells<2>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
+    test2cells<3>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
+    test2cells<3>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
+  }
+  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/hp/hp_constraints_neither_dominate_01.output b/tests/hp/hp_constraints_neither_dominate_01.output
new file mode 100644
index 0000000000..10f1da2b1b
--- /dev/null
+++ b/tests/hp/hp_constraints_neither_dominate_01.output
@@ -0,0 +1,101 @@
+
+DEAL::2cells: FESystem<2>[FE_Q<2>(1)-FE_Nothing<2>()] vs FESystem<2>[FE_Q<2>(2)-FE_Q<2>(1)]
+    2 = 0
+    5 = 0
+    7 11:  0.500000
+    7 13:  0.500000
+DEAL::cell=0.0 face=1
+DEAL::u[0.0000 -1.0000]={23.5000,0}
+DEAL::u[0.0000 -0.5000]={23.0000,0}
+DEAL::u[0.0000 0.0000]={22.5000,0}
+DEAL::cell=0.1 face=0
+DEAL::u[0.0000 -1.0000]={23.5000,0}
+DEAL::u[0.0000 -0.5000]={23.0000,0}
+DEAL::u[0.0000 0.0000]={22.5000,0}
+DEAL::2cells: FESystem<2>[FE_Q<2>(2)-FE_Nothing<2>()] vs FESystem<2>[FE_Q<2>(1)-FE_Q<2>(1)]
+    2 = 0
+    5 = 0
+    10 6:  0.500000
+    10 8:  0.500000
+DEAL::cell=0.0 face=1
+DEAL::u[0.0000 -1.0000]={0.5000,0}
+DEAL::u[0.0000 -0.5000]={1.5000,0}
+DEAL::u[0.0000 0.0000]={2.5000,0}
+DEAL::cell=0.1 face=0
+DEAL::u[0.0000 -1.0000]={0.5000,0}
+DEAL::u[0.0000 -0.5000]={1.5000,0}
+DEAL::u[0.0000 0.0000]={2.5000,0}
+DEAL::2cells: FESystem<3>[FE_Q<3>(1)-FE_Nothing<3>()] vs FESystem<3>[FE_Q<3>(2)-FE_Q<3>(1)]
+    2 = 0
+    5 = 0
+    8 = 0
+    11 = 0
+    13 31:  0.500000
+    13 33:  0.500000
+    17 35:  0.500000
+    17 37:  0.500000
+    21 31:  0.500000
+    21 35:  0.500000
+    23 33:  0.500000
+    23 37:  0.500000
+    25 31:  0.250000
+    25 33:  0.250000
+    25 35:  0.250000
+    25 37:  0.250000
+DEAL::cell=0.0 face=1
+DEAL::u[0.0000 -1.0000 -1.0000]={22.5000,0}
+DEAL::u[0.0000 -0.5000 -1.0000]={22.0000,0}
+DEAL::u[0.0000 0.0000 -1.0000]={21.5000,0}
+DEAL::u[0.0000 -1.0000 -0.5000]={23.0000,0}
+DEAL::u[0.0000 -0.5000 -0.5000]={22.5000,0}
+DEAL::u[0.0000 0.0000 -0.5000]={22.0000,0}
+DEAL::u[0.0000 -1.0000 0.0000]={23.5000,0}
+DEAL::u[0.0000 -0.5000 0.0000]={23.0000,0}
+DEAL::u[0.0000 0.0000 0.0000]={22.5000,0}
+DEAL::cell=0.1 face=0
+DEAL::u[0.0000 -1.0000 -1.0000]={22.5000,0}
+DEAL::u[0.0000 -0.5000 -1.0000]={22.0000,0}
+DEAL::u[0.0000 0.0000 -1.0000]={21.5000,0}
+DEAL::u[0.0000 -1.0000 -0.5000]={23.0000,0}
+DEAL::u[0.0000 -0.5000 -0.5000]={22.5000,0}
+DEAL::u[0.0000 0.0000 -0.5000]={22.0000,0}
+DEAL::u[0.0000 -1.0000 0.0000]={23.5000,0}
+DEAL::u[0.0000 -0.5000 0.0000]={23.0000,0}
+DEAL::u[0.0000 0.0000 0.0000]={22.5000,0}
+DEAL::2cells: FESystem<3>[FE_Q<3>(2)-FE_Nothing<3>()] vs FESystem<3>[FE_Q<3>(1)-FE_Q<3>(1)]
+    2 = 0
+    5 = 0
+    8 = 0
+    11 = 0
+    20 12:  0.500000
+    20 14:  0.500000
+    24 16:  0.500000
+    24 18:  0.500000
+    28 12:  0.500000
+    28 16:  0.500000
+    30 14:  0.500000
+    30 18:  0.500000
+    32 12:  0.250000
+    32 14:  0.250000
+    32 16:  0.250000
+    32 18:  0.250000
+DEAL::cell=0.0 face=1
+DEAL::u[0.0000 -1.0000 -1.0000]={0.5000,0}
+DEAL::u[0.0000 -0.5000 -1.0000]={1.5000,0}
+DEAL::u[0.0000 0.0000 -1.0000]={2.5000,0}
+DEAL::u[0.0000 -1.0000 -0.5000]={1.0000,0}
+DEAL::u[0.0000 -0.5000 -0.5000]={1.2500,0}
+DEAL::u[0.0000 0.0000 -0.5000]={1.5000,0}
+DEAL::u[0.0000 -1.0000 0.0000]={1.5000,0}
+DEAL::u[0.0000 -0.5000 0.0000]={1.0000,0}
+DEAL::u[0.0000 0.0000 0.0000]={0.5000,0}
+DEAL::cell=0.1 face=0
+DEAL::u[0.0000 -1.0000 -1.0000]={0.5000,0}
+DEAL::u[0.0000 -0.5000 -1.0000]={1.5000,0}
+DEAL::u[0.0000 0.0000 -1.0000]={2.5000,0}
+DEAL::u[0.0000 -1.0000 -0.5000]={1.0000,0}
+DEAL::u[0.0000 -0.5000 -0.5000]={1.2500,0}
+DEAL::u[0.0000 0.0000 -0.5000]={1.5000,0}
+DEAL::u[0.0000 -1.0000 0.0000]={1.5000,0}
+DEAL::u[0.0000 -0.5000 0.0000]={1.0000,0}
+DEAL::u[0.0000 0.0000 0.0000]={0.5000,0}
-- 
2.39.5