]> https://gitweb.dealii.org/ - dealii.git/commitdiff
use new logic in multigrid assemblers and switch step 16 to MeshWorker
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 31 Jul 2014 12:35:40 +0000 (14:35 +0200)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Thu, 31 Jul 2014 17:33:16 +0000 (19:33 +0200)
examples/step-16/step-16.cc
include/deal.II/meshworker/assembler.h
include/deal.II/meshworker/simple.h

index 32aa7d2bdeaecdd9a5b2817ea6833da844d28165..9f4e4c67f12912ecdd66563f4b6c5548c1e8867f 100644 (file)
 #include <deal.II/multigrid/mg_smoother.h>
 #include <deal.II/multigrid/mg_matrix.h>
 
+// Finally we include the MeshWorker framework. Since we have to build
+// several matrices and have to be aware of several sets of
+// constraints, we do not program loops over cells ourselves, but
+// rather leave the actual logic to MeshWorker::loop().
+#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>
+
+// In order to save effort, we use the pre-implemented Laplacian found in
+#include <deal.II/integrators/laplace.h>
+#include <deal.II/integrators/l2.h>
+
 // This is C++:
 #include <fstream>
 #include <sstream>
 
-// The last step is as in all previous programs:
+// Lazy as we are, we avoid typing namespace names
+
+using namespace dealii;
+using namespace LocalIntegrators;
+
 namespace Step16
 {
-  using namespace dealii;
+  // @{sect3}{The integrator on each cell}
+
+  // MeshWorker::integration_loop() expects a class that provides
+  // functions for integration on cells and boundary and interior
+  // faces. This is done by the following class. In the constructor,
+  // we tell the loop that cell integrals should be computed (the
+  // 'true'), but integrals should not be computed on boundary and
+  // interior faces (the two 'false').
+
+  template <int dim>
+  class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
+  {
+public:
+LaplaceMatrix();
+virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+};
+
+
+template <int dim>
+LaplaceMatrix<dim>::LaplaceMatrix()
+               :
+               MeshWorker::LocalIntegrator<dim>(true, false, false)
+{}
+
+
+template <int dim>
+void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
+{
+  AssertDimension (dinfo.n_matrices(), 1);
+const double coefficient = (dinfo.cell->center()(0) > 0.)
+                          ? .1 : 1.;
+ Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), coefficient);
+
+  if (dinfo.n_vectors() > 0)
+    {
+std::vector<double> rhs(info.fe_values(0).n_quadrature_points, 1.);
+ L2::L2(dinfo.vector(0).block(0), info.fe_values(0), rhs);
+}
+}
 
 
   // @sect3{The <code>LaplaceProblem</code> class template}
@@ -151,63 +208,12 @@ namespace Step16
     // to build one; we choose the one from coarse to fine.
     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;
   };
 
 
-
-  // @sect3{Nonconstant coefficients}
-
-  // The implementation of nonconstant coefficients is copied verbatim from
-  // step-5 and step-6:
-
-  template <int dim>
-  class Coefficient : public Function<dim>
-  {
-  public:
-    Coefficient () : Function<dim>() {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                             std::vector<double>            &values,
-                             const unsigned int              component = 0) const;
-  };
-
-
-
-  template <int dim>
-  double Coefficient<dim>::value (const Point<dim> &p,
-                                  const unsigned int) const
-  {
-    if (p.square() < 0.5*0.5)
-      return 20;
-    else
-      return 1;
-  }
-
-
-
-  template <int dim>
-  void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
-                                     std::vector<double>            &values,
-                                     const unsigned int              component) const
-  {
-    const unsigned int n_points = points.size();
-
-    Assert (values.size() == n_points,
-            ExcDimensionMismatch (values.size(), n_points));
-
-    Assert (component == 0,
-            ExcIndexRange (component, 0, 1));
-
-    for (unsigned int i=0; i<n_points; ++i)
-      values[i] = Coefficient<dim>::value (points[i]);
-  }
-
-
   // @sect3{The <code>LaplaceProblem</code> class implementation}
 
   // @sect4{LaplaceProblem::LaplaceProblem}
@@ -326,8 +332,10 @@ namespace Step16
     // upon resizing.
     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);
@@ -357,7 +365,8 @@ namespace Step16
         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]);
       }
   }
 
