From 0bb3408d3fc1948eba2874933d578f51f9249fbf Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 21 Apr 2003 16:23:52 +0000
Subject: [PATCH] Tests for Raviart-Thomas elements.

git-svn-id: https://svn.dealii.org/trunk@7427 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/fe/rt_1.cc |  76 ++++++++++++++++++++++
 tests/fe/rt_2.cc | 161 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 237 insertions(+)
 create mode 100644 tests/fe/rt_1.cc
 create mode 100644 tests/fe/rt_2.cc

diff --git a/tests/fe/rt_1.cc b/tests/fe/rt_1.cc
new file mode 100644
index 0000000000..1911800c15
--- /dev/null
+++ b/tests/fe/rt_1.cc
@@ -0,0 +1,76 @@
+// $Id$
+// (c) Wolfgang Bangerth
+//
+// Show the shape functions of the Raviart-Thomas element on the unit cell
+
+#include <base/logstream.h>
+#include <fe/fe_raviart_thomas.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 2
+
+
+
+template<int dim>
+inline void
+plot_shape_functions(const unsigned int degree)
+{
+  FE_RaviartThomas<dim> fe_rt(degree);
+
+  const unsigned int div=2;
+  for (unsigned int mz=0;mz<=((dim>2) ? div : 0) ;++mz)
+    for (unsigned int my=0;my<=((dim>1) ? div : 0) ;++my)
+      for (unsigned int mx=0;mx<=div;++mx)
+        {
+          const Point<dim> p = (dim==2 ?
+                                Point<dim>(1.*mx/div, 1.*my/div) :
+                                Point<dim>(1.*mx/div, 1.*my/div, 1.*mz/div));
+          deallog << "q_point = " << p << std::endl;
+	      
+          for (unsigned int i=0;i<fe_rt.dofs_per_cell;++i)
+            {
+              deallog << "  phi(" << i << ") = [";
+              for (unsigned int c=0; c<dim; ++c)
+                deallog << " " << fe_rt.shape_value_component(i,p,c);
+              deallog << " ]" << std::endl;
+            };
+          for (unsigned int i=0;i<fe_rt.dofs_per_cell;++i)
+            {
+              deallog << "  grad phi(" << i << ") = ";
+              for (unsigned int c=0; c<dim; ++c)
+                {
+                  deallog << "[";
+                  for (unsigned int d=0; d<dim; ++d)
+                    deallog << " " << fe_rt.shape_grad_component(i,p,c)[d];
+                  deallog << " ]" << std::endl;
+                  if (c != dim-1)
+                    deallog << "                ";
+                };
+            };
+        }
+  
+  deallog << std::endl;
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("rt_1.output");
+  logfile.precision (PRECISION);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  for (unsigned int degree=0; degree<2; ++degree)
+    plot_shape_functions<2>(degree);
+//  plot_shape_functions<3>(degree);
+  
+  return 0;
+}
+
+
+
diff --git a/tests/fe/rt_2.cc b/tests/fe/rt_2.cc
new file mode 100644
index 0000000000..1cb6e336db
--- /dev/null
+++ b/tests/fe/rt_2.cc
@@ -0,0 +1,161 @@
+// $Id$
+// (c) Wolfgang Bangerth
+//
+// Show the shape functions of the Raviart-Thomas element on a grid
+// with only one cell. This cell is rotated, stretched, scaled, etc,
+// and on each of these cells each time we evaluate the shape
+// functions.
+
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_values.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 2
+
+
+Point<2> stretch_coordinates (const Point<2> p)
+{
+  return Point<2>(2*p(0), p(1));
+}
+
+
+
+Point<2> tilt_coordinates (const Point<2> p)
+{
+  return Point<2>(p(0)+p(1), p(1));
+}
+
+
+void
+transform_grid (Triangulation<2> &tria,
+		const unsigned int  transform)
+{
+  switch (transform)
+    {
+                                       // first round: take
+                                       // original grid
+      case 0:
+	    break;
+            
+					     // second round: rotate
+					     // triangulation
+      case 1:
+	    GridTools::rotate (3.14159265358/2, tria);
+	    break;
+
+					     // third round: inflate
+					     // by a factor of 2
+      case 2:
+	    GridTools::scale (2, tria);
+	    break;
+
+                                             // third round: scale
+                                             // back, rotate back,
+                                             // stretch
+      case 3:
+            GridTools::scale (.5, tria);
+	    GridTools::rotate (-3.14159265358/2, tria);
+            GridTools::transform (&stretch_coordinates, tria);
+            
+	    break;
+            
+      default:
+	    Assert (false, ExcNotImplemented());
+    };
+}
+
+	    
+
+
+template<int dim>
+inline void
+plot_shape_functions()
+{
+  FE_RaviartThomas<dim> fe_rt(0);
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr, 0., 1.);
+
+				   // check the following with a
+				   // number of transformed
+				   // triangulations
+  for (unsigned int transform=0; transform<4; ++transform)
+    {
+      deallog << "GRID TRANSFORMATION " << transform << std::endl;
+      
+      transform_grid (tr, transform);
+
+      DoFHandler<dim> dof(tr);
+      typename DoFHandler<dim>::cell_iterator c = dof.begin();
+      dof.distribute_dofs(fe_rt);
+      
+      QTrapez<1> q_trapez;
+      const unsigned int div=2;
+      QIterated<dim> q(q_trapez, div);
+      FEValues<dim> fe(fe_rt, q, update_values|update_gradients|update_q_points);
+      fe.reinit(c);
+      
+      unsigned int q_point=0;
+      for (unsigned int mz=0;mz<=((dim>2) ? div : 0) ;++mz)
+	for (unsigned int my=0;my<=((dim>1) ? div : 0) ;++my)
+	  for (unsigned int mx=0;mx<=div;++mx)
+	    {
+	      deallog << "q_point(" << q_point << ")=" << fe.quadrature_point(q_point)
+		      << std::endl;
+	      
+	      for (unsigned int i=0;i<fe_rt.dofs_per_cell;++i)
+		{
+		  deallog << "  phi(" << i << ") = [";
+		  for (unsigned int c=0; c<fe.get_fe().n_components(); ++c)
+		    deallog << " " << fe.shape_value_component(i,q_point,c);
+		  deallog << " ]" << std::endl;
+		};
+	      for (unsigned int i=0;i<fe_rt.dofs_per_cell;++i)
+		{
+		  deallog << "  grad phi(" << i << ") = ";
+		  for (unsigned int c=0; c<dim; ++c)
+                    {
+                      deallog << "[";
+                      for (unsigned int d=0; d<dim; ++d)
+                        deallog << " " << fe.shape_grad_component(i,q_point,c)[d];
+                      deallog << " ]" << std::endl;
+                      if (c != dim-1)
+                        deallog << "                ";
+                    };
+		};
+
+              q_point++;
+	    }
+
+      deallog << std::endl;
+    }
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("rt_2.output");
+  logfile.precision (PRECISION);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  plot_shape_functions<2>();
+//  plot_shape_functions<3>();
+  
+  return 0;
+}
+
+
+
-- 
2.39.5