]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
still debugging
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 May 2005 05:36:35 +0000 (05:36 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 May 2005 05:36:35 +0000 (05:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@10737 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_poly_tensor.cc
deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc
deal.II/deal.II/source/fe/fe_tools.cc
tests/fe/Makefile
tests/fe/rtn_1.cc [new file with mode: 0644]

index b6fb9f3a8f2a7f84e176b0afcf4ea077220eeb97..1d40c915f155ec0e676ce7961b26aee7ea9f41a7 100644 (file)
@@ -16,6 +16,7 @@
 #include <fe/fe_poly_tensor.h>
 #include <fe/fe_values.h>
 
+#include <base/logstream.h>
 
 template <class POLY, int dim>
 FE_PolyTensor<POLY,dim>::FE_PolyTensor (
@@ -64,7 +65,7 @@ FE_PolyTensor<POLY,dim>::shape_value_component (
     return cached_values[i][component];
   else
     for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
-      s += inverse_node_matrix(i,j) * cached_values[j][component];
+      s += inverse_node_matrix(j,i) * cached_values[j][component];
   return s;
 }
 
@@ -103,7 +104,7 @@ FE_PolyTensor<POLY,dim>::shape_grad_component (
     return cached_grads[i][component];
   else
     for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
-      s += inverse_node_matrix(i,j) * cached_grads[j][component];
+      s += inverse_node_matrix(j,i) * cached_grads[j][component];
   
   return s;
 }
@@ -174,7 +175,9 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
 
   const UpdateFlags flags(data->update_flags);
   const unsigned int n_q_points = quadrature.n_quadrature_points;
-
+  double d;
+  d = floor (1.3);
+  
                                   // some scratch arrays
   std::vector<Tensor<1,dim> > values(0);
   std::vector<Tensor<2,dim> > grads(0);
@@ -228,7 +231,7 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
          else
            for (unsigned int i=0; i<this->dofs_per_cell; ++i)
              for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-               data->shape_values[i][k] = inverse_node_matrix(i,j) * values[j];
+               data->shape_values[i][k] = inverse_node_matrix(j,i) * values[j];
        
        if (flags & update_gradients)
          if (inverse_node_matrix.n_cols() == 0)
@@ -237,7 +240,7 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
          else
            for (unsigned int i=0; i<this->dofs_per_cell; ++i)
              for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-               data->shape_grads[i][k] = inverse_node_matrix(i,j) * grads[j];
+               data->shape_grads[i][k] = inverse_node_matrix(j,i) * grads[j];
       }
   return data;
 }
@@ -279,7 +282,7 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
-      const unsigned int first = data.shape_function_to_row_table[i];
+      const unsigned int first = i*dim;//data.shape_function_to_row_table[i];
       
       if (flags & update_values)
        for (unsigned int k=0; k<n_quad; ++k)
@@ -302,6 +305,14 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
   if (flags & update_second_derivatives)
     this->compute_2nd (mapping, cell, dsd.cell(),
                       mapping_data, fe_data, data);
+  deallog << "Shape data " << get_name() << std::endl;
+  
+  for (unsigned int i=0;i<data.shape_values.n_rows();++i)
+    {
+      for (unsigned int j=0;j<n_quad;++j)
+       deallog << ' ' << data.shape_values(i,j);
+      deallog << std::endl; 
+    }
 }
 
 
index 25ff458fa3ba284f1349dfb3e0718fd10a176af3..8f57d435a254a72d58f7c00c9823f2a80c6cb3f5 100644 (file)
@@ -45,15 +45,16 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   
   this->mapping_type = this->independent_on_cartesian;
-  
+                                  // These must be done first, since
+                                  // they change the evaluation of
+                                  // basis functions
+  initialize_unit_support_points (deg);
+  initialize_node_matrix();
+
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
     this->prolongation[i].reinit (this->dofs_per_cell,
                                  this->dofs_per_cell);
   FETools::compute_embedding_matrices (*this, &this->prolongation[0]);
-                                  // finally fill in support points
-                                  // on cell and face
-  initialize_unit_support_points (deg);
-  initialize_node_matrix();
 }
 
 
index faed071f6d9186364ecd17b55f285f589993c2d4..76cb76474603c046c24da5f9affc9b9300942cd9 100644 (file)
@@ -602,6 +602,88 @@ FETools::compute_embedding_matrices(const FiniteElement<dim>& fe,
 }
 
 