@@ -374,55 +383,26 @@ namespace Step16
   template <int dim>
   void LaplaceProblem<dim>::assemble_system ()
   {
-    const QGauss<dim>  quadrature_formula(degree+1);
-
-    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();
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       cell_rhs (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);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      {
-        cell_matrix = 0;
-        cell_rhs = 0;
-
-        fe_values.reinit (cell);
-
-        coefficient.value_list (fe_values.get_quadrature_points(),
-                                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));
-
-              cell_rhs(i) += (fe_values.shape_value(i,q_point) *
-                              1.0 *
-                              fe_values.JxW(q_point));
-            }
-
-        cell->get_dof_indices (local_dof_indices);
-        constraints.distribute_local_to_global (cell_matrix, cell_rhs,
-                                                local_dof_indices,
-                                                system_matrix, system_rhs);
-      }
+    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);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> > assembler;
+    assembler.initialize(constraints);
+    assembler.initialize(system_matrix, system_rhs);
+    LaplaceMatrix<dim> matrix_integrator;
+    MeshWorker::integration_loop<dim, dim> (
+        dof_handler.begin_active(), dof_handler.end(),
+        dof_info, info_box, matrix_integrator, assembler);
+
+    for(unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+      if(constraints.is_constrained(i))
+       system_matrix.set(i,i,1.);
   }
 
 
@@ -441,152 +421,31 @@ namespace Step16
   template <int dim>
   void LaplaceProblem<dim>::assemble_multigrid ()
   {
-    QGauss<dim>  quadrature_formula(1+degree);
-
-    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();
-
-    FullMatrix<double>   cell_matrix (dofs_per_cell, 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);
-
-    // Next a few things that are specific to building the multigrid data
-    // structures (since we only need them in the current function, rather
-    // than also elsewhere, we build them here instead of the
-    // <code>setup_system</code> function). Some of the following may be a bit
-    // obscure if you're not familiar with the algorithm actually implemented
-    // in deal.II to support multilevel algorithms on adaptive meshes; if some
-    // of the things below seem strange, take a look at the @ref mg_paper.
-
-    // Our first job is to identify the boundary conditions for the levels. On
-    // each level, we impose Dirichlet boundary conditions on the exterior
-    // boundary of the domain as well as on interfaces between adaptively
-    // refined levels. As in many other parts of the library, we do this by
-    // using a mask described by an IndexSet. The <code>MGConstraints</code>
-    // already computed the information for us when we called initialize in
-    // <code>setup_system()</code>. So we simply ask for them by calling
-    // <code>get_boundary_indices()</code> and
-    // <code>get_refinement_edge_indices(level)</code>. Moreover, we have to
-    // identify the subset of the refinement edge indices which are also
-    // located on the boundary as they require special treatment in the
-    // algorithm further down.
-
-    // These three masks are used to fill a ConstraintMatrix objects for each
-    // level that we use during the assembly of the matrix: the value of the
-    // associated degrees of freedom should be zero after each application of
-    // the level operators. Due to the way the ConstraintMatrix stores its
-    // data, the function to add a constraint on a single degree of freedom
-    // and force it to be zero is called Constraintmatrix::add_line(); doing
-    // so for several degrees of freedom at once can be done using
-    // Constraintmatrix::add_lines():
-    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 (mg_constrained_dofs.get_boundary_indices()[level]);
-        boundary_constraints[level].
-        add_lines (mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].close ();
-
-        boundary_interface_constraints[level]
-        .add_lines (mg_constrained_dofs.get_refinement_edge_boundary_indices(level));
-        boundary_interface_constraints[level].close ();
-      }
-
-    // Now that we're done with most of our preliminaries, let's start the
-    // integration loop. It looks mostly like the loop in
-    // <code>assemble_system</code>, with two exceptions: (i) we don't need a
-    // right hand side, and more significantly (ii) we don't just loop over
-    // all active cells, but in fact all cells, active or not. Consequently,
-    // the correct iterator to use is DoFHandler::cell_iterator rather than
-    // DoFHandler::active_cell_iterator. Let's go about it:
-    typename DoFHandler<dim>::cell_iterator cell = dof_handler.begin(),
-                                            endc = dof_handler.end();
-
-    for (; cell!=endc; ++cell)
-      {
-        cell_matrix = 0;
-        fe_values.reinit (cell);
-
-        coefficient.value_list (fe_values.get_quadrature_points(),
-                                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));
-
-        // The rest of the assembly is again slightly different. This starts
-        // with a gotcha that is easily forgotten: The indices of global
-        // degrees of freedom we want here are the ones for current level, not
-        // for the global matrix. We therefore need the function
-        // DoFAccessor::get_mg_dof_indices, not
-        // DoFAccessor::get_dof_indices as used in the assembly of the
-        // global system:
-        cell->get_mg_dof_indices (local_dof_indices);
-
-        // Next, we need to copy local contributions into the level
-        // objects. We can do this in the same way as in the global assembly,
-        // using a constraint object that takes care of constrained degrees
-        // (which here are only boundary nodes, as the individual levels have
-        // no hanging node constraints). Note that the
-        // <code>boundary_constraints</code> object makes sure that the level
-        // matrices contains no contributions from degrees of freedom at the
-        // interface between cells of different refinement level.
-        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 they 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.
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          for (unsigned int j=0; j<dofs_per_cell; ++j)
-            if ( !(mg_constrained_dofs.get_refinement_edge_indices(cell->level()).
-                   is_element(local_dof_indices[i])==true &&
-                   mg_constrained_dofs.get_refinement_edge_indices(cell->level()).
-                   is_element(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()]);
-      }
+    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);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+    MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
+    assembler.initialize(mg_constrained_dofs);
+    assembler.initialize(mg_matrices);
+    assembler.initialize_interfaces(mg_interface_in, mg_interface_out);
+
+    LaplaceMatrix<dim> matrix_integrator;
+    MeshWorker::integration_loop<dim, dim> (
+        dof_handler.begin_mg(), dof_handler.end_mg(),
+        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<dof_handler.n_dofs(level); ++i)
+        if(mg_matrices[level].diag_element(i)==0)
+          mg_matrices[level].set(i,i,1.);
+    }
   }
 
 
