]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add tests for new nodal RT element
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 30 May 2005 09:08:21 +0000 (09:08 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 30 May 2005 09:08:21 +0000 (09:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@10775 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fe/Makefile
tests/fe/fe_data_test.cc
tests/fe/rtdiff.cc [new file with mode: 0644]
tests/fe/rtn_1.cc

index 150a22d78c8bf026284d0c45fedf750c45419767..02b1a13a1e25accb5f2138f69b090e6d0a870b8e 100644 (file)
@@ -49,6 +49,7 @@ rt_4.exe                : rt_4.g.$(OBJEXT)                $(libraries)
 rt_5.exe                : rt_5.g.$(OBJEXT)                $(libraries)
 rt_6.exe                : rt_6.g.$(OBJEXT)                $(libraries)
 rtn_1.exe               : rtn_1.g.$(OBJEXT)               $(libraries)
+rtdiff.exe              : rtdiff.g.$(OBJEXT)              $(libraries)
 dgp_monomial_1.exe     : dgp_monomial_1.g.$(OBJEXT)      $(libraries)
 dgp_monomial_2.exe     : dgp_monomial_2.g.$(OBJEXT)      $(libraries)
 system_1.exe            : system_1.g.$(OBJEXT)            $(libraries)
@@ -61,7 +62,7 @@ tests = fe_data_test traits fe_tools fe_tools_test mapping \
         transfer internals \
        dgq_1 \
        q_1 q_2 q_3 q_4 \
-       rt_1 rt_2 rt_3 rt_4 rt_5 rt_6 \
+       rt_1 rt_2 rt_3 rt_4 rt_5 rt_6 rtdiff rtn_1 \
        dgp_monomial_1 dgp_monomial_2 \
        system_1 \
         non_primitive_1 non_primitive_2 \
index add70c3e632a99ba024b9a1d3ef8c6b30d56b880..072a957536a1e3265623a969ee1d00c214369d1d 100644 (file)
@@ -59,12 +59,29 @@ void test_fe_datas()
                                   // Check Raviart-Thomas in 2d only
   if (dim==2)
     {
-      fe_datas.push_back(new FE_RaviartThomas<dim>(0));
+      FE_RaviartThomas<dim>* rt0 = new FE_RaviartThomas<dim>(0);
+      FE_RaviartThomas<dim>* rt1 = new FE_RaviartThomas<dim>(1);
+      fe_datas.push_back(rt0);
       deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-      fe_datas.push_back(new FE_RaviartThomas<dim>(1));
+      fe_datas.push_back(rt1);
       deallog << (*fe_datas.rbegin())->get_name() << std::endl;
       fe_datas.push_back(new FE_RaviartThomas<dim>(2));
       deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+      fe_datas.push_back(new FESystem<dim>(*rt1, 1,
+                                          FE_DGQ<dim> (1), 1));
+      deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+    }
+  if (dim>1)
+    {
+      FE_RaviartThomasNodal<dim>* rt0 = new FE_RaviartThomasNodal<dim>(0);
+      FE_RaviartThomasNodal<dim>* rt1 = new FE_RaviartThomasNodal<dim>(1);
+      fe_datas.push_back(rt0);
+      deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+      fe_datas.push_back(rt1);
+      deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+      fe_datas.push_back(new FESystem<dim>(*rt1, 1,
+                                          FE_DGQ<dim> (1), 1));
+      deallog << (*fe_datas.rbegin())->get_name() << std::endl;
     }
   
 
diff --git a/tests/fe/rtdiff.cc b/tests/fe/rtdiff.cc
new file mode 100644 (file)
index 0000000..d8bd721
--- /dev/null
@@ -0,0 +1,150 @@
+//----------------------------------------------------------------------
+//    rt_1.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
+//    Version: 
+//
+//    Copyright (C) 2005 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.
+//
+//----------------------------------------------------------------------
+
+// FE_RaviartThomasNodal and FE_RaviartThomas
+
+// compare the shape funcions and shape values after converting to the
+// same basis.
+
+// Summary: the different Raviart-Thomas implementations use the same
+// polynomial spaces, but different basis functions. Here, we convert
+// between the bases and test if the resulting functions are the same
+// point-wise.
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/mapping_cartesian.h>
+#include <fe/fe_values.h>
+
+#include <fstream>
+
+// This function copied from FERaviartThomasNodal. nodes is the
+// element having the support points and the value of other in these
+// points is computed.
+template <int dim>
+void
+initialize_node_matrix (const FiniteElement<dim>& other,
+                       const FiniteElement<dim>& nodes,
+                       FullMatrix<double>& N)
+{
+  const unsigned int n_dofs = other.dofs_per_cell;
+  Assert (n_dofs == nodes.dofs_per_cell,
+         ExcDimensionMismatch(n_dofs, nodes.dofs_per_cell));
+  
+  N.reinit(n_dofs, n_dofs);
+
+  const std::vector<Point<dim> >& unit_support_points = nodes.get_unit_support_points();
+  
+                                  // The curent node functional index
+  unsigned int current = 0;
+
+                                  // For each face and all quadrature
+                                  // points on it, the node value is
+                                  // the normal component of the
+                                  // shape function, possibly
+                                  // pointing in negative direction.
+  for (unsigned int face = 0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int k=0;k<other.dofs_per_face;++k)
+      {
+       for (unsigned int i=0;i<n_dofs;++i)
+         N(current,i) = other.shape_value_component(
+           i, unit_support_points[current],
+           GeometryInfo< dim >::unit_normal_direction[face]);
+       ++current;
+      }
+                                  // Interior degrees of freedom in each direction
+  const unsigned int n_cell = (n_dofs - current) / dim;
+  
+  for (unsigned int d=0;d<dim;++d)
+    for (unsigned int k=0;k<n_cell;++k)
+      {
+       for (unsigned int i=0;i<n_dofs;++i)
+         N(current,i) = other.shape_value_component(i, unit_support_points[current], d);
+       ++current;
+      }
+  Assert (current == n_dofs, ExcInternalError());
+}
+
+
+template <int dim>
+void
+compare_shapes (const FiniteElement<dim>& other,
+               const FiniteElement<dim>& nodes,
+               FullMatrix<double>& M)
+{
+  QGauss<dim> quadrature(other.degree+1);
+  Table<3,double> other_values(quadrature.n_quadrature_points, other.dofs_per_cell, dim);
+  Table<3,double> nodes_values(quadrature.n_quadrature_points, other.dofs_per_cell, dim);
+  Table<3,Tensor<1,dim> > other_grads(quadrature.n_quadrature_points, other.dofs_per_cell, dim);
+  Table<3,Tensor<1,dim> > nodes_grads(quadrature.n_quadrature_points, other.dofs_per_cell, dim);
+  for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+    for (unsigned int i=0;i<other.dofs_per_cell;++i)
+      for (unsigned int d=0;d<dim;++d)
+       {
+         other_values[k][i][d] = other.shape_value_component(i,quadrature.point(k),d);
+         nodes_values[k][i][d] = nodes.shape_value_component(i,quadrature.point(k),d);
+         other_grads[k][i][d] = other.shape_grad_component(i,quadrature.point(k),d);
+         nodes_grads[k][i][d] = nodes.shape_grad_component(i,quadrature.point(k),d);
+       }
+
+  for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+    {
+      for (unsigned int i=0;i<other.dofs_per_cell;++i)
+       for (unsigned int d=0;d<dim;++d)
+         {
+           double value = other_values[k][i][d];
+           Tensor<1,dim> grad = other_grads[k][i][d];
+           for (unsigned int j=0;j<other.dofs_per_cell;++j)
+             {
+               value -= M(j,i) * nodes_values[k][j][d];              
+               grad -= M(j,i) * nodes_grads[k][j][d];
+             }
+           deallog << '.';
+           if (std::fabs(value) > 1.e-12)
+             deallog << "Error value\t" << k << '\t' << i << '\t' << d << '\t' << value
+                     << std::endl;
+           if (grad.norm() > 1.e-12)
+             deallog << "Error grad\t" << k << '\t' << i << '\t' << d << '\t' << grad
+                     << '\t' << other_grads[k][i][d]
+                     << std::endl;
+         }
+      deallog << std::endl;
+    }
+}
+
+
+template<int dim>
+void
+test (unsigned int degree)
+{
+  FE_RaviartThomas<dim> rt1(degree);
+  FE_RaviartThomasNodal<dim> rtn1(degree);
+  FullMatrix<double> N;
+  initialize_node_matrix(rt1, rtn1, N);
+  compare_shapes(rt1, rtn1, N);
+}
+
+int main()
+{
+  std::ofstream logfile ("rtdiff.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<2>(0);
+  test<2>(1);
+  test<2>(2);  
+}
index 75f7ce463a8c213e093c6dec075f1efb275057f7..2fbb1fba7a94cff390cd5fa9aa988b2912feb375 100644 (file)
@@ -138,11 +138,13 @@ main()
 {
   std::ofstream logfile ("rtn_1.output");
   deallog.attach(logfile);
-  deallog.depth_console(10);
+  deallog.depth_console(0);
   deallog.threshold_double(1.e-10);
 
   FE_RaviartThomasNodal<2> e20(0);
   check_support_points(e20);
+  FE_RaviartThomasNodal<3> e30(0);
+  check_support_points(e30);
   check_face_support_points(e20);
   FE_RaviartThomasNodal<2> e21(1);
   check_support_points(e21);

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.