From 3b37071b97f2e88d4fab627dd9f34daf2a00da75 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 31 May 2005 14:26:15 +0000
Subject: [PATCH] New test

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

diff --git a/tests/bits/rt_2.cc b/tests/bits/rt_2.cc
new file mode 100644
index 0000000000..68b8b5a9d5
--- /dev/null
+++ b/tests/bits/rt_2.cc
@@ -0,0 +1,164 @@
+//----------------------------  rt_1.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 by the deal.II authors and Oliver Kayser-Herold
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  rt_1.cc  ---------------------------
+
+// Small test to analyse the equivalence from both sides of a face of
+// the normal component for the Raviart-Thomas element
+//
+// Test adapted from a small program by Oliver Kayser-Herold.
+
+#include "../tests.h"
+#include <fstream>
+
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+#include <fe/mapping_q.h>
+#include <fe/mapping_q1_eulerian.h>
+
+#include <numerics/vectors.h>
+#include <numerics/data_out.h>
+#include <lac/vector.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_dgq.h>
+
+
+void evaluate_normal (DoFHandler<2>  &dof_handler,
+		      Vector<double> &solution)
+{
+				   // This quadrature rule determines
+				   // the points, where the continuity
+				   // will be tested.
+  QGauss<1> quad (6);
+  FEFaceValues<2> fe_v_face (dof_handler.get_fe (), quad, 
+			     UpdateFlags(update_values    |
+					 update_q_points  |
+					 update_gradients |
+					 update_normal_vectors |
+					 update_JxW_values));
+
+  FEFaceValues<2> fe_v_face_n (dof_handler.get_fe (), quad, 
+			       UpdateFlags(update_values    |
+					   update_q_points  |
+					   update_gradients |
+					   update_normal_vectors |
+					   update_JxW_values));
+
+  const unsigned int   n_q_face    = quad.n_quadrature_points;
+  const unsigned int   n_components   = dof_handler.get_fe().n_components();
+
+				   // Cell iterators
+  DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
+				      endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+	{
+	  if (!cell->face(f)->at_boundary ())
+	    {
+	      deallog << "cell=" << cell
+		      << ", face=" << f
+		      << std::endl;
+	      
+	      fe_v_face.reinit (cell, f);
+		
+	      const unsigned int neighbor=cell->neighbor_of_neighbor(f);
+	      fe_v_face_n.reinit (cell->neighbor (f), neighbor); 
+
+					       // Get values from
+					       // solution vector (For
+					       // Trap.Rule)
+	      std::vector<Vector<double> > this_value
+		(n_q_face, Vector<double>(n_components));
+	      fe_v_face.get_function_values (solution, this_value);
+
+					       // Same for neighbor cell
+	      std::vector<Vector<double> > this_value_n
+		(n_q_face, Vector<double>(n_components));
+	      fe_v_face_n.get_function_values (solution, this_value_n);
+
+	      for (unsigned int q_point=0; q_point<n_q_face; ++q_point)
+		{
+		  Point<2> vn = fe_v_face.normal_vector (q_point);
+		  double nx = vn(0);
+		  double ny = vn(1);
+		    
+		  double u = this_value[q_point](0);
+		  double v = this_value[q_point](1);
+
+		  double u_n = this_value_n[q_point](0);
+		  double v_n = this_value_n[q_point](1);
+
+		  deallog << "quadrature point " << q_point
+			  << " (u-u_n)*nx+(v-v_n)*ny="
+			  << (u-u_n) * nx + (v-v_n) * ny
+			  << ", u*nx+v*ny="
+			  << u * nx + v * ny
+			  << ", u_n*nx+v_n*ny="
+			  << u_n * nx + v_n * ny
+			  << std::endl;
+		}
+	    }
+	}
+    }
+}
+
+
+int main ()
+{
+  std::ofstream logfile("rt_2.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Triangulation<2> tria_test;
+  Point<2> p1 (0,0),
+    p2 (1, 1);
+  std::vector<unsigned int> sub_div;
+
+  sub_div.push_back (1);
+  sub_div.push_back (4);
+
+  GridGenerator::subdivided_hyper_rectangle (tria_test, sub_div, p1, p2);
+  tria_test.refine_global (2);
+
+  FE_RaviartThomas<2> fe (1);
+  DoFHandler<2> dof_handler(tria_test);
+  dof_handler.distribute_dofs (fe);
+
+				   // Fill solution vector with random
+				   // values between 0 and 1.
+  Vector<double> solution(dof_handler.n_dofs ());
+  for (unsigned int i = 0; i < dof_handler.n_dofs (); ++i)
+    solution(i) = (double) rand () / (double) RAND_MAX;
+
+				   // Now check if the function is
+				   // continuous in normal direction
+				   // on uniform mesh:
+  deallog << "Uniform mesh test" << std::endl;
+  evaluate_normal (dof_handler, solution);
+
+
+				   // Then test same on distorted mesh
+  tria_test.distort_random (0.05);
+  deallog << "Distorted mesh test" << std::endl;
+  evaluate_normal (dof_handler, solution);
+}
-- 
2.39.5