]> https://gitweb.dealii.org/ - dealii.git/commitdiff
clean code in Assembler::MGMatrixSimple
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 23 Aug 2013 20:20:53 +0000 (20:20 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 23 Aug 2013 20:20:53 +0000 (20:20 +0000)
add mg::SparseMatrixCollection

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

deal.II/include/deal.II/meshworker/simple.h
deal.II/include/deal.II/multigrid/sparse_matrix_collection.h [new file with mode: 0644]

index e7a9c73544054ffdd205bc279ff46961c2cfa948..e43257df1566aa2f821c7c5ff30dd8cb1c58b987 100644 (file)
@@ -256,29 +256,24 @@ namespace MeshWorker
     {
     public:
       /**
-       * Constructor, initializing
-       * the #threshold, which
-       * limits how small numbers
-       * may be to be entered into
-       * the matrix.
+       * Constructor, initializing the #threshold, which limits how
+       * small numbers may be to be entered into the matrix.
        */
       MGMatrixSimple(double threshold = 1.e-12);
 
       /**
-       * Store the result matrix
-       * for later assembling.
+       * Store the result matrix for later assembling.
        */
       void initialize(MGLevelObject<MATRIX> &m);
 
       /**
-       * Initialize the multilevel
-       * constraints.
+       * Initialize the multilevel constraints.
        */
       void initialize(const MGConstrainedDoFs &mg_constrained_dofs);
 
       /**
-       * @deprecated This function is of no effect. Only the block info
-       * structure in DoFInfo is being used.
+       * @deprecated This function is of no effect. Only the block
+       * info structure in DoFInfo is being used.
        *
        * Store information on the local block structure. If the
        * assembler is inititialized with this function,
@@ -920,8 +915,7 @@ namespace MeshWorker
           for (unsigned int j=0; j<i1.size(); ++j)
             for (unsigned int k=0; k<i2.size(); ++k)
               if (std::fabs(M(k,j)) >= threshold)
-                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(level, i2[k]))
                   G.add(i1[j], i2[k], M(k,j));
         }
     }
@@ -950,8 +944,7 @@ namespace MeshWorker
           for (unsigned int j=0; j<i1.size(); ++j)
             for (unsigned int k=0; k<i2.size(); ++k)
               if (std::fabs(M(j,k)) >= threshold)
-                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(level, i2[k]))
                   G.add(i1[j], i2[k], M(j,k));
         }
     }
@@ -967,49 +960,41 @@ namespace MeshWorker
     {
       AssertDimension(M.m(), i1.size());
       AssertDimension(M.n(), i2.size());
-
-      if (mg_constrained_dofs == 0)
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                G.add(i1[j], i2[k], M(j,k));
-        }
-      else
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                // Enter values into matrix only if j corresponds to a
-                // degree of freedom on the refinemenent edge, k does
-                // not, and both are not on the boundary. This is part
-                // the difference between the complete matrix with no
-                // boundary condition at the refinement edge and and
-                // the matrix assembled above by assemble().
-
-                // Thus the logic is: enter the row if it is
-                // constrained by hanging node constraints (actually,
-                // the whole refinement edge), but not if it is
-                // constrained by a boundary constraint.
-                if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                  {
-                    if (mg_constrained_dofs->set_boundary_values())
-                      {
-                        if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                             !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
-                            ||
-                            (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                             mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                             i1[j] == i2[k]))
-                          G.add(i1[j], i2[k], M(j,k));
-                      }
-                    else
-                      G.add(i1[j], i2[k], M(j,k));
-                  }
-        }
+      Assert(mg_constrained_dofs != 0, ExcInternalError());
+      
+      for (unsigned int j=0; j<i1.size(); ++j)
+       for (unsigned int k=0; k<i2.size(); ++k)
+         if (std::fabs(M(j,k)) >= threshold)
+                                            // Enter values into matrix only if j corresponds to a
+                                            // degree of freedom on the refinemenent edge, k does
+                                            // not, and both are not on the boundary. This is part
+                                            // the difference between the complete matrix with no
+                                            // boundary condition at the refinement edge and and
+                                            // the matrix assembled above by assemble().
+           
+                                            // Thus the logic is: enter the row if it is
+                                            // constrained by hanging node constraints (actually,
+                                            // the whole refinement edge), but not if it is
+                                            // constrained by a boundary constraint.
+           if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+               !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+             {
+               if (mg_constrained_dofs->set_boundary_values())
+                 {
+                   if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                       ||
+                       (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                        i1[j] == i2[k]))
+                     G.add(i1[j], i2[k], M(j,k));
+                 }
+               else
+                 G.add(i1[j], i2[k], M(j,k));
+             }
     }
-
+    
+    
     template <class MATRIX>
     inline void
     MGMatrixSimple<MATRIX>::assemble_out(
@@ -1021,38 +1006,29 @@ namespace MeshWorker
     {
       AssertDimension(M.n(), i1.size());
       AssertDimension(M.m(), i2.size());
-
-      if (mg_constrained_dofs == 0)
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(k,j)) >= threshold)
-                G.add(i1[j], i2[k], M(k,j));
-        }
-      else
-        {
-          for (unsigned int j=0; j<i1.size(); ++j)
-            for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(k,j)) >= threshold)
-                if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                  {
-                    if (mg_constrained_dofs->set_boundary_values())
-                      {
-                        if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                             !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
-                            ||
-                            (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                             mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                             i1[j] == i2[k]))
-                          G.add(i1[j], i2[k], M(k,j));
-                      }
-                    else
-                      G.add(i1[j], i2[k], M(k,j));
-                  }
-        }
+      Assert(mg_constrained_dofs != 0, ExcInternalError());
+      
+      for (unsigned int j=0; j<i1.size(); ++j)
+       for (unsigned int k=0; k<i2.size(); ++k)
+         if (std::fabs(M(k,j)) >= threshold)
+           if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+               !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+             {
+               if (mg_constrained_dofs->set_boundary_values())
+                 {
+                   if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                       ||
+                       (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                        i1[j] == i2[k]))
+                     G.add(i1[j], i2[k], M(k,j));
+                 }
+               else
+                 G.add(i1[j], i2[k], M(k,j));
+             }
     }
-
+  
 
     template <class MATRIX>
     template <class DOFINFO>
diff --git a/deal.II/include/deal.II/multigrid/sparse_matrix_collection.h b/deal.II/include/deal.II/multigrid/sparse_matrix_collection.h
new file mode 100644 (file)
index 0000000..a7dd038
--- /dev/null
@@ -0,0 +1,124 @@
+// ---------------------------------------------------------------------
+// $Id: mg_matrix.h 30036 2013-07-18 16:55:32Z maier $
+//
+// Copyright (C) 2003 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__mg_sparse_matrix_collection_h
+#define __deal2__mg_sparse_matrix_collection_h
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/pointer_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/multigrid/mg_base.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/base/mg_level_object.h>
+#include <deal.II/base/std_cxx1x/shared_ptr.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace mg
+{
+/**
+ * Handler and storage for all five SparseMatrix object involved in
+ * using multigrid with local refinement.
+ *
+ * @author Baerbel Janssen, Guido Kanschat
+ * @date 2013
+ */
+  template <typename number>
+  class SparseMatrixCollection : public Subscriptor
+  {
+    public:
+      void resize(const unsigned int minlevel, const unsigned  int maxlevel);
+      
+      template <class DH>
+      void reinit(const DH& dof_handler);
+
+      void set_zero();
+      
+      MGLevelObject<SparsityPattern> sparsity;
+      MGLevelObject<SparsityPattern> sparsity_edge;
+      
+      MGLevelObject<SparseMatrix<number> > matrix;
+      MGLevelObject<SparseMatrix<number> > matrix_down;
+      MGLevelObject<SparseMatrix<number> > matrix_up;
+      MGLevelObject<SparseMatrix<number> > matrix_in;
+      MGLevelObject<SparseMatrix<number> > matrix_out;
+  };
+
+
+  template <typename number>
+  void
+  SparseMatrixCollection<number>::resize(const unsigned int minlevel, const unsigned  int maxlevel)
+  {
+    matrix.resize(minlevel, maxlevel);
+    matrix.clear();
+    matrix_up.resize(minlevel+1, maxlevel);
+    matrix_up.clear();
+    matrix_down.resize(minlevel+1, maxlevel);
+    matrix_down.clear();
+    matrix_in.resize(minlevel, maxlevel);
+    matrix_in.clear();
+    matrix_out.resize(minlevel, maxlevel);
+    matrix_out.clear();
+    sparsity.resize(minlevel, maxlevel);
+    sparsity_edge.resize(minlevel, maxlevel);
+  }
+
+
+  template <typename number>
+  template <class DH>
+  void
+  SparseMatrixCollection<number>::reinit(const DH& dof_handler)
+  {
+    AssertIndexRange(sparsity.max_level(), dof_handler.get_tria().n_levels());
+    
+    for (unsigned int level=sparsity.min_level();
+        level<=sparsity.max_level();++level)
+      {
+       CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level));      
+       MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level);
+       sparsity[level].copy_from(c_sparsity);
+       matrix[level].reinit(sparsity[level]);
+       matrix_in[level].reinit(sparsity[level]);
+       matrix_out[level].reinit(sparsity[level]);
+       if (level>0)
+         {
+           CompressedSparsityPattern ci_sparsity;
+           ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level));
+           MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level);
+           sparsity_edge[level].copy_from(ci_sparsity);
+           matrix_up[level].reinit(sparsity_edge[level]);
+           matrix_down[level].reinit(sparsity_edge[level]);
+         }
+      }
+  }
+
+  template <typename number>
+  void
+  SparseMatrixCollection<number>::set_zero()
+  {
+    matrix = 0.;
+    matrix_in = 0.;
+    matrix_out = 0.;
+    matrix_up = 0.;
+    matrix_down = 0.;
+  }
+  
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.