]> https://gitweb.dealii.org/ - dealii.git/commitdiff
local blocks are in DoFInfo now only
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 14 Mar 2013 22:24:11 +0000 (22:24 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 14 Mar 2013 22:24:11 +0000 (22:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@28905 0785d39b-7218-0410-832d-ea1e28bc413d

tests/integrators/assembler_simple_matrix_01.cc
tests/integrators/assembler_simple_matrix_01/cmp/generic
tests/integrators/assembler_simple_mgmatrix_01.cc
tests/integrators/assembler_simple_mgmatrix_01/cmp/generic
tests/integrators/elasticity_01.cc [new file with mode: 0644]

index 50e180a96176b5322571cb9c210d0bf4027e7a32..43de184d5c74f5973459e99b04cae630943baa07 100644 (file)
@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2012 by the deal.II authors
+//    Copyright (C) 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 
 /**
  * @file Test initialization of Assembler::MatrixSimple and
- * LocalResults
+ * DoFInfo
  */
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/block_indices.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/meshworker/local_results.h>
+#include <deal.II/meshworker/dof_info.h>
 #include <deal.II/meshworker/simple.h>
 
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_system.h>
+
 using namespace dealii;
 
-template <class MATRIX>
-void test(MeshWorker::Assembler::MatrixSimple<MATRIX>& ass)
-{
-  MeshWorker::LocalResults<double> info1;
-  ass.initialize_info(info1, false);
+template <class DOFINFO, class MATRIX>
+void test(DOFINFO& info, MeshWorker::Assembler::MatrixSimple<MATRIX>& ass)
+{  
+  ass.initialize_info(info, false);
   deallog << "No faces" << std::endl;
-  info1.print_debug(deallog);
+  info.print_debug(deallog);
   
-  MeshWorker::LocalResults<double> info2;
-  ass.initialize_info(info2, true);
+  ass.initialize_info(info, true);
   deallog << "With faces" << std::endl;
-  info2.print_debug(deallog);
+  info.print_debug(deallog);
 }
 
 int main()
@@ -45,23 +50,32 @@ int main()
   deallog.attach(logfile);
   deallog.depth_console (0);
   
+  Triangulation<2,2> tr;
+  GridGenerator::hyper_cube(tr);
+  FE_DGP<2,2> fe1(1);
+  FE_DGP<2,2> fe2(2);
+  FE_DGP<2,2> fe3(3);
+  FE_DGP<2,2> fe5(5);
+  FESystem<2,2> fes1(fe3,1,fe5,1,fe1,1);
+  FESystem<2,2> fes2(fe3,1,fe5,1,fe1,1, fe2,1);
+
+  DoFHandler<2,2> dof1(tr);
+  dof1.distribute_dofs(fes1);
+  DoFHandler<2,2> dof2(tr);
+  dof2.distribute_dofs(fes2);    
+  dof1.initialize_local_block_info();
+  dof2.initialize_local_block_info();
+  MeshWorker::DoFInfo<2,2,double> info1(dof1);
+  MeshWorker::DoFInfo<2,2,double> info1b(dof1.block_info());
+  MeshWorker::DoFInfo<2,2,double> info2b(dof2.block_info());
+
   MeshWorker::Assembler::MatrixSimple<FullMatrix<double> > ass1;
   deallog.push("Single block");
-  test(ass1);
-  BlockIndices ind;
-  ind.push_back(3);
-  ind.push_back(5);
-  ind.push_back(1);
-  ass1.initialize_local_blocks(ind);
+  test(info1, ass1);
   deallog.pop();
   deallog.push("Multiple blocks");
-  test(ass1);
-  ind.push_back(2);
-  deallog.pop();
-  deallog.push("Same blocks");
-  test(ass1);
-  ass1.initialize_local_blocks(ind);
+  test(info1b, ass1);
   deallog.pop();
   deallog.push("More blocks");
-  test(ass1);  
+  test(info2b, ass1);  
 }
index 72033b0d24408c000aeb066442874c4ce0dfade1..6852e6de9056d6a3080b8e0c0f0c758cba5a37c4 100644 (file)
@@ -35,32 +35,6 @@ DEAL:Multiple blocks::  1,2 0x0 face 1,2 0x0
 DEAL:Multiple blocks::  2,0 0x0 face 2,0 0x0
 DEAL:Multiple blocks::  2,1 0x0 face 2,1 0x0
 DEAL:Multiple blocks::  2,2 0x0 face 2,2 0x0
-DEAL:Same blocks::No faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 0
-DEAL:Same blocks::  0,0 0x0
-DEAL:Same blocks::  0,1 0x0
-DEAL:Same blocks::  0,2 0x0
-DEAL:Same blocks::  1,0 0x0
-DEAL:Same blocks::  1,1 0x0
-DEAL:Same blocks::  1,2 0x0
-DEAL:Same blocks::  2,0 0x0
-DEAL:Same blocks::  2,1 0x0
-DEAL:Same blocks::  2,2 0x0
-DEAL:Same blocks::With faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 9
-DEAL:Same blocks::  0,0 0x0 face 0,0 0x0
-DEAL:Same blocks::  0,1 0x0 face 0,1 0x0
-DEAL:Same blocks::  0,2 0x0 face 0,2 0x0
-DEAL:Same blocks::  1,0 0x0 face 1,0 0x0
-DEAL:Same blocks::  1,1 0x0 face 1,1 0x0
-DEAL:Same blocks::  1,2 0x0 face 1,2 0x0
-DEAL:Same blocks::  2,0 0x0 face 2,0 0x0
-DEAL:Same blocks::  2,1 0x0 face 2,1 0x0
-DEAL:Same blocks::  2,2 0x0 face 2,2 0x0
 DEAL:More blocks::No faces
 DEAL:More blocks::J: 0
 DEAL:More blocks::R: 0
index 9a62ca644abe2d8ffcd58a3abd6282b3aa4d00fe..d20bd8b56238adfa7e7588978def4574e3547b1a 100644 (file)
@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2012 by the deal.II authors
+//    Copyright (C) 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 
 /**
  * @file Test initialization of Assembler::MatrixSimple and
- * LocalResults
+ * DoFInfo
  */
 
 #include "../tests.h"
 #include <deal.II/base/logstream.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/block_indices.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/meshworker/local_results.h>
+#include <deal.II/meshworker/dof_info.h>
 #include <deal.II/meshworker/simple.h>
 
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_system.h>
+
 using namespace dealii;
 
-template <class MATRIX>
-void test(MeshWorker::Assembler::MGMatrixSimple<MATRIX>& ass)
-{
-  MeshWorker::LocalResults<double> info1;
-  ass.initialize_info(info1, false);
+template <class DOFINFO, class MATRIX>
+void test(DOFINFO& info, MeshWorker::Assembler::MGMatrixSimple<MATRIX>& ass)
+{  
+  ass.initialize_info(info, false);
   deallog << "No faces" << std::endl;
-  info1.print_debug(deallog);
+  info.print_debug(deallog);
   
-  MeshWorker::LocalResults<double> info2;
-  ass.initialize_info(info2, true);
+  ass.initialize_info(info, true);
   deallog << "With faces" << std::endl;
-  info2.print_debug(deallog);
+  info.print_debug(deallog);
 }
 
 int main()
@@ -45,23 +50,32 @@ int main()
   deallog.attach(logfile);
   deallog.depth_console (0);
   
+  Triangulation<2,2> tr;
+  GridGenerator::hyper_cube(tr);
+  FE_DGP<2,2> fe1(1);
+  FE_DGP<2,2> fe2(2);
+  FE_DGP<2,2> fe3(3);
+  FE_DGP<2,2> fe5(5);
+  FESystem<2,2> fes1(fe3,1,fe5,1,fe1,1);
+  FESystem<2,2> fes2(fe3,1,fe5,1,fe1,1, fe2,1);
+
+  DoFHandler<2,2> dof1(tr);
+  dof1.distribute_dofs(fes1);
+  DoFHandler<2,2> dof2(tr);
+  dof2.distribute_dofs(fes2);    
+  dof1.initialize_local_block_info();
+  dof2.initialize_local_block_info();
+  MeshWorker::DoFInfo<2,2,double> info1(dof1);
+  MeshWorker::DoFInfo<2,2,double> info1b(dof1.block_info());
+  MeshWorker::DoFInfo<2,2,double> info2b(dof2.block_info());
+
   MeshWorker::Assembler::MGMatrixSimple<FullMatrix<double> > ass1;
   deallog.push("Single block");
-  test(ass1);
-  BlockIndices ind;
-  ind.push_back(3);
-  ind.push_back(5);
-  ind.push_back(1);
-  ass1.initialize_local_blocks(ind);
+  test(info1, ass1);
   deallog.pop();
   deallog.push("Multiple blocks");
-  test(ass1);
-  ind.push_back(2);
-  deallog.pop();
-  deallog.push("Same blocks");
-  test(ass1);
-  ass1.initialize_local_blocks(ind);
+  test(info1b, ass1);
   deallog.pop();
   deallog.push("More blocks");
-  test(ass1);  
+  test(info2b, ass1);  
 }
index 72033b0d24408c000aeb066442874c4ce0dfade1..6852e6de9056d6a3080b8e0c0f0c758cba5a37c4 100644 (file)
@@ -35,32 +35,6 @@ DEAL:Multiple blocks::  1,2 0x0 face 1,2 0x0
 DEAL:Multiple blocks::  2,0 0x0 face 2,0 0x0
 DEAL:Multiple blocks::  2,1 0x0 face 2,1 0x0
 DEAL:Multiple blocks::  2,2 0x0 face 2,2 0x0
-DEAL:Same blocks::No faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 0
-DEAL:Same blocks::  0,0 0x0
-DEAL:Same blocks::  0,1 0x0
-DEAL:Same blocks::  0,2 0x0
-DEAL:Same blocks::  1,0 0x0
-DEAL:Same blocks::  1,1 0x0
-DEAL:Same blocks::  1,2 0x0
-DEAL:Same blocks::  2,0 0x0
-DEAL:Same blocks::  2,1 0x0
-DEAL:Same blocks::  2,2 0x0
-DEAL:Same blocks::With faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 9
-DEAL:Same blocks::  0,0 0x0 face 0,0 0x0
-DEAL:Same blocks::  0,1 0x0 face 0,1 0x0
-DEAL:Same blocks::  0,2 0x0 face 0,2 0x0
-DEAL:Same blocks::  1,0 0x0 face 1,0 0x0
-DEAL:Same blocks::  1,1 0x0 face 1,1 0x0
-DEAL:Same blocks::  1,2 0x0 face 1,2 0x0
-DEAL:Same blocks::  2,0 0x0 face 2,0 0x0
-DEAL:Same blocks::  2,1 0x0 face 2,1 0x0
-DEAL:Same blocks::  2,2 0x0 face 2,2 0x0
 DEAL:More blocks::No faces
 DEAL:More blocks::J: 0
 DEAL:More blocks::R: 0
diff --git a/tests/integrators/elasticity_01.cc b/tests/integrators/elasticity_01.cc
new file mode 100644 (file)
index 0000000..b3234e4
--- /dev/null
@@ -0,0 +1,242 @@
+//--------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2012, 2013 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 the functions in integrators/elasticity.h
+// Output matrices and assert consistency of residuals
+
+#include "../tests.h"
+#include "../lib/test_grids.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+
+using namespace LocalIntegrators::Laplace;
+
+template <int dim>
+void test_cell(const FEValuesBase<dim>& fev)
+{
+  const unsigned int n = fev.dofs_per_cell;
+  FullMatrix<double> M(n,n);
+  cell_matrix(M,fev);
+  M.print(deallog,8);
+  
+  Vector<double> u(n), v(n), w(n);
+  std::vector<std::vector<Tensor<1,dim> > >
+    ugrad(dim,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+  
+  std::vector<unsigned int> indices(n);
+  for (unsigned int i=0;i<n;++i)
+    indices[i] = i;
+  
+  deallog << "Residuals";
+  for (unsigned int i=0;i<n;++i)
+    {
+      u = 0.;
+      u(i) = 1.;
+      w = 0.;
+      fev.get_function_gradients(u, indices, ugrad, true);
+      cell_residual(w, fev, make_slice(ugrad));
+      M.vmult(v,u);
+      w.add(-1., v);
+      deallog << ' ' << w.l2_norm();
+    }
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void test_boundary(const FEValuesBase<dim>& fev)
+{
+  const unsigned int n = fev.dofs_per_cell;
+  unsigned int d=fev.get_fe().n_components();
+  FullMatrix<double> M(n,n);
+  nitsche_matrix(M, fev, 17);
+  M.print(deallog,8);
+  
+  Vector<double> u(n), v(n), w(n);
+  std::vector<std::vector<double> >
+    uval    (d,std::vector<double>(fev.n_quadrature_points)),
+    null_val(d,std::vector<double>(fev.n_quadrature_points, 0.));
+  std::vector<std::vector<Tensor<1,dim> > >
+    ugrad   (d,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+  
+  std::vector<unsigned int> indices(n);
+  for (unsigned int i=0;i<n;++i)
+    indices[i] = i;
+  
+  deallog << "Residuals";
+  for (unsigned int i=0;i<n;++i)
+    {
+      u = 0.;
+      u(i) = 1.;
+      w = 0.;
+      fev.get_function_values(u, indices, uval, true);
+      fev.get_function_gradients(u, indices, ugrad, true);
+      nitsche_residual(w, fev, make_slice(uval), make_slice(ugrad), make_slice(null_val), 17);
+      M.vmult(v,u);
+      w.add(-1., v);
+      deallog << ' ' << w.l2_norm();
+      if (d==1)
+       {
+         nitsche_residual(w, fev, uval[0], ugrad[0], null_val[0], 17);
+         M.vmult(v,u);
+         w.add(-1., v);
+         deallog << " e" << w.l2_norm();         
+       }
+    }
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void test_face(const FEValuesBase<dim>& fev1,
+              const FEValuesBase<dim>& fev2)
+{
+  const unsigned int n1 = fev1.dofs_per_cell;
+  const unsigned int n2 = fev1.dofs_per_cell;
+  unsigned int d=fev1.get_fe().n_components();
+  FullMatrix<double> M11(n1,n1);
+  FullMatrix<double> M12(n1,n2);
+  FullMatrix<double> M21(n2,n1);
+  FullMatrix<double> M22(n2,n2);
+  
+  ip_matrix(M11, M12, M21, M22, fev1, fev2, 17);
+  deallog << "M11" << std::endl;
+  M11.print(deallog,8);
+  deallog << "M22" << std::endl;
+  M22.print(deallog,8);
+  deallog << "M12" << std::endl;
+  M12.print(deallog,8);
+  deallog << "M21" << std::endl;
+  M21.print(deallog,8);
+  
+  Vector<double> u1(n1), v1(n1), w1(n1);
+  Vector<double> u2(n2), v2(n2), w2(n2);
+  std::vector<std::vector<double> >
+    u1val (d,std::vector<double>(fev1.n_quadrature_points)),
+    u2val (d,std::vector<double>(fev2.n_quadrature_points));
+  std::vector<std::vector<Tensor<1,dim> > >
+    u1grad(d,std::vector<Tensor<1,dim> >(fev1.n_quadrature_points)),
+    u2grad(d,std::vector<Tensor<1,dim> >(fev2.n_quadrature_points));
+  
+  std::vector<unsigned int> indices1(n1), indices2(n2);
+  for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+  for (unsigned int i=0;i<n2;++i) indices2[i] = i;
+  
+  deallog << "Residuals";  for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+
+  for (unsigned int i1=0;i1<n1;++i1)
+    for (unsigned int i2=0;i2<n2;++i2)
+      {
+       u1 = 0.;
+       u1(i1) = 1.;
+       w1 = 0.;
+       fev1.get_function_values(u1, indices1, u1val, true);
+       fev1.get_function_gradients(u1, indices1, u1grad, true);
+       u2 = 0.;
+       u2(i2) = 1.;
+       w2 = 0.;
+       fev2.get_function_values(u2, indices2, u2val, true);
+       fev2.get_function_gradients(u2, indices2, u2grad, true);
+       ip_residual(w1, w2, fev1, fev2,
+                   make_slice(u1val), make_slice(u1grad),
+                   make_slice(u2val), make_slice(u2grad),
+                   17);
+       M11.vmult(v1,u1); w1.add(-1., v1);
+       M12.vmult(v1,u2); w1.add(-1., v1);
+       M21.vmult(v2,u1); w2.add(-1., v2);
+       M22.vmult(v2,u2); w2.add(-1., v2);
+       deallog << ' ' << w1.l2_norm() + w2.l2_norm();
+       if (d==1)
+         {
+           ip_residual(w1, w2, fev1, fev2, u1val[0], u1grad[0], u2val[0], u2grad[0], 17);
+           M11.vmult(v1,u1); w1.add(-1., v1);
+           M12.vmult(v1,u2); w1.add(-1., v1);
+           M21.vmult(v2,u1); w2.add(-1., v2);
+           M22.vmult(v2,u2); w2.add(-1., v2);
+           deallog << " e" << w1.l2_norm() + w2.l2_norm();       
+         }
+      }
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void
+test_fe(Triangulation<dim>& tr, FiniteElement<dim>& fe)
+{
+  deallog << fe.get_name() << std::endl << "cell matrix" << std::endl;
+  QGauss<dim> quadrature(fe.tensor_degree()+1);
+  FEValues<dim> fev(fe, quadrature, update_gradients);
+
+  typename Triangulation<dim>::cell_iterator cell1 = tr.begin(1);
+  fev.reinit(cell1);
+  test_cell(fev);
+  
+  QGauss<dim-1> face_quadrature(fe.tensor_degree()+1);
+  FEFaceValues<dim> fef1(fe, face_quadrature, update_values | update_gradients | update_normal_vectors);
+  for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+    {
+      deallog << "boundary_matrix " << i << std::endl;
+      fef1.reinit(cell1, i);
+      test_boundary(fef1);
+    }
+
+  FEFaceValues<dim> fef2(fe, face_quadrature, update_values | update_gradients);
+  typename Triangulation<dim>::cell_iterator cell2 = cell1->neighbor(1);
+  
+  deallog << "face_matrix " << 0 << std::endl;
+  cell1 = tr.begin(1);
+  fef1.reinit(cell1, 1);
+  fef2.reinit(cell2, 0);
+  test_face(fef1, fef2);  
+}
+
+
+template <int dim>
+void
+test(Triangulation<dim>& tr)
+{
+  FE_DGQ<dim> q1(1);
+  FESystem<dim> fe1(q1,dim);
+  test_fe(tr, fe1);
+  
+  FE_DGQ<dim> q2(2);
+  FESystem<dim> fe2(q2,dim);
+  test_fe(tr, fe2);
+  
+  FE_Nedelec<dim> n1(1);  
+  test_fe(tr, n1);  
+  
+  FE_RaviartThomas<dim> rt1(1);  
+  test_fe(tr, rt1);  
+}
+
+
+int main()
+{
+  initlog(__FILE__);
+  deallog.threshold_double(1.e-10);
+  
+  Triangulation<2> tr2;
+  TestGrids::hypercube(tr2, 1);
+  test(tr2);
+  
+ Triangulation<3> tr3;
+ TestGrids::hypercube(tr3, 1);
+ test(tr3);  
+}

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.