Tests for Raviart-Thomas elements.

git-svn-id: https://svn.dealii.org/trunk@7427 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/fe/rt_1.cc b/tests/fe/rt_1.cc
new file mode 100644 (file)
index 0000000..1911800
--- /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 (file)
index 0000000..1cb6e33
--- /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;
+}
+
+
+