run into trouble

git-svn-id: https://svn.dealii.org/trunk@30385 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/multigrid/step-16-02.cc b/tests/multigrid/step-16-02.cc
index 01bd6949bb11e12c8302ecf571d5baeb3b5ebc16..2c3283c5ae285bb689a2e66aedc15209eabd2b80 100644 (file)
--- a/tests/multigrid/step-16-02.cc
+++ b/tests/multigrid/step-16-02.cc
@@ -42,6 +42,8 @@
 #include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/tria_boundary_lib.h>
 
+#include <deal.II/integrators/laplace.h>
+
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_tools.h>
 
@@ -60,10 +62,64 @@
 #include <deal.II/multigrid/mg_smoother.h>
 #include <deal.II/multigrid/mg_matrix.h>
 
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/simple.h>
+#include <deal.II/meshworker/output.h>
+#include <deal.II/meshworker/loop.h>
+
 #include <fstream>
 #include <sstream>
 
 using namespace dealii;
+using namespace LocalIntegrators;
+
+template <int dim>
+class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  LaplaceMatrix();
+  virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+  virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+  virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
+                   MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const;
+};
+
+
+template <int dim>
+LaplaceMatrix<dim>::LaplaceMatrix()
+{}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
+{
+  AssertDimension (dinfo.n_matrices(), 1);  
+  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0));
+}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::boundary(MeshWorker::DoFInfo<dim>& dinfo,
+                                      typename MeshWorker::IntegrationInfo<dim>& info) const
+{
+  const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
+  Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0),
+                         Laplace::compute_penalty(dinfo, dinfo, deg, deg));
+}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::face(
+  MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
+  MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const
+{
+  const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
+  Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, 
+                    dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
+                    info1.fe_values(0), info2.fe_values(0),
+                    Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+}
 
 template <int dim>
 class LaplaceProblem
@@ -75,7 +131,7 @@ public:
 private:
   void setup_system ();
   void assemble_system ();
-  void assemble_multigrid ();
+  void assemble_multigrid (const bool &use_mw);
   void solve ();
   void refine_grid ();
   void output_results (const unsigned int cycle) const;
@@ -94,10 +150,12 @@ private:
   Vector<double>       system_rhs;
 
   const unsigned int degree;
+  LaplaceMatrix<dim> matrix_integrator;
 
   MGLevelObject<SparsityPattern>       mg_sparsity_patterns;
   MGLevelObject<SparseMatrix<double> > mg_matrices;
-  MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
+  MGLevelObject<SparseMatrix<double> > mg_interface_in;
+  MGLevelObject<SparseMatrix<double> > mg_interface_out;
   MGConstrainedDoFs                    mg_constrained_dofs;
 };
 
@@ -155,7 +213,8 @@ LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
                  limit_level_difference_at_vertices),
   fe (degree),
   mg_dof_handler (triangulation),
-  degree(degree)
+  degree(degree),
+  matrix_integrator()
 {}
 
 
@@ -203,8 +262,10 @@ void LaplaceProblem<dim>::setup_system ()
   mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary);
   const unsigned int n_levels = triangulation.n_levels();
 
-  mg_interface_matrices.resize(0, n_levels-1);
-  mg_interface_matrices.clear ();
+  mg_interface_in.resize(0, n_levels-1);
+  mg_interface_in.clear ();
+  mg_interface_out.resize(0, n_levels-1);
+  mg_interface_out.clear ();
   mg_matrices.resize(0, n_levels-1);
   mg_matrices.clear ();
   mg_sparsity_patterns.resize(0, n_levels-1);
@@ -219,7 +280,8 @@ void LaplaceProblem<dim>::setup_system ()
       mg_sparsity_patterns[level].copy_from (csp);
 
       mg_matrices[level].reinit(mg_sparsity_patterns[level]);
-      mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]);
+      mg_interface_in[level].reinit(mg_sparsity_patterns[level]);
+      mg_interface_out[level].reinit(mg_sparsity_patterns[level]);
     }
 }
 
