]> https://gitweb.dealii.org/ - dealii.git/commitdiff
First test for integrators... failing due to program crashing in constructor of FE_Ne...
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 9 Nov 2010 15:45:13 +0000 (15:45 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 9 Nov 2010 15:45:13 +0000 (15:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@22645 0785d39b-7218-0410-832d-ea1e28bc413d

tests/integrators/Makefile [new file with mode: 0644]
tests/integrators/cochain_01.cc [new file with mode: 0644]
tests/lib/test_grids.h

diff --git a/tests/integrators/Makefile b/tests/integrators/Makefile
new file mode 100644 (file)
index 0000000..6476e92
--- /dev/null
@@ -0,0 +1,10 @@
+############################################################
+# $Id$
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+############################################################
+
+include ../Makefile.paths
+include $D/common/Make.global_options
+include ../Makefile.rules
+include Makefile.depend
+include Makefile.tests
diff --git a/tests/integrators/cochain_01.cc b/tests/integrators/cochain_01.cc
new file mode 100644 (file)
index 0000000..f410187
--- /dev/null
@@ -0,0 +1,273 @@
+//--------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2005, 2006, 2010 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, whether differential operators produce a cochain complex on
+// the standard Hilbert space sequence
+
+#include "../tests.h"
+#include "../lib/test_grids.h"
+
+#include <base/logstream.h>
+
+#include <lac/matrix_block.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <numerics/mesh_worker_info.h>
+#include <numerics/mesh_worker_assembler.h>
+#include <numerics/mesh_worker_loop.h>
+
+#include <integrators/l2.h>
+#include <integrators/differential.h>
+#include <integrators/laplace.h>
+#include <integrators/maxwell.h>
+#include <integrators/elasticity.h>
+
+using namespace LocalIntegrators;
+
+const bool debugging = false;
+
+template <int dim>
+void cell_matrix(
+  MeshWorker::DoFInfo<dim>& dinfo,
+  typename MeshWorker::IntegrationInfo<dim,dim>& info)
+{
+  unsigned int dm=0; // Matrix index
+  unsigned int de=0; // Element index
+  
+  L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+  Laplace::cell_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+  Differential::gradient_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1));
+
+  
+  ++de;
+  L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+  Maxwell::curl_curl_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+  Maxwell::curl_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1));
+
+  if (dim>2)
+    {
+      ++de;
+      L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+      Elasticity::grad_div_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+      Differential::divergence_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1));
+    }
+  
+  ++de;
+  L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+}
+
+
+
+template <int dim>
+void
+test_cochain(const Triangulation<dim>& tr, const FiniteElement<dim>& fe)
+{
+  MappingQ1<dim> mapping;
+                                  // Initialize DofHandler for a
+                                  // block system with local blocks
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+  DoFRenumbering::component_wise(dof);
+  dof.initialize_local_block_info();
+
+                                  // Initialize hanging node constraints
+  ConstraintMatrix constraints;
+//  DoFTools::make_zero_boundary_constraints(dof, constraints);
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  constraints.close();
+                                  // Setup sparsity pattern for the
+                                  // whole system
+  BlockSparsityPattern sparsity;
+  sparsity.reinit(dof.block_info().global().size(),
+                 dof.block_info().global().size());
+  BlockCompressedSparsityPattern c_sparsity(dof.block_info().global(),
+                                           dof.block_info().global());
+  DoFTools::make_flux_sparsity_pattern(dof, c_sparsity, constraints, false);
+  sparsity.copy_from(c_sparsity);
+
+                                  // Setup matrices
+  MatrixBlockVector<SparseMatrix<double> > matrices;
+
+  unsigned int d=0;
+  matrices.add(d,d,"mass-H1");
+  matrices.add(d,d,"Laplacian");
+  matrices.add(d+1,d,"grad");
+  
+  ++d;
+  matrices.add(d,d,"mass-Hcurl");
+  matrices.add(d,d,"curl-curl");
+  matrices.add(d+1,d,"curl");
+  
+  if (dim>2)
+    {
+      ++d;
+      matrices.add(d,d,"mass-Hdiv");
+      matrices.add(d,d,"grad-div");
+      matrices.add(d+1,d,"div");
+    }
+
+  ++d;
+  matrices.add(d,d,"mass-L2");
+
+  matrices.reinit(sparsity);
+
+  if (debugging)
+    for (unsigned int i=0;i<matrices.size();++i)
+      deallog << "Block " << '(' << matrices.block(i).row
+             << ',' << matrices.block(i).column << ") "
+             << std::setw(3) << std::right << matrices.matrix(i).m() << std::left
+             << 'x' << std::setw(3) << matrices.matrix(i).n()
+             << ' ' << matrices.name(i) << std::endl;
+  
+                                  // Build matrices
+  
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  UpdateFlags update_flags = update_values | update_gradients;
+  info_box.add_update_flags(update_flags, true, true, true, true);
+  info_box.initialize(fe, mapping, &dof.block_info());
+  
+  MeshWorker::DoFInfo<dim> dof_info(dof.block_info());
+  
+  MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks<SparseMatrix<double> > assembler;
+  assembler.initialize(&dof.block_info(), matrices);
+  assembler.initialize(constraints);
+  
+  MeshWorker::integration_loop<dim, dim>(
+    dof.begin_active(), dof.end(), dof_info, info_box,
+    cell_matrix<dim>, 0, 0, assembler);
+  
+  for (unsigned int b=0;b<matrices.size();++b)
+    if (b%3 == 0)
+      for (unsigned int i=0;i<matrices.block(b).matrix.m();++i)
+       if (matrices.block(b).matrix.diag_element(i) == 0.)
+         matrices.block(b).matrix.diag_element(i) = 1.;
+
+                                  // Set up vectors
+  BlockVector<double> source(dof.block_info().global());
+  BlockVector<double> result1(dof.block_info().global());
+  BlockVector<double> result2(dof.block_info().global());
+  BlockVector<double> aux(dof.block_info().global());
+  for (unsigned int i=0;i<source.size();++i)
+    source(i) = i%5;
+
+                                  // now check, whether d*d =
+                                  // D^TM^{-1}D
+  SolverControl control(100, 1.e-13, false, false);
+  SolverCG<Vector<double> > solver(control);
+  
+  for (unsigned int d=0;d<dim;++d)
+    {
+      deallog << "Form " << d << std::endl;
+      const unsigned int m=3*d;
+
+      if (d>0)
+       {
+         matrices.matrix(m+2).vmult(result2.block(d+1), aux.block(d));
+         deallog << "d^2            " << result2.block(d+1).l2_norm() << std::endl;  
+         matrices.matrix(m+1).vmult(result2.block(d), aux.block(d));
+         deallog << "d*d^2          " << result2.block(d).l2_norm() << std::endl;  
+       }
+      
+      matrices.matrix(m+2).vmult(result1.block(d+1), source.block(d));      
+
+      solver.solve(matrices.matrix(m+3), aux.block(d+1), result1.block(d+1), PreconditionIdentity());
+      
+      matrices.matrix(m+2).Tvmult(result1.block(d), aux.block(d+1));
+      matrices.matrix(m+1).vmult(result2.block(d), source.block(d));
+
+      if (debugging)
+       deallog << "u " << source.block(d).l2_norm() << ' '
+               << matrices.name(m+2) << ' ' << result1.block(d+1).l2_norm()
+               << ' ' << matrices.name(m+3) << ' ' << aux.block(d+1).l2_norm()
+               << ' ' << matrices.name(m+2) << "^T " << result1.block(d).l2_norm()
+               << ' ' << matrices.name(m+1) << ' ' << result2.block(d).l2_norm() << std::endl;
+      result2.block(d) -= result1.block(d);
+      deallog << "Difference d*d " << result2.block(d).l2_norm() << std::endl;
+    }
+}
+
+void run (const Triangulation<2>& tr, unsigned int degree)
+{
+  std::ostringstream prefix;
+  prefix << "d2-p" << degree;
+  deallog.push(prefix.str());
+  deallog << "Setup" << std::endl;
+  
+  FE_Q<2> h1(degree+1);
+  FE_Nedelec<2> hdiv(degree);
+  FE_DGQ<2> l2(degree);
+
+  FESystem<2> fe(h1,1,hdiv,1,l2,1);
+
+  test_cochain(tr, fe);
+  deallog.pop();
+}
+
+
+void run (const Triangulation<3>& tr, unsigned int degree)
+{
+  std::ostringstream prefix;
+  prefix << "d3-p" << degree;
+  deallog.push(prefix.str());
+  deallog << "Setup" << std::endl;
+  
+  FE_Q<3> h1(degree+1);
+  if (debugging) deallog << "H1" << std::endl;
+  FE_Nedelec<3> hcurl(degree);
+  if (debugging) deallog << "Hcurl" << std::endl;
+  FE_RaviartThomas<3> hdiv(degree);
+  if (debugging) deallog << "Hdiv" << std::endl;
+  FE_DGQ<3> l2(degree);
+  if (debugging) deallog << "L2" << std::endl;
+
+  FESystem<3> fe(h1,1,hcurl,1,hdiv,1,l2,1);
+
+  test_cochain(tr, fe);
+  deallog.pop();
+}
+
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+//  deallog.depth_console(0);
+//  deallog.threshold_double(1.e-10);
+
+  Triangulation<2> tr2;
+  TestGrids::hypercube(tr2, 1);
+  run(tr2, 0);
+  run(tr2, 1);
+  run(tr2, 2);
+  
+  Triangulation<3> tr3;
+  TestGrids::hypercube(tr3, 1);
+  run(tr3, 0);
+  run(tr3, 1);
+}
index dcd1564c3c1b68dbfeb39edc683511d27f12fbee..717435e10a2245ca80638d07e9938686cfed9f36 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2006 by the deal.II authors
+//    Copyright (C) 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -12,6 +12,7 @@
 //----------------------------------------------------------------------
 
 #include <grid/tria.h>
