]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new test
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Mar 2012 23:20:38 +0000 (23:20 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 9 Mar 2012 23:20:38 +0000 (23:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@25231 0785d39b-7218-0410-832d-ea1e28bc413d

tests/integrators/mesh_worker_01.cc
tests/integrators/mesh_worker_02.cc
tests/integrators/mesh_worker_03.cc [new file with mode: 0644]
tests/integrators/mesh_worker_03/cmp/generic [new file with mode: 0644]

index bd25b2be3781bca0aa4cacde104ff594aa25b6c8..653379dc9c1c61144d0028be56705423c2d43432 100644 (file)
@@ -1,8 +1,7 @@
 //----------------------------------------------------------------------
 //    $Id$
-//    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index d87b1c20c8d4d80968ec66d24e5358175cafd2f5..2e859cc7590536d6f44fa66af79255a5319d531d 100644 (file)
@@ -1,8 +1,7 @@
 //----------------------------------------------------------------------
 //    $Id$
-//    Version: $Name$ 
 //
-//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
diff --git a/tests/integrators/mesh_worker_03.cc b/tests/integrators/mesh_worker_03.cc
new file mode 100644 (file)
index 0000000..109298a
--- /dev/null
@@ -0,0 +1,317 @@
+//----------------------------------------------------------------------
+//    $Id: mesh_worker_02.cc 23710 2011-05-17 04:50:10Z bangerth $
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010, 2012 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.
+//
+//----------------------------------------------------------------------
+
+// Test consistency of assemblers MatrixSimple with and without local blocks
+
+#include "../tests.h"
+#include <deal.II/numerics/mesh_worker.h>
+#include <deal.II/numerics/mesh_worker_assembler.h>
+#include <deal.II/numerics/mesh_worker_loop.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/multigrid/mg_tools.h>
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+
+// Local integrators fill every matrix entry with 10*block_row + block_col
+template <int dim>
+class MatrixIntegrator : public Subscriptor
+{
+  public:
+    static void cell(MeshWorker::DoFInfo<dim>& dinfo,
+                    MeshWorker::IntegrationInfo<dim>& info);
+    static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+                    MeshWorker::DoFInfo<dim>& dinfo2,
+                    MeshWorker::IntegrationInfo<dim>& info1,
+                    MeshWorker::IntegrationInfo<dim>& info2);
+    static void block_cell(MeshWorker::DoFInfo<dim>& dinfo,
+                          MeshWorker::IntegrationInfo<dim>& info);
+    static void block_face(MeshWorker::DoFInfo<dim>& dinfo1,
+                          MeshWorker::DoFInfo<dim>& dinfo2,
+                          MeshWorker::IntegrationInfo<dim>& info1,
+                          MeshWorker::IntegrationInfo<dim>& info2);
+};
+
+
+template <int dim>
+void MatrixIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo,
+                                MeshWorker::IntegrationInfo<dim>& info)
+{
+  const FiniteElement<dim>& fe = info.fe_values().get_fe();
+  FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
+  
+  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+    for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+      {
+       local_matrix(i,j) = 10 * fe.system_to_block_index(i).first
+                           + fe.system_to_block_index(j).first;
+      }
+}
+
+
+template <int dim>
+void MatrixIntegrator<dim>::face(MeshWorker::DoFInfo<dim>& dinfo1,
+                                MeshWorker::DoFInfo<dim>& dinfo2,
+                                MeshWorker::IntegrationInfo<dim>& info1,
+                                MeshWorker::IntegrationInfo<dim>& info2)
+{
+  const FiniteElement<dim>& fe1 = info1.fe_values().get_fe();
+  const FiniteElement<dim>& fe2 = info2.fe_values().get_fe();
+  FullMatrix<double>& matrix_v1u2 = dinfo1.matrix(0, true).matrix;
+  FullMatrix<double>& matrix_v2u1 = dinfo2.matrix(0, true).matrix;
+  
+  for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
+    for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
+      {
+       matrix_v1u2(i,j) = 10 * fe1.system_to_block_index(i).first
+                           + fe2.system_to_block_index(j).first;
+       matrix_v2u1(i,j) = 10 * fe2.system_to_block_index(i).first
+                           + fe1.system_to_block_index(j).first;
+      }
+}
+
+
+template <int dim>
+void MatrixIntegrator<dim>::block_cell(
+  MeshWorker::DoFInfo<dim>& dinfo,
+  MeshWorker::IntegrationInfo<dim>&)
+{
+  for (unsigned int m=0;m<dinfo.n_matrices();++m)
+    {
+      MatrixBlock<FullMatrix<double> >& loc = dinfo.matrix(m);
+      
+      for (unsigned int i=0; i<loc.matrix.m(); ++i)
+       for (unsigned int j=0; j<loc.matrix.n(); ++j)
+         {
+           loc.matrix(i,j) = 10 * loc.row + loc.column;
+         }
+    }
+}
+
+
+template <int dim>
+void MatrixIntegrator<dim>::block_face(
+  MeshWorker::DoFInfo<dim>& dinfo1,
+  MeshWorker::DoFInfo<dim>& dinfo2,
+  MeshWorker::IntegrationInfo<dim>&,
+  MeshWorker::IntegrationInfo<dim>&)
+{
+  for (unsigned int m=0;m<dinfo1.n_matrices();++m)
+    {
+      MatrixBlock<FullMatrix<double> >& loc = dinfo1.matrix(m, true);
+      
+      for (unsigned int i=0; i<loc.matrix.m(); ++i)
+       for (unsigned int j=0; j<loc.matrix.n(); ++j)
+         {
+           loc.matrix(i,j) = 10 * loc.row + loc.column;
+         }
+    }
+  for (unsigned int m=0;m<dinfo2.n_matrices();++m)
+    {
+      MatrixBlock<FullMatrix<double> >& loc = dinfo2.matrix(m, true);
+      
+      for (unsigned int i=0; i<loc.matrix.m(); ++i)
+       for (unsigned int j=0; j<loc.matrix.n(); ++j)
+         {
+           loc.matrix(i,j) = 10 * loc.row + loc.column;
+         }
+    }
+}
+
+
+template <int dim>
+void
+assemble(const DoFHandler<dim>& dof_handler, SparseMatrix<double>& matrix)
+{
+  const FiniteElement<dim>& fe = dof_handler.get_fe();
+  MappingQ1<dim> mapping;  
+  
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+  info_box.initialize_update_flags();
+  UpdateFlags update_flags = update_values | update_gradients;
+  info_box.add_update_flags(update_flags, true, true, true, true);
+  info_box.initialize(fe, mapping);
+  
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  
+  MeshWorker::Assembler::MatrixSimple<SparseMatrix<double> > assembler;
+  assembler.initialize(matrix);
+  
+  MeshWorker::integration_loop<dim, dim>
+    (dof_handler.begin_active(),
+     dof_handler.end(),
+     dof_info,
+     info_box,
+     &MatrixIntegrator<dim>::cell,
+     0,//&MatrixIntegrator<dim>::bdry,
+     &MatrixIntegrator<dim>::face,
+     assembler);
+}
+
+
+template <int dim>
+void
+assemble(const MGDoFHandler<dim>& dof_handler,
+        MGLevelObject<SparseMatrix<double> > matrix,
+        MGLevelObject<SparseMatrix<double> > dg_up,
+        MGLevelObject<SparseMatrix<double> > dg_down)
+{
+  const FiniteElement<dim>& fe = dof_handler.get_fe();
+  MappingQ1<dim> mapping;
+  
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+  info_box.initialize_update_flags();
+  UpdateFlags update_flags = update_values | update_gradients;
+  info_box.add_update_flags(update_flags, true, true, true, true);
+  info_box.initialize(fe, mapping);
+  
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  
+  MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
+  assembler.initialize(matrix);
+  
+  MeshWorker::loop<MeshWorker::DoFInfo<dim>, MeshWorker::IntegrationInfoBox<dim> >
+    (dof_handler.begin(),
+     dof_handler.end(),
+     dof_info,
+     info_box,
+     &MatrixIntegrator<dim>::cell,
+     &MatrixIntegrator<dim>::bdry,
+     &MatrixIntegrator<dim>::face,
+     assembler);
+}
+
+
+template <int dim>
+void
+test_simple(MGDoFHandler<dim>& mgdofs)
+{
+  SparsityPattern pattern;
+  SparseMatrix<double> matrix;
+  Vector<double> v;
+
+  const DoFHandler<dim>& dofs = mgdofs;
+  const FiniteElement<dim>& fe = dofs.get_fe();
+  pattern.reinit (dofs.n_dofs(), dofs.n_dofs(),
+                 (GeometryInfo<dim>::faces_per_cell
+                  *GeometryInfo<dim>::max_children_per_face+1)*fe.dofs_per_cell);
+  DoFTools::make_flux_sparsity_pattern (dofs, pattern);
+  pattern.compress();
+  matrix.reinit (pattern);
+  
+  assemble(dofs, matrix);
+  deallog << std::setprecision(3);
+  matrix.print(deallog);
+
+  MGLevelObject<SparsityPattern> mg_sparsity;
+  MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
+  MGLevelObject<SparseMatrix<double> > mg_matrix;
+  MGLevelObject<SparseMatrix<double> > mg_matrix_dg_up;
+  MGLevelObject<SparseMatrix<double> > mg_matrix_dg_down;
+
+  const unsigned int n_levels = mgdofs.get_tria().n_levels();
+  
+  mg_sparsity.resize(0, n_levels-1);
+  mg_sparsity_dg_interface.resize(0, n_levels-1);
+  mg_matrix.resize(0, n_levels-1);
+  mg_matrix_dg_up.resize(0, n_levels-1);
+  mg_matrix_dg_down.resize(0, n_levels-1);
+  
+  for (unsigned int level=mg_sparsity.get_minlevel();
+       level<=mg_sparsity.get_maxlevel();++level)
+    {
+      CompressedSparsityPattern c_sparsity(mgdofs.n_dofs(level));
+      CompressedSparsityPattern ci_sparsity;
+      if (level>0)
+       ci_sparsity.reinit(mgdofs.n_dofs(level-1), mgdofs.n_dofs(level));
+      
+      MGTools::make_flux_sparsity_pattern(mgdofs, c_sparsity, level);
+      if (level>0)
+       MGTools::make_flux_sparsity_pattern_edge(mgdofs, ci_sparsity, level);
+      
+      mg_sparsity[level].copy_from(c_sparsity);
+      mg_matrix[level].reinit(mg_sparsity[level]);
+      if (level>0)
+       {
+         mg_sparsity_dg_interface[level].copy_from(ci_sparsity);
+         mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
+         mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
+       }
+    }
+  
+}
+
+
+template<int dim>
+void
+test(const FiniteElement<dim>& fe)
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_L(tr);
+  tr.begin()->set_refine_flag();
+  tr.execute_coarsening_and_refinement();
+//  tr.begin(2)->set_refine_flag();
+//  tr.execute_coarsening_and_refinement();
+//  tr.refine_global(1);
+  deallog << "Triangulation levels";
+  for (unsigned int l=0;l<tr.n_levels();++l)
+    deallog << ' ' << l << ':' << tr.n_cells(l);
+  deallog << std::endl;
+  
+  unsigned int cn = 0;
+  for (typename Triangulation<dim>::cell_iterator cell = tr.begin();
+       cell != tr.end(); ++cell, ++cn)
+    cell->set_user_index(cn);
+
+  MGDoFHandler<dim> dofs(tr);
+  dofs.distribute_dofs(fe);
+  deallog << "DoFHandler " << dofs.n_dofs() << " levels";
+  for (unsigned int l=0;l<tr.n_levels();++l)
+    deallog << ' ' << l << ':' << dofs.n_dofs(l);
+  deallog << std::endl;
+  
+  test_simple(dofs);
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("mesh_worker_02/output");
+  logfile << std::setprecision (2);
+  logfile << std::fixed;  
+  deallog << std::setprecision (2);
+  deallog << std::fixed;  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  std::vector<std_cxx1x::shared_ptr<FiniteElement<2> > > fe2;
+  fe2.push_back(std_cxx1x::shared_ptr<FiniteElement<2> >(new  FE_DGP<2>(1)));
+//  fe2.push_back(std_cxx1x::shared_ptr<FiniteElement<2> >(new  FE_Q<2>(1)));
+  
+  for (unsigned int i=0;i<fe2.size();++i)
+    test(*fe2[i]);
+}
diff --git a/tests/integrators/mesh_worker_03/cmp/generic b/tests/integrators/mesh_worker_03/cmp/generic
new file mode 100644 (file)
index 0000000..85ba2f0
--- /dev/null
@@ -0,0 +1,201 @@
+
+DEAL::Triangulation levels 0:3 1:4
+DEAL::DoFHandler 18 levels 0:9 1:12
+DEAL::(0,0) 15.000
+DEAL::(0,1) 0.000
+DEAL::(0,2) 0
+DEAL::(0,9) -1.500
+DEAL::(0,10) -0.866
+DEAL::(0,11) 0
+DEAL::(0,15) -1.500
+DEAL::(0,16) -0.866
+DEAL::(0,17) 0
+DEAL::(1,1) 23.000
+DEAL::(1,0) 0.000
+DEAL::(1,2) 0
+DEAL::(1,9) 1.732
+DEAL::(1,10) 0
+DEAL::(1,11) 0
+DEAL::(1,15) 1.732
+DEAL::(1,16) 0
+DEAL::(1,17) 0
+DEAL::(2,2) 19.000
+DEAL::(2,0) 0
+DEAL::(2,1) 0
+DEAL::(2,9) 1.299
+DEAL::(2,10) 0.750
+DEAL::(2,11) -0.750
+DEAL::(2,15) -1.299
+DEAL::(2,16) -0.750
+DEAL::(2,17) -0.750
+DEAL::(3,3) 15.000
+DEAL::(3,4) 0
+DEAL::(3,5) 0.000
+DEAL::(3,12) -1.500
+DEAL::(3,13) 0
+DEAL::(3,14) -0.866
+DEAL::(3,15) -1.500
+DEAL::(3,16) 0
+DEAL::(3,17) -0.866
+DEAL::(4,4) 19.000
+DEAL::(4,3) 0
+DEAL::(4,5) 0
+DEAL::(4,12) 1.299
+DEAL::(4,13) -0.750
+DEAL::(4,14) 0.750
+DEAL::(4,15) -1.299
+DEAL::(4,16) -0.750
+DEAL::(4,17) -0.750
+DEAL::(5,5) 23.000
+DEAL::(5,3) 0.000
+DEAL::(5,4) 0
+DEAL::(5,12) 1.732
+DEAL::(5,13) 0
+DEAL::(5,14) 0
+DEAL::(5,15) 1.732
+DEAL::(5,16) 0
+DEAL::(5,17) 0
+DEAL::(6,6) 12.000
+DEAL::(6,7) -1.732
+DEAL::(6,8) -1.732
+DEAL::(6,9) -2.000
+DEAL::(6,10) 1.732
+DEAL::(6,11) 0
+DEAL::(6,12) -2.000
+DEAL::(6,13) 0
+DEAL::(6,14) 1.732
+DEAL::(7,7) 18.000
+DEAL::(7,6) -1.732
+DEAL::(7,8) 0
+DEAL::(7,9) -1.732
+DEAL::(7,10) 0
+DEAL::(7,11) 0
+DEAL::(7,12) 0
+DEAL::(7,13) -2.000
+DEAL::(7,14) 0
+DEAL::(8,8) 18.000
+DEAL::(8,6) -1.732
+DEAL::(8,7) 0
+DEAL::(8,9) 0
+DEAL::(8,10) 0
+DEAL::(8,11) -2.000
+DEAL::(8,12) -1.732
+DEAL::(8,13) 0
+DEAL::(8,14) 0
+DEAL::(9,9) 9.500
+DEAL::(9,0) -1.500
+DEAL::(9,1) 1.732
+DEAL::(9,2) 1.299
+DEAL::(9,6) -2.000
+DEAL::(9,7) -1.732
+DEAL::(9,8) 0
+DEAL::(9,10) -0.866
+DEAL::(9,11) -1.732
+DEAL::(9,15) -2.000
+DEAL::(9,16) 0
+DEAL::(9,17) 1.732
+DEAL::(10,10) 16.500
+DEAL::(10,0) -0.866
+DEAL::(10,1) 0
+DEAL::(10,2) 0.750
+DEAL::(10,6) 1.732
+DEAL::(10,7) 0
+DEAL::(10,8) 0
+DEAL::(10,9) -0.866
+DEAL::(10,11) 0
+DEAL::(10,15) 0
+DEAL::(10,16) -2.000
+DEAL::(10,17) 0
+DEAL::(11,11) 15.500
+DEAL::(11,0) 0
+DEAL::(11,1) 0
+DEAL::(11,2) -0.750
+DEAL::(11,6) 0
+DEAL::(11,7) 0
+DEAL::(11,8) -2.000
+DEAL::(11,9) -1.732
+DEAL::(11,10) 0
+DEAL::(11,15) -1.732
+DEAL::(11,16) 0
+DEAL::(11,17) 0
+DEAL::(12,12) 9.500
+DEAL::(12,3) -1.500
+DEAL::(12,4) 1.299
+DEAL::(12,5) 1.732
+DEAL::(12,6) -2.000
+DEAL::(12,7) 0
+DEAL::(12,8) -1.732
+DEAL::(12,13) -1.732
+DEAL::(12,14) -0.866
+DEAL::(12,15) -2.000
+DEAL::(12,16) 1.732
+DEAL::(12,17) 0
+DEAL::(13,13) 15.500
+DEAL::(13,3) 0
+DEAL::(13,4) -0.750
+DEAL::(13,5) 0
+DEAL::(13,6) 0
+DEAL::(13,7) -2.000
+DEAL::(13,8) 0
+DEAL::(13,12) -1.732
+DEAL::(13,14) 0
+DEAL::(13,15) -1.732
+DEAL::(13,16) 0
+DEAL::(13,17) 0
+DEAL::(14,14) 16.500
+DEAL::(14,3) -0.866
+DEAL::(14,4) 0.750
+DEAL::(14,5) 0
+DEAL::(14,6) 1.732
+DEAL::(14,7) 0
+DEAL::(14,8) 0
+DEAL::(14,12) -0.866
+DEAL::(14,13) 0
+DEAL::(14,15) 0
+DEAL::(14,16) 0
+DEAL::(14,17) -2.000
+DEAL::(15,15) 7.000
+DEAL::(15,0) -1.500
+DEAL::(15,1) 1.732
+DEAL::(15,2) -1.299
+DEAL::(15,3) -1.500
+DEAL::(15,4) -1.299
+DEAL::(15,5) 1.732
+DEAL::(15,9) -2.000
+DEAL::(15,10) 0
+DEAL::(15,11) -1.732
+DEAL::(15,12) -2.000
+DEAL::(15,13) -1.732
+DEAL::(15,14) 0
+DEAL::(15,16) -0.866
+DEAL::(15,17) -0.866
+DEAL::(16,16) 14.000
+DEAL::(16,0) -0.866
+DEAL::(16,1) 0
+DEAL::(16,2) -0.750
+DEAL::(16,3) 0
+DEAL::(16,4) -0.750
+DEAL::(16,5) 0
+DEAL::(16,9) 0
+DEAL::(16,10) -2.000
+DEAL::(16,11) 0
+DEAL::(16,12) 1.732
+DEAL::(16,13) 0
+DEAL::(16,14) 0
+DEAL::(16,15) -0.866
+DEAL::(16,17) 0
+DEAL::(17,17) 14.000
+DEAL::(17,0) 0
+DEAL::(17,1) 0
+DEAL::(17,2) -0.750
+DEAL::(17,3) -0.866
+DEAL::(17,4) -0.750
+DEAL::(17,5) 0
+DEAL::(17,9) 1.732
+DEAL::(17,10) 0
+DEAL::(17,11) 0
+DEAL::(17,12) 0
+DEAL::(17,13) 0
+DEAL::(17,14) -2.000
+DEAL::(17,15) -0.866
+DEAL::(17,16) 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.