]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
local blocks for Assembler::MGMatrixSimple
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 14 Mar 2012 17:30:50 +0000 (17:30 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 14 Mar 2012 17:30:50 +0000 (17:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@25280 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/meshworker/dof_info.h
deal.II/include/deal.II/meshworker/dof_info.templates.h
deal.II/include/deal.II/meshworker/simple.h

index 40addebb6c6ab8061e34f15d5b72edc7f5879ac9..f09b6ab2d30bf190893afb9b07087d3cf207c28e 100644 (file)
@@ -187,7 +187,9 @@ namespace MeshWorker
                                        * DoFInfoBox.
                                        */
       DoFInfo ();
-
+      
+                                      /// Set up local block indices
+      void set_block_indices ();
                                       /// Fill index vector with active indices
       void get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator& c);
 
index 76faf577b6d6227c616c25fb7c975f194cc72ad2..3ca29e4bc2f323c2d2ead362015d30a68f781cac 100644 (file)
@@ -33,7 +33,24 @@ namespace MeshWorker
   DoFInfo<dim,spacedim,number>::DoFInfo()
   {}
 
+  
+  template <int dim, int spacedim, typename number>
+  void
+  DoFInfo<dim,spacedim,number>::set_block_indices()
+  {
+    for (unsigned int i=0;i<indices.size();++i)
+      {
+       const std::pair<unsigned int, unsigned int>
+         bi = block_info->local().global_to_local(this->block_info->renumber(i));
+       indices_by_block[bi.first][bi.second] = indices_org[i];
+      }
+                                    // Remove this after
+                                    // changing block codes
+    for (unsigned int i=0;i<indices.size();++i)
+      indices[this->block_info->renumber(i)] = indices_org[i];
+  }
 
+  
   template <int dim, int spacedim, typename number>
   void
   DoFInfo<dim,spacedim,number>::get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator& c)
@@ -48,16 +65,7 @@ namespace MeshWorker
          {
            indices_org.resize(c->get_fe().dofs_per_cell);
            c->get_dof_indices(indices_org);
-           for (unsigned int i=0;i<indices.size();++i)
-             {
-               const std::pair<unsigned int, unsigned int>
-                 bi = block_info->local().global_to_local(this->block_info->renumber(i));
-               indices_by_block[bi.first][bi.second] = indices_org[i];
-             }
-                                            // Remove this after
-                                            // changing block codes
-           for (unsigned int i=0;i<indices.size();++i)
-             indices[this->block_info->renumber(i)] = indices_org[i];
+           set_block_indices();
          }
       }
     else
@@ -79,8 +87,7 @@ namespace MeshWorker
       {
        indices_org.resize(c->get_fe().dofs_per_cell);
        c->get_mg_dof_indices(indices_org);
-       for (unsigned int i=0;i<indices.size();++i)
-         indices[this->block_info->renumber(i)] = indices_org[i];
+       set_block_indices();
       }
     level_cell = true;
   }
index f3c5e6b4c2212d9e5930cef2e3096225511c3cb7..6aec7dcd698f33084fd001155a291acf5736635f 100644 (file)
@@ -301,6 +301,38 @@ namespace MeshWorker
                                          */
         void initialize(const MGConstrainedDoFs& mg_constrained_dofs);
 
+                                        /**
+                                         * Store information on the
+                                         * local block structure. If
+                                         * the assembler is
+                                         * inititialized with this
+                                         * function,
+                                         * initialize_info() will
+                                         * generate one local matrix
+                                         * for each block row and
+                                         * column, which will be
+                                         * numbered
+                                         * lexicographically, row by
+                                         * row.
+                                         *
+                                         * In spite of using local
+                                         * block structure, all
+                                         * blocks will be enteres
+                                         * into the same global
+                                         * matrix, disregarding any
+                                         * global block structure.
+                                         *
+                                         * @note The argument of this
+                                         * function will be copied
+                                         * into the member object
+                                         * #local_indices. Thus,
+                                         * every subsequent change in
+                                         * the block structure must
+                                         * be initialzied or will not
+                                         * be used by the assembler.
+                                         */    
+       void initialize_local_blocks(const BlockIndices& local_indices);
+       
                                         /**
                                          * Initialize the matrices
                                          * #flux_up and #flux_down
@@ -461,6 +493,16 @@ namespace MeshWorker
                                          */
         SmartPointer<const MGConstrainedDoFs,MGMatrixSimple<MATRIX> > mg_constrained_dofs;
 