+#include <grid/grid_generator.h>
 
 /**
  * A set of test meshes for the deal.II test suite.
@@ -34,7 +35,8 @@
  * level</td></tr>
  * </table>
  *
- * @author Guido Kanschat, 2006
+ * @author Guido Kanschat
+ * @date 2006, 2010
  */
 namespace TestGrids
 {
@@ -55,7 +57,24 @@ namespace TestGrids
   template <int dim>
   void hypercube(Triangulation<dim>& tr,
                 unsigned int refinement = 0,
-                bool local = false);
+                bool local = false)
+  {
+    GridGenerator::hyper_cube(tr);
+    if (refinement && !local)
+      tr.refine_global(refinement);
+    if (refinement && local)
+      {
+       for (typename Triangulation<dim>::active_cell_iterator
+              cell = tr.begin_active(); cell != tr.end(); ++cell)
+         {
+           const Point<dim>& p = cell->center();
+           bool negative = true;
+           for (unsigned int d=0;d<dim;++d)
+             if (p(d) >= 0.)negative = false;
+         }
+      }
+  }
+  
                                   /**
                                    * Create a star-shaped mesh,
                                    * having more than the average
@@ -91,26 +110,4 @@ namespace TestGrids
                                    */
   template <int dim>
   void laguna(Triangulation<dim>& tr);
-
-
-  template <int dim>
-  void hypercube(Triangulation<dim>& tr,
-                unsigned int refinement,
-                bool local)
-  {
-    GridGenerator::hyper_cube(tr);
-    if (!local)
-      {
-       if (refinement > 0)
-         tr.refine_global(refinement);
-      }
-    else
-      {
-       for (unsigned int i=0;i<refinement;++i)
-         {
-           tr.begin_active()->set_refine_flag();
-           tr.execute_coarsening_and_refinement();
-         }
-      }
-  }
 }

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.