+template<int dim, typename number>
+void
+FETools::compute_projection_matrices(const FiniteElement<dim>& fe,
+                                    FullMatrix<number>* matrices)
+{
+  const unsigned int nc = GeometryInfo<dim>::children_per_cell;
+  const unsigned int n  = fe.dofs_per_cell;
+  const unsigned int nd = fe.n_components();
+  const unsigned int degree = fe.degree;
+  
+  for (unsigned int i=0;i<nc;++i)
+    {
+      Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n(),n));
+      Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m(),n));
+    }
+  
+                                  /*
+                                   * Set up two meshes, one with a
+                                   * single reference cell and the
+                                   * other refined, together with
+                                   * DoFHandler and finite elements.
+                                   */
+  Triangulation<dim> tr_coarse;
+  Triangulation<dim> tr_fine;
+  GridGenerator::hyper_cube (tr_coarse, 0, 1);
+  GridGenerator::hyper_cube (tr_fine, 0, 1);
+  tr_fine.refine_global(1);
+  DoFHandler<dim> dof_coarse(tr_coarse);
+  dof_coarse.distribute_dofs(fe);
+  DoFHandler<dim> dof_fine(tr_fine);
+  dof_fine.distribute_dofs(fe);
+
+  MappingCartesian<dim> mapping;
+  QGauss<dim> q_fine(degree+1);
+  const unsigned int nq = q_fine.n_quadrature_points;
+  
+  FEValues<dim> coarse (mapping, fe, q_fine,
+                       update_q_points | update_JxW_values | update_values);
+  FEValues<dim> fine (mapping, fe, q_fine,
+                     update_q_points | update_JxW_values | update_values);
+  
+  typename DoFHandler<dim>::active_cell_iterator coarse_cell
+    = dof_coarse.begin_active();
+  typename DoFHandler<dim>::active_cell_iterator fine_cell;
+
+                                  // Compute the coarse level mass
+                                  // matrix
+  coarse.reinit(dof_coarse);
+  FullMatrix<number> A(n, n);
+  for (unsigned int k=0;k<nq;++k)
+    for (unsigned int i=0;i<n;++j)
+      for (unsigned int j=0;j<n;++j)
+       A(i,j) = coarse.JxW(k)
+                * coarse.shape_value(i,k)
+                * coarse.shape_value(j,k);
+  
+  Householder<double> H(A);
+  
+  Vector<number> v_coarse(n);
+  Vector<number> v_fine(n);
+  
+  unsigned int cell_number = 0;
+  for (fine_cell = dof_fine.begin_active();
+       fine_cell != dof_fine.end();
+       ++fine_cell, ++cell_number)
+    {
+                                      // Compute right hand side,
+                                      // which is a fine level basis
+                                      // function tested with the
+                                      // coarse level functions.
+      Assert(false, ExcNotImplemented());
+      
+                                      // Remove small entries from
+                                      // the matrix
+      for (unsigned int i=0; i<matrix.m(); ++i)
+       for (unsigned int j=0; j<matrix.n(); ++j)
+         if (std::fabs(matrix(i,j)) < 1e-12)
+           matrix(i,j) = 0.;
+    }
+}
+
+
 template <int dim,
           class InVector, class OutVector>
 void
index 8a70014ed71ffa4dbcbbd52670d251ae486812c0..150a22d78c8bf026284d0c45fedf750c45419767 100644 (file)
@@ -1,7 +1,6 @@
-# Generated automatically from Makefile.in by configure.
 ############################################################
 # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp
-# Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -49,6 +48,7 @@ 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)
+rtn_1.exe               : rtn_1.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)
diff --git a/tests/fe/rtn_1.cc b/tests/fe/rtn_1.cc
new file mode 100644 (file)
index 0000000..170ef72
--- /dev/null
@@ -0,0 +1,92 @@
+//----------------------------------------------------------------------
+//    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
+// Show the shape functions in support points
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/mapping_cartesian.h>
+#include <fe/fe_values.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#include <cmath>
+extern "C" 
+{
+#include <math.h>
+};
+
+template <int dim>
+void
+check_support_points (const FiniteElement<dim>& fe)
+{
+  deallog << fe.get_name() << std::endl;
+  
+  const std::vector< Point<dim > > & points = fe.get_unit_support_points();
+  std::vector<double> weights(points.size());
+
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+  
+  MappingCartesian<dim> mapping;
+  
+  Quadrature<dim> support_quadrature(points, weights);
+  UpdateFlags flags = update_values;
+  FEValues<dim> vals(mapping, fe, support_quadrature, flags);
+  vals.reinit(dof.begin_active());
+  
+  for (unsigned int k=0;k<points.size();++k)
+    {
+      const Point<dim>& p = points[k];
+      deallog << p;
+      for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+       {
+         deallog << '\t';
+         for (unsigned int d=0;d<dim;++d)
+           {
+             const double sf = fe.shape_value_component(i,p,d);
+             const double sv = vals.shape_value_component(i,k,d);
+             deallog << ' ' << (int) round(60*sf)
+                     << '/' << (int) round(60*sv);
+             
+//           const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
+//           if (diff > 1.e-12)
+//             deallog << "Error" << vals.shape_value_component(i,k,d) << std::endl;
+           }
+       }
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("rtn_1.output");
+  deallog.attach(logfile);
+  deallog.depth_console(10);
+  deallog.threshold_double(1.e-10);
+
+  FE_RaviartThomasNodal<2> e20(0);
+  check_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.