@@ -280,124 +342,122 @@ void LaplaceProblem<dim>::assemble_system ()
 
 
 template <int dim>
-void LaplaceProblem<dim>::assemble_multigrid ()
+void LaplaceProblem<dim>::assemble_multigrid (const bool& use_mw)
 {
-  QGauss<dim>  quadrature_formula(1+degree);
+  if(use_mw == true)
+  {
+    mg_matrices = 0.;
+
+    MappingQ1<dim> mapping;
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags update_flags = update_values | update_gradients | update_hessians;
+    info_box.add_update_flags_all(update_flags);
+    info_box.initialize(fe, mapping);//, &dof_handler.block_info());
+
+    MeshWorker::DoFInfo<dim> dof_info(mg_dof_handler);//.block_info());
+
+    MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
+    assembler.initialize(mg_constrained_dofs);
+    //  assembler.initialize_local_blocks(dof_handler.block_info().local());
+    assembler.initialize(mg_matrices);
+    assembler.initialize_interfaces(mg_interface_in, mg_interface_out);
+//    assembler.initialize_fluxes(mg_matrix_up, mg_matrix_down);
+
+    MeshWorker::integration_loop<dim, dim> (
+        mg_dof_handler.begin(), mg_dof_handler.end(),
+        dof_info, info_box, matrix_integrator, assembler);
+
+    const unsigned int nlevels = triangulation.n_levels();
+    for (unsigned int level=0;level<nlevels;++level)
+    {
+      for(unsigned int i=0; i<mg_dof_handler.n_dofs(level); ++i)
+        if(mg_matrices[level].diag_element(i)==0)
+          mg_matrices[level].set(i,i,1.);
+    }
+  }
+  else
+  {
+    QGauss<dim>  quadrature_formula(1+degree);
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                           update_values   | update_gradients |
-                           update_quadrature_points | update_JxW_values);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+        update_values   | update_gradients |
+        update_quadrature_points | update_JxW_values);
 
-  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-  const unsigned int   n_q_points      = quadrature_formula.size();
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int   n_q_points      = quadrature_formula.size();
 
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
-  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values (n_q_points);
+    const Coefficient<dim> coefficient;
+    std::vector<double>    coefficient_values (n_q_points);
 
-  std::vector<std::vector<bool> > interface_dofs
-    = mg_constrained_dofs.get_refinement_edge_indices ();
-  std::vector<std::vector<bool> > boundary_interface_dofs
-    = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
+    std::vector<std::vector<bool> > interface_dofs
+      = mg_constrained_dofs.get_refinement_edge_indices ();
+    std::vector<std::vector<bool> > boundary_interface_dofs
+      = mg_constrained_dofs.get_refinement_edge_boundary_indices ();
 
-  std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
-  std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
-  for (unsigned int level=0; level<triangulation.n_levels(); ++level)
+    std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_levels());
+    std::vector<ConstraintMatrix> boundary_interface_constraints (triangulation.n_levels());
+    for (unsigned int level=0; level<triangulation.n_levels(); ++level)
     {
       boundary_constraints[level].add_lines (interface_dofs[level]);
       boundary_constraints[level].add_lines (mg_constrained_dofs.get_boundary_indices()[level]);
       boundary_constraints[level].close ();
 
       boundary_interface_constraints[level]
-      .add_lines (boundary_interface_dofs[level]);
+        .add_lines (boundary_interface_dofs[level]);
       boundary_interface_constraints[level].close ();
     }
 
-  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                            endc = mg_dof_handler.end();
+    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+             endc = mg_dof_handler.end();
 
-  for (; cell!=endc; ++cell)
+    for (; cell!=endc; ++cell)
     {
       cell_matrix = 0;
       fe_values.reinit (cell);
 
       coefficient.value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
+          coefficient_values);
 
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           for (unsigned int j=0; j<dofs_per_cell; ++j)
             cell_matrix(i,j) += (coefficient_values[q_point] *
-                                 fe_values.shape_grad(i,q_point) *
-                                 fe_values.shape_grad(j,q_point) *
-                                 fe_values.JxW(q_point));
+                fe_values.shape_grad(i,q_point) *
+                fe_values.shape_grad(j,q_point) *
+                fe_values.JxW(q_point));
 
       cell->get_mg_dof_indices (local_dof_indices);
 
       boundary_constraints[cell->level()]
