]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add RT tests.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 9 Jun 2003 15:59:07 +0000 (15:59 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 9 Jun 2003 15:59:07 +0000 (15:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@7771 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fe/Makefile
tests/fe/rt_5.cc [new file with mode: 0644]
tests/fe/rt_6.cc [new file with mode: 0644]

index 16b2535e5b77cdfd5120afa31b55196120fd2d39..2ba9d037c328f754c5490b14d376d2c22b297fbd 100644 (file)
@@ -45,6 +45,8 @@ rt_1.exe                : rt_1.g.$(OBJEXT)                $(libraries)
 rt_2.exe                : rt_2.g.$(OBJEXT)                $(libraries)
 rt_3.exe                : rt_3.g.$(OBJEXT)                $(libraries)
 rt_4.exe                : rt_4.g.$(OBJEXT)                $(libraries)
+rt_5.exe                : rt_5.g.$(OBJEXT)                $(libraries)
+rt_6.exe                : rt_6.g.$(OBJEXT)                $(libraries)
 system_1.exe            : system_1.g.$(OBJEXT)            $(libraries)
 shapes.exe              : shapes.g.$(OBJEXT)              $(libraries)
 transfer.exe           : transfer.g.$(OBJEXT)            $(libraries)
diff --git a/tests/fe/rt_5.cc b/tests/fe/rt_5.cc
new file mode 100644 (file)
index 0000000..d4c9ff7
--- /dev/null
@@ -0,0 +1,69 @@
+//----------------------------  rt_5.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003 by the deal.II authors
+//
+//    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_5.cc  ---------------------------
+
+
+// Just output the restriction matrices of the RT element
+
+#include <base/logstream.h>
+#include <fe/fe_raviart_thomas.h>
+
+#include <fstream>
+#include <string>
+
+#define PRECISION 5
+
+
+
+template<int dim>
+void
+test(const unsigned int degree)
+{
+  deallog << "FE_RaviartThomas<" << dim << "> (" << degree << ")"
+         << std::endl;
+  
+  FE_RaviartThomas<dim> fe_rt(degree);
+
+  for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+    {
+      const FullMatrix<double> & m = fe_rt.restrict(c);
+
+      for (unsigned int i=0; i<m.m(); ++i)
+       {
+         for (unsigned int j=0; j<m.n(); ++j)
+           deallog << m(i,j) << ' ';
+         deallog << std::endl;
+       }
+      
+      deallog << std::endl;
+    }
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("rt_5.output");
+  logfile.precision (PRECISION);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  for (unsigned int degree=0; degree<1; ++degree)
+    test<2>(degree);
+//  test<3>(degree);
+  
+  return 0;
+}
+
+
+
diff --git a/tests/fe/rt_6.cc b/tests/fe/rt_6.cc
new file mode 100644 (file)
index 0000000..0abac9a
--- /dev/null
@@ -0,0 +1,174 @@
+//----------------------------  rt_6.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003 by the deal.II authors
+//
+//    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_6.cc  ---------------------------
+
+// adaptation of up_and_down for RT elements
+
+#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_nedelec.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_values.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 4
+
+
+template <int dim>
+Point<dim> transform (const Point<dim> p)
+{
+  switch (dim) 
+    {
+      case 1:
+            return p;
+      case 2:
+            return Point<dim>(p(0)*(1+p(1)), p(1)*(1+p(0)));
+      case 3:
+            return Point<dim>(p(0)*(1+p(1))*(1+p(2)),
+                              p(1)*(1+p(0))*(1+p(2)),
+                              p(2)*(1+p(0))*(1+p(1)));
+      default:
+            Assert (false, ExcNotImplemented());
+            return Point<dim>();
+    };
+}
+
+
+template <int dim>
+void check_element (const Triangulation<dim> &tr,
+                    const FiniteElement<dim> &fe)
+{
+  DoFHandler<dim> dof_handler(tr);
+  dof_handler.distribute_dofs (fe);
+
+                                   // create a mostly arbitrary
+                                   // function plus a trend on this
+                                   // grid
+  Vector<double> tmp(dof_handler.n_dofs());
+  for (unsigned int i=0; i<tmp.size(); ++i)
+    tmp(i) = (i + 13*i%17);
+
+                                   // restrict this function to the
+                                   // next coarser level and
+                                   // distribute it again to the
+                                   // higher level
+  Vector<double> x(tmp.size());
+  Vector<double> v(fe.dofs_per_cell);
+  for (typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+       cell!=dof_handler.end(); ++cell)
+    if (cell->has_children() &&
+        cell->child(0)->active())
+      {
+                                         // first make sure that what
+                                         // we do is reasonable. for
+                                         // this, _all_ children have
+                                         // to be active, not only
+                                         // some of them
+        for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+          Assert (cell->child(c)->active(), ExcInternalError());
+
+                                         // then restrict and prolong
+        cell->get_interpolated_dof_values (tmp, v);
+        cell->set_dof_values_by_interpolation (v, x);
+      };
+
+                                   // now x is a function on the fine
+                                   // grid that is representable on
+                                   // the coarse grid. so another
+                                   // cycle should not alter it any
+                                   // more:
+  Vector<double> x2(x.size());
+  for (typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+       cell!=dof_handler.end(); ++cell)
+    if (cell->has_children() &&
+        cell->child(0)->active())
+      {
+        cell->get_interpolated_dof_values (x, v);
+        cell->set_dof_values_by_interpolation (v, x2);
+      };
+
+                                   // then check that this is so:
+  x2 -= x;
+  const double relative_residual = (x2.l2_norm() / x.l2_norm());
+
+  const double threshold = 1e-6;
+  deallog << ", dofs_per_cell=" << fe.dofs_per_cell
+          << "; relative residual: "
+          << (relative_residual<threshold ? "ok" : "botched up!")
+         << " (relative residual=" << relative_residual << ")"
+          << std::endl;
+  
+  Assert (relative_residual < threshold*x.l2_norm(), ExcInternalError());
+}
+
+
+template <int dim>
+void test () 
+{
+                                   // make a coarse triangulation as a
+                                   // hypercube. if in more than 1d,
+                                   // distort it so that it is no more
+                                   // an affine image of the
+                                   // hypercube, to make things more
+                                   // difficult. then refine it twice
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr, 0., 1.);
+  Point<dim> (*p) (Point<dim>) = &transform<dim>;
+  GridTools::transform (p, tr);
+  tr.refine_global (2);
+
+                                   // now for a list of finite
+                                   // elements, for which we want to
+                                   // test. we happily waste tons of
+                                   // memory here, but who cares...
+  const FiniteElement<dim> *fe_list[]
+    = 
+    {
+          new FE_RaviartThomas<dim>(0),
+          new FE_RaviartThomas<dim>(1)
+    };
+  
+  for (unsigned int i=0; i<sizeof(fe_list)/sizeof(fe_list[0]); ++i)
+    if (fe_list[i] != 0)
+      {
+        deallog << dim << "d, uniform grid, fe #" << i;
+        check_element (tr, *fe_list[i]);
+      }
+}
+
+
+
+int
+main()
+{
+  std::ofstream logfile ("rt_6.output");
+  logfile.precision (PRECISION);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  test<2>();
+  
+  return 0;
+}
+
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.