@@ -671,8 +530,8 @@ namespace Step16
     // initialize both up and down versions of the operator with the matrices
     // we already built:
     mg::Matrix<Vector<double> > mg_matrix(mg_matrices);
-    mg::Matrix<Vector<double> > mg_interface_up(mg_interface_matrices);
-    mg::Matrix<Vector<double> > mg_interface_down(mg_interface_matrices);
+    mg::Matrix<Vector<double> > mg_interface_up(mg_interface_in);
+    mg::Matrix<Vector<double> > mg_interface_down(mg_interface_out);
 
     // Now, we are ready to set up the V-cycle operator and the multilevel
     // preconditioner.
@@ -775,6 +634,7 @@ namespace Step16
           }
         else
           refine_grid ();
+//            triangulation.refine_global (1);
 
 
         std::cout << "   Number of active cells:       "
@@ -800,7 +660,6 @@ int main ()
 {
   try
     {
-      using namespace dealii;
       using namespace Step16;
 
       deallog.depth_console (0);
index 6135a031f4fb65745327ddc17fd3886c85413abb..a567b5775cb4cee1c8e022ad2bb65da3463c942b 100644 (file)
@@ -1140,11 +1140,11 @@ namespace MeshWorker
                     {
                       if (mg_constrained_dofs->set_boundary_values())
                         {
-                          if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
-                               !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
+                          if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
+                               !mg_constrained_dofs->is_boundary_index(level2, kglobal))
                               ||
-                              (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
-                               mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
+                              (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
+                               mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
                                jglobal == kglobal))
                             global.add(jglobal, kglobal, local(j,k));
                         }
@@ -1205,11 +1205,11 @@ namespace MeshWorker
                     {
                       if (mg_constrained_dofs->set_boundary_values())
                         {
-                          if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
-                               !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
+                          if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
+                               !mg_constrained_dofs->is_boundary_index(level2, kglobal))
                               ||
-                              (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
-                               mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
+                              (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
+                               mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
                                jglobal == kglobal))
                             global.add(jglobal, kglobal, local(k,j));
                         }
index 297e2421a68bedfe2f0870a25d2f52977ceba9c2..aa304eeda8cbdf4785c05c9e11f1127f25565d4f 100644 (file)
@@ -995,11 +995,11 @@ namespace MeshWorker
             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
               {
-               if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+               if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                    !mg_constrained_dofs->is_boundary_index(level, i2[k]))
                    ||
-                   (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                    mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                   (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                    mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
                     i1[j] == i2[k]))
                  G.add(i1[j], i2[k], M(j,k));
               }
@@ -1025,11 +1025,11 @@ namespace MeshWorker
             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
               {
-               if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+               if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                    !mg_constrained_dofs->is_boundary_index(level, i2[k]))
                    ||
-                   (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                    mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                   (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+                    mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
                     i1[j] == i2[k]))
                  G.add(i1[j], i2[k], M(k,j));
               }

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.