-      .distribute_local_to_global (cell_matrix,
-                                   local_dof_indices,
-                                   mg_matrices[cell->level()]);
-
-      // The next step is again slightly more
-      // obscure (but explained in the @ref
-      // mg_paper): We need the remainder of
-      // the operator that we just copied
-      // into the <code>mg_matrices</code>
-      // object, namely the part on the
-      // interface between cells at the
-      // current level and cells one level
-      // coarser. This matrix exists in two
-      // directions: for interior DoFs (index
-      // $i$) of the current level to those
-      // sitting on the interface (index
-      // $j$), and the other way around. Of
-      // course, since we have a symmetric
-      // operator, one of these matrices is
-      // the transpose of the other.
-      //
-      // The way we assemble these matrices
-      // is as follows: since the are formed
-      // from parts of the local
-      // contributions, we first delete all
-      // those parts of the local
-      // contributions that we are not
-      // interested in, namely all those
-      // elements of the local matrix for
-      // which not $i$ is an interface DoF
-      // and $j$ is not. The result is one of
-      // the two matrices that we are
-      // interested in, and we then copy it
-      // into the
-      // <code>mg_interface_matrices</code>
-      // object. The
-      // <code>boundary_interface_constraints</code>
-      // object at the same time makes sure
-      // that we delete contributions from
-      // all degrees of freedom that are not
-      // only on the interface but also on
-      // the external boundary of the domain.
-      //
-      // The last part to remember is how to
-      // get the other matrix. Since it is
-      // only the transpose, we will later
-      // (in the <code>solve()</code>
-      // function) be able to just pass the
-      // transpose matrix where necessary.
+        .distribute_local_to_global (cell_matrix,
+            local_dof_indices,
+            mg_matrices[cell->level()]);
+
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           if ( !(interface_dofs[cell->level()][local_dof_indices[i]]==true &&
-                 interface_dofs[cell->level()][local_dof_indices[j]]==false))
+                interface_dofs[cell->level()][local_dof_indices[j]]==false))
             cell_matrix(i,j) = 0;
 
       boundary_interface_constraints[cell->level()]
-      .distribute_local_to_global (cell_matrix,
-                                   local_dof_indices,
-                                   mg_interface_matrices[cell->level()]);
+        .distribute_local_to_global (cell_matrix,
+            local_dof_indices,
+            mg_interface_in[cell->level()]);
     }
+  }
+
+  const unsigned int nlevels = triangulation.n_levels();
+  for (unsigned int level=1;level<nlevels;++level)
+  {
+//    deallog << "dG up " << mg_matrix_up[level].l1_norm();
+//    deallog << " dG down " << mg_matrix_down[level].l1_norm() << std::endl;
+    deallog << "cG in " << mg_interface_in[level].l1_norm();
+    deallog << " cG out " << mg_interface_out[level].l1_norm() << std::endl;
+  }
 }
 
 
@@ -455,8 +515,8 @@ void LaplaceProblem<dim>::solve ()
   mg_smoother.set_symmetric(true);
 
   MGMatrix<> mg_matrix(&mg_matrices);
-  MGMatrix<> mg_interface_up(&mg_interface_matrices);
-  MGMatrix<> mg_interface_down(&mg_interface_matrices);
+  MGMatrix<> mg_interface_up(&mg_interface_in);
+  MGMatrix<> mg_interface_down(&mg_interface_out);
 
   Multigrid<Vector<double> > mg(mg_dof_handler,
                                 mg_matrix,
@@ -586,7 +646,7 @@ void LaplaceProblem<dim>::run ()
       deallog << std::endl;
 
       assemble_system ();
-      assemble_multigrid ();
+      assemble_multigrid (true);
 
       solve ();
 //      output_results (cycle);