From 11c9f16691c70dd912550f46c2c524d01e41fe7b Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 28 Aug 2012 07:38:02 +0000
Subject: [PATCH] failing constraints

git-svn-id: https://svn.dealii.org/trunk@26138 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/fail/constraints_c1.cc | 157 +++++++++++++++++++++++++++++++++++
 1 file changed, 157 insertions(+)
 create mode 100644 tests/fail/constraints_c1.cc

diff --git a/tests/fail/constraints_c1.cc b/tests/fail/constraints_c1.cc
new file mode 100644
index 0000000000..dcec5fdb84
--- /dev/null
+++ b/tests/fail/constraints_c1.cc
@@ -0,0 +1,157 @@
+/**********************************************************************
+ * $Id: constraint_test.cc -1   $
+ *
+ * Copyright Guido Kanschat, 2010, 2012
+ *
+ **********************************************************************/
+
+// H2-elements implemented through constraints on the degrees of
+// freedom. After adding all constrained lines, the program hangs in
+// constraints.close().
+
+
+#include <deal.II/base/job_identifier.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+
+using namespace dealii;
+
+template <int dim>
+void
+setup_constraints(const DoFHandler<dim>& dof_handler)
+{
+  ConstraintMatrix     constraints;
+  constraints.clear();
+  const FiniteElement<dim>& fe = dof_handler.get_fe();
+  
+				   // Set up derivative constraints to
+				   // make element C1
+  std::vector<std::vector<std::vector<double> > > vertex_constraints(
+    GeometryInfo<dim>::vertices_per_cell, std::vector<std::vector<double> >(
+      dim+1, std::vector<double>(fe.dofs_per_cell, 0.)));
+  
+  for (unsigned int vertex=0;vertex<GeometryInfo<dim>::vertices_per_cell;++vertex)
+    {
+      Point<dim> v;
+      for (unsigned int d=0;d<dim;++d)
+	{
+	  const unsigned int ds = 1<<d;
+	  const unsigned int id = vertex/ds;
+	  v[d] = (id%2 != 0) ? 1. : 0.;
+	}
+      
+      for (unsigned int j=0;j<fe.dofs_per_cell;++j)
+	{
+	  for (unsigned int d=0;d<dim;++d)
+	    {
+	      const double f = fe.shape_grad(j,v)[d];
+	      if (std::fabs(f) > 1.e-12)
+		vertex_constraints[vertex][d][j] = f;
+	    }
+	  const double f = fe.shape_grad_grad(j,v)[0][1];
+	  if (std::fabs(f) > 1.e-12)
+	    vertex_constraints[vertex][dim][j] = f;
+	}      
+    }
+  
+  std::vector<unsigned int> cell_indices(fe.dofs_per_cell);
+  std::vector<unsigned int> neighbor_indices(fe.dofs_per_cell);
+  
+  for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+       cell != dof_handler.end();++cell)
+    {
+      deallog << "New cell" << std::endl;
+      
+      cell->get_dof_indices(cell_indices);
+      
+				       // Do lower left and upper
+				       // right vertex
+      for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;face+=2)
+      {
+	if (cell->at_boundary(face)) continue;
+	for (unsigned int fvertex=0;fvertex<GeometryInfo<dim>::vertices_per_face;++fvertex)
+	    {
+	      const unsigned int other_face = (face==0) ? (2+fvertex) : fvertex;
+	      if (cell->at_boundary(other_face)) continue;
+	      
+	      const unsigned int vertex = GeometryInfo<dim>::face_to_cell_vertices(face, fvertex);
+	      typename DoFHandler<dim>::active_cell_iterator neighbor = cell->neighbor(other_face);
+	      const unsigned int neighbor_face = cell->neighbor_of_neighbor(other_face);
+	      neighbor->get_dof_indices(neighbor_indices);
+	      unsigned int neighbor_vertex = GeometryInfo<dim>::face_to_cell_vertices(neighbor_face, fvertex);
+	      
+	      const unsigned int d = GeometryInfo<dim>::unit_normal_direction[other_face];
+	      
+	      std::vector<std::pair<unsigned int, double> > rhs;
+	      const unsigned int constrained = fe.face_to_cell_index(
+		(fvertex == 0) ? 2 : fe.dofs_per_face-1, face);
+	      double constrained_weight = 0.;
+	      
+	      for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+		{
+		  if (i == constrained)
+		    constrained_weight = vertex_constraints[vertex][d][i];
+		  else
+		    if (vertex_constraints[vertex][d][i] != 0.)
+		      rhs.push_back(std::make_pair(cell_indices[i], -vertex_constraints[vertex][d][i]));
+		  if (vertex_constraints[neighbor_vertex][d][i] != 0.)
+		    rhs.push_back(std::make_pair(neighbor_indices[i], vertex_constraints[neighbor_vertex][d][i]));
+		}
+	      deallog << " v" << vertex
+		      << " f" << face
+		      << " o" << other_face
+		      << " d" << d
+		      << " l" << constrained
+		      << " g" << cell_indices[constrained]
+//		      << " w " << constrained_weight
+		      << " rhs ";
+	      for (unsigned int i=0;i<rhs.size();++i)
+		deallog << ' ' << rhs[i].first;
+	      deallog << std::endl;
+	      for (unsigned int i=0;i<rhs.size();++i)
+		rhs[i].second /= constrained_weight;
+	      constraints.add_line(cell_indices[constrained]);
+	      constraints.add_entries(cell_indices[constrained], rhs);
+	    }
+	}
+    }
+
+  deallog << "Closing" << std::endl;
+  constraints.close();
+  deallog << "Closed" << std::endl;
+}
+
+
+template <int dim>
+void
+run(const FiniteElement<dim>& fe)
+{
+  deallog << "Element: " << fe.get_name() << std::endl;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global(2);
+  deallog << "Triangulation "
+	  << triangulation.n_active_cells() << " cells, "
+	  << triangulation.n_levels() << " levels" << std::endl;
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+  setup_constraints (dof_handler);  
+}
+
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  
+  FE_Q<2> fe(3);
+  run(fe);
+}
-- 
2.39.5