+                                        /**
+                                         * The object containing the
+                                         * local block structure. Set
+                                         * by
+                                         * initialize_local_blocks()
+                                         * and used by assembling
+                                         * functions.
+                                         */
+       BlockIndices local_indices;
+
                                         /**
                                          * The smallest positive
                                          * number that will be
@@ -768,6 +810,14 @@ namespace MeshWorker
       mg_constrained_dofs = &c;
     }
 
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::initialize_local_blocks(const BlockIndices& b)
+    {
+      local_indices = b;
+    }
+
+    
     template <class MATRIX>
     inline void
     MGMatrixSimple<MATRIX>::initialize_fluxes(
@@ -793,7 +843,26 @@ namespace MeshWorker
     inline void
     MGMatrixSimple<MATRIX>::initialize_info(DOFINFO& info, bool face) const
     {
-      info.initialize_matrices(1, face);
+      const unsigned int n = local_indices.size();
+      
+      if (n == 0)
+       info.initialize_matrices(1, face);
+      else
+       {
+         info.initialize_matrices(n*n, face);
+         unsigned int k=0;
+         for (unsigned int i=0;i<n;++i)
+           for (unsigned int j=0;j<n;++j,++k)
+             {
+               info.matrix(k,false).row = i;
+               info.matrix(k,false).column = j;
+               if (face)
+                 {
+                   info.matrix(k,true).row = i;
+                   info.matrix(k,true).column = j;
+                 }
+             }
+       }
     }
 
 
@@ -1025,16 +1094,39 @@ namespace MeshWorker
     MGMatrixSimple<MATRIX>::assemble(const DOFINFO& info)
     {
       const unsigned int level = info.cell->level();
-      assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
 
-      if(mg_constrained_dofs != 0)
+      if (local_indices.size() == 0)
        {
-         assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
-         assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
+         assemble((*matrix)[level], info.matrix(0,false).matrix,
+                  info.indices, info.indices, level);
+         if(mg_constrained_dofs != 0)
+           {
+             assemble_in((*interface_in)[level], info.matrix(0,false).matrix,
+                         info.indices, info.indices, level);
+             assemble_out((*interface_out)[level],info.matrix(0,false).matrix,
+                          info.indices, info.indices, level);
+           }
        }
+      else
+       for (unsigned int k=0;k<info.n_matrices();++k)
+         {
+           const unsigned int row = info.matrix(k,false).row;
+           const unsigned int column = info.matrix(k,false).column;
+           
+           assemble((*matrix)[level], info.matrix(k,false).matrix,
+                    info.indices_by_block[row], info.indices_by_block[column], level);
+           
+           if(mg_constrained_dofs != 0)
+             {
+               assemble_in((*interface_in)[level], info.matrix(k,false).matrix,
+                           info.indices, info.indices, level);
+               assemble_out((*interface_out)[level],info.matrix(k,false).matrix,
+                            info.indices_by_block[row], info.indices_by_block[column], level); 
+             }
+         }
     }
-
-
+    
+    
     template <class MATRIX>
     template <class DOFINFO>
     inline void
@@ -1044,39 +1136,79 @@ namespace MeshWorker
       const unsigned int level1 = info1.cell->level();
       const unsigned int level2 = info2.cell->level();
 
-      if (level1 == level2)
+      
+      if (local_indices.size() == 0)
        {
-          if(mg_constrained_dofs == 0)
+         if (level1 == level2)
            {
-             assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
-             assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
-             assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
-             assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
+             if(mg_constrained_dofs == 0)
+               {
+                 assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+                 assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
+                 assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
+                 assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
+               }
+             else
+               {
+                 assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
+                 assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
+                 assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
+                 assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+               }
            }
-          else
+         else
            {
-             assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
-             assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
-             assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
-             assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+             Assert(level1 > level2, ExcInternalError());
+                                              // Do not add info2.M1,
+                                              // which is done by
+                                              // the coarser cell
+             assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+             if(level1>0)
+               {
+                 assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+                 assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+               }
            }
        }
       else
+       for (unsigned int k=0;k<info1.n_matrices();++k)
        {
-         Assert(level1 > level2, ExcInternalError());
-                                          // Do not add info2.M1,
-                                          // which is done by
-                                          // the coarser cell
-          assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
-                                          //assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices);
-                                          //assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
-          if(level1>0)
+         const unsigned int row = info1.matrix(k,false).row;
+         const unsigned int column = info1.matrix(k,false).column;
+         
+         if (level1 == level2)
+           {
+             if(mg_constrained_dofs == 0)
+               {
+                 assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
+                 assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]);
+                 assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]);
+                 assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]);
+               }
+             else
+               {
+                 assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
+                 assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
+                 assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
+                 assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
+               }
+           }
+         else
            {
-             assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
-             assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+             Assert(level1 > level2, ExcInternalError());
+                                              // Do not add info2.M1,
+                                              // which is done by
+                                              // the coarser cell
+             assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
+             if(level1>0)
+               {
+                 assemble_up((*flux_up)[level1],info1.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
+                 assemble_down((*flux_down)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
+               }
            }
        }
     }
+    
 //----------------------------------------------------------------------//
 
     template <class MATRIX, class VECTOR>

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.