]> https://gitweb.dealii.org/ - dealii.git/commitdiff
prepared MGTransferSelect for adaptive refinement, cleaned up mg_dof_tools.cc
authorBaerbel Jannsen <baerbel.janssen@gmail.com>
Fri, 5 Feb 2010 13:56:31 +0000 (13:56 +0000)
committerBaerbel Jannsen <baerbel.janssen@gmail.com>
Fri, 5 Feb 2010 13:56:31 +0000 (13:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@20506 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/multigrid/mg_dof_tools.cc
deal.II/deal.II/source/multigrid/mg_transfer_component.cc
deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc

index 46c2d251349640038afe28281978009419ee8529..aa11f8b5d652d0a15cd6b46fbc5c984a83a37653 100644 (file)
@@ -30,7 +30,6 @@
 #include <dofs/dof_tools.h>
 #include <fe/fe.h>
 
-#include <iostream>
 #include <vector>
 #include <algorithm>
 #include <numeric>
@@ -1072,7 +1071,10 @@ MGTools::count_dofs_per_block (
       std::vector<std::vector<bool> >
        dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
       std::vector<std::vector<bool> >
-       block_select (n_blocks, std::vector<bool>(n_blocks, false));
+       block_select;
+     block_select.resize(n_blocks);
+     for (int iii=0;iii<n_blocks;++iii)
+       block_select[iii].resize(n_blocks, false);
       Threads::TaskGroup<> tasks;
       for (unsigned int i=0; i<n_blocks; ++i)
        {
@@ -1301,6 +1303,68 @@ MGTools::make_boundary_list(
 
 
 
+template <int dim, int spacedim>
+void
+MGTools::
+extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
+                             std::vector<std::vector<bool> >  &interface_dofs)
+{
+  Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_levels(),
+         ExcDimensionMismatch (interface_dofs.size(),
+                               mg_dof_handler.get_tria().n_levels()));
+
+  for (unsigned int l=0; l<mg_dof_handler.get_tria().n_levels(); ++l)
+    {
+      Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
+             ExcDimensionMismatch (interface_dofs[l].size(),
+                                   mg_dof_handler.n_dofs(l)));
+
+      std::fill (interface_dofs[l].begin(),
+                interface_dofs[l].end(),
+                false);
+    }
+
+  const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  const unsigned int   dofs_per_face   = fe.dofs_per_face;
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  std::vector<bool> cell_dofs(dofs_per_cell, false);
+
+  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+                                           endc = mg_dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      std::fill (cell_dofs.begin(), cell_dofs.end(), false);
+
+      for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+       {
+         const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
+         if (!face->at_boundary())
+           {
+                                              //interior face
+             const typename MGDoFHandler<dim>::cell_iterator
+               neighbor = cell->neighbor(face_nr);
+
+             if (neighbor->level() < cell->level())
+               {
+                 for (unsigned int j=0; j<dofs_per_face; ++j)
+                   cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+
+               }
+           }
+       }
+
+      const unsigned int level = cell->level();
+      cell->get_mg_dof_indices (local_dof_indices);
+
+      for(unsigned int i=0; i<dofs_per_cell; ++i)
+         if (cell_dofs[i])
+           interface_dofs[level][local_dof_indices[i]] = true;
+    }
+}
 
 
 template <int dim, int spacedim>
@@ -1380,7 +1444,7 @@ extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
        for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
          if(cell->at_boundary(face_nr))
            for(unsigned int j=0; j<dofs_per_face; ++j)
-             if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true)
+             if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) //is this necessary?
                boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 
 
@@ -1476,6 +1540,11 @@ MGTools::
 extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
                              std::vector<std::vector<bool> >  &interface_dofs,
                              std::vector<std::vector<bool> >  &boundary_interface_dofs);
+template
+void
+MGTools::
+extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
+                             std::vector<std::vector<bool> >  &interface_dofs);
 #endif
 
 DEAL_II_NAMESPACE_CLOSE
index 7d2ae176eb4ada772521cc55ef3b8cf4977473f7..26d2891c864eda7709396bc27d61cf6495b4dad3 100644 (file)
 //---------------------------------------------------------------------------
 
 #include <base/logstream.h>
+#include <base/function.h>
 
 #include <lac/vector.h>
 #include <lac/block_vector.h>
+#include <lac/block_indices.h>
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparse_matrix.h>
 #include <grid/tria.h>
@@ -29,7 +31,7 @@
 
 #include <algorithm>
 #include <numeric>
-
+#include <iostream>
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -123,8 +125,7 @@ namespace
          new_dofs(mg_dof.get_tria().n_levels(),
                   std::vector<unsigned int>(target_component.size()));
        std::swap(ndofs, new_dofs);
-       MGTools::count_dofs_per_component (mg_dof, ndofs,
-                                          true, target_component);
+       MGTools::count_dofs_per_block (mg_dof, ndofs, target_component);
       }
 
     for (unsigned int level=v.get_minlevel();
@@ -184,25 +185,26 @@ namespace
 
                                     // Compute the number of blocks needed
 #ifdef DEBUG
-    const unsigned int n_selected
-      = std::accumulate(selected.begin(),
-                       selected.end(),
-                       0U);
-    Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1));
+//    const unsigned int n_selected
+//      = std::accumulate(selected.begin(),
+//                     selected.end(),
+//                     0U);
+//    Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1));
 #endif
 
     unsigned int selected_block = 0;
-    while (!selected[selected_block])
-      ++selected_block;
-
+    for (unsigned int i=0; i<target_component.size(); ++i)
+      if(selected[i])
+        selected_block = target_component[i];
+    
     if (ndofs.size() == 0)
       {
        std::vector<std::vector<unsigned int> >
          new_dofs(mg_dof.get_tria().n_levels(),
                   std::vector<unsigned int>(target_component.size()));
        std::swap(ndofs, new_dofs);
-       MGTools::count_dofs_per_component (mg_dof, ndofs,
-                                          true, target_component);
+       MGTools::count_dofs_per_block (mg_dof, ndofs,
+                                          target_component);
       }
 
     for (unsigned int level=v.get_minlevel();
@@ -220,18 +222,8 @@ void
 MGTransferSelect<number>::do_copy_to_mg (
   const MGDoFHandler<dim,spacedim>&        mg_dof_handler,
   MGLevelObject<Vector<number> >& dst,
-  const InVector&                 src,
-  const unsigned int              offset) const
+  const InVector&                 src) const
 {
-  const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
-
-  const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
-                                  // set the elements of the vectors
-                                  // on all levels to zero
-  unsigned int minlevel = dst.get_minlevel();
-  unsigned int maxlevel = dst.get_maxlevel();
-
   dst=0;
 
   Assert(sizes.size()==mg_dof_handler.get_tria().n_levels(),
@@ -240,21 +232,6 @@ MGTransferSelect<number>::do_copy_to_mg (
   reinit_vector_by_components(mg_dof_handler, dst, mg_selected,
                              mg_target_component, sizes);
 
-  std::vector<unsigned int> global_dof_indices (dofs_per_cell);
-  std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
-
-                                  // Build a vector of the selected
-                                  // indices, since traversing all
-                                  // indices on each cell is too
-                                  // slow.
-  std::vector<unsigned int> selected_indices;
-  selected_indices.reserve(dofs_per_cell);
-  for (unsigned int i=0; i<dofs_per_cell; ++i)
-    if (mg_selected[mg_target_component[fe.system_to_component_index(i).first]])
-      selected_indices.push_back(i);
-  unsigned int selected_component
-    = mg_target_component[fe.system_to_component_index(selected_indices[0]).first];
-
                                   // traverse the grid top-down
                                   // (i.e. starting with the most
                                   // refined grid). this way, we can
@@ -265,50 +242,24 @@ MGTransferSelect<number>::do_copy_to_mg (
                                   // the respective vector on the
                                   // next finer level, which we then
                                   // already have built.
-  for (int level=maxlevel; level>=static_cast<signed int>(minlevel); --level)
-    {
 
-      typename MGDoFHandler<dim,spacedim>::active_cell_iterator
-       level_cell = mg_dof_handler.begin_active(level);
-      const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
-       level_end  = mg_dof_handler.end_active(level);
-                                      // Compute coarse level right hand side
-                                      // by restricting from fine level.
-      for (; level_cell!=level_end; ++level_cell)
-       {
-                                          // get the dof numbers of
-                                          // this cell for the global
-                                          // and the level-wise
-                                          // numbering
-         level_cell->get_dof_indices(global_dof_indices);
-         level_cell->get_mg_dof_indices (level_dof_indices);
-
-                                          // transfer the global
-                                          // defect in the vector
-                                          // into the level-wise one
-         const unsigned int level_start
-           = mg_component_start[level][selected_component];
-         const typename std::vector<unsigned int>::const_iterator
-           end = selected_indices.end();
-
-         for (typename std::vector<unsigned int>::const_iterator
-                i=selected_indices.begin();
-              i != end; ++i)
-           {
-             dst[level](level_dof_indices[*i] - level_start)
-               = src(global_dof_indices[*i] - offset);
-           }
-       }
+  bool first = true;
+  for (unsigned int level=mg_dof_handler.get_tria().n_levels(); level!=0;)
+    {
+      --level;
 
+      typedef std::map<unsigned int, unsigned int>::const_iterator IT;
+         for (IT i=copy_to_and_from_indices[level].begin();
+              i != copy_to_and_from_indices[level].end(); ++i)
+              dst[level](i->second) = src(i->first);
                                       // for that part of the level
                                       // which is further refined:
                                       // get the defect by
                                       // restriction of the defect on
                                       // one level higher
-      if (static_cast<unsigned int>(level) < maxlevel)
-       {
-         restrict_and_add (level+1, dst[level], dst[level+1]);
-       }
+          if(!first)
+            restrict_and_add (level+1, dst[level], dst[level+1]);
+          first = false;
     }
 }
 
@@ -375,7 +326,10 @@ void MGTransferComponentBase::build_matrices (
 
                                   // Compute the lengths of all blocks
   sizes.resize(n_levels);
-  MGTools::count_dofs_per_component(mg_dof, sizes, true, mg_target_component);
+  for(unsigned int l=0; l<n_levels; ++l)
+    sizes[l].resize(n_components);
+
+  MGTools::count_dofs_per_block(mg_dof, sizes, mg_target_component);
 
                                   // Fill some index vectors
                                   // for later use.
@@ -393,10 +347,9 @@ void MGTransferComponentBase::build_matrices (
     }
 
   component_start.resize(*std::max_element (target_component.begin(),
-                                           target_component.end()) + 1);
+                                       target_component.end()) + 1);
   DoFTools::
-    count_dofs_per_component (static_cast<const DoFHandler<dim,spacedim>&>(mg_dof),
-                              component_start, true, target_component);
+    count_dofs_per_block (mg_dof, component_start, target_component);
 
   unsigned int k=0;
   for (unsigned int i=0;i<component_start.size();++i)
@@ -504,7 +457,7 @@ void MGTransferComponentBase::build_matrices (
                          = fe.system_to_component_index(i).first;
                        const unsigned int jcomp
                          = fe.system_to_component_index(j).first;
-                       if ((icomp==jcomp) && mg_selected[mg_target_component[icomp]])
+                       if ((icomp==jcomp) && mg_selected[icomp])
                          prolongation_sparsities[level]->add(dof_indices_child[i],
                                                              dof_indices_parent[j]);
                      };
@@ -538,7 +491,7 @@ void MGTransferComponentBase::build_matrices (
                      {
                        const unsigned int icomp = fe.system_to_component_index(i).first;
                        const unsigned int jcomp = fe.system_to_component_index(j).first;
-                       if ((icomp==jcomp) && mg_selected[mg_target_component[icomp]])
+                       if ((icomp==jcomp) && mg_selected[icomp])
                          prolongation_matrices[level]->set(dof_indices_child[i],
                                                            dof_indices_parent[j],
                                                            prolongation(i,j));
@@ -546,6 +499,56 @@ void MGTransferComponentBase::build_matrices (
              }
          }
     }
+  // impose boundary conditions
+  // but only in the column of
+  // the prolongation matrix
+  //TODO: this way is not very efficient
+       
+  if(boundary_indices.size() != 0)
+  {
+    std::vector<std::vector<unsigned int> > 
+      dofs_per_component(mg_dof.get_tria().n_levels(),
+          std::vector<unsigned int>(n_components));
+
+    MGTools::count_dofs_per_block (mg_dof, dofs_per_component, mg_target_component);
+    for (unsigned int level=0; level<n_levels-1; ++level)
+    {
+      if (boundary_indices[level].size() == 0)
+        continue;
+
+      for (unsigned int iblock=0; iblock<n_components; ++iblock)
+        for (unsigned int jblock=0; jblock<n_components; ++jblock)
+          if (iblock==jblock)
+          {
+            const unsigned int n_dofs = prolongation_matrices[level]->block(iblock,jblock).m();
+            for (unsigned int i=0; i<n_dofs; ++i)
+            {
+              SparseMatrix<double>::iterator anfang = prolongation_matrices[level]->block(iblock,jblock).begin(i),
+                ende = prolongation_matrices[level]->block(iblock,jblock).end(i);
+              for(; anfang != ende; ++anfang)
+              {
+                const unsigned int column_number = anfang->column();
+
+                //convert global indices into local ones
+                const BlockIndices block_indices_coarse (dofs_per_component[level]);
+                const unsigned int global_j = block_indices_coarse.local_to_global(iblock, column_number);
+
+                std::set<unsigned int>::const_iterator found_dof = 
+                  boundary_indices[level].find(global_j);
+
+                const bool is_boundary_index =
+                  (found_dof != boundary_indices[level].end());
+
+                if(is_boundary_index)
+                {
+                  prolongation_matrices[level]->block(iblock,jblock)
+                    .set(i,column_number,0);
+                }
+              }
+            }
+          }
+    }
+  }
 }
 
 
@@ -557,20 +560,27 @@ void MGTransferSelect<number>::build_matrices (
   unsigned int select,
   unsigned int mg_select,
   const std::vector<unsigned int>& t_component,
-  const std::vector<unsigned int>& mg_t_component)
+  const std::vector<unsigned int>& mg_t_component,
+  const std::vector<std::set<unsigned int> >& bdry_indices)
 {
   const FiniteElement<dim>& fe = mg_dof.get_fe();
   unsigned int ncomp = mg_dof.get_fe().n_components();
 
   target_component = t_component;
   mg_target_component = mg_t_component;
+  boundary_indices = bdry_indices;
 
   selected_component = select;
   mg_selected_component = mg_select;
   selected.resize(ncomp, false);
-  selected[select] = true;
+  for(unsigned int c=0; c<ncomp; ++c)
+    if(t_component[c] == selected_component)
+      selected[c] = true;
+
   mg_selected.resize(ncomp, false);
-  mg_selected[mg_select] = true;
+  for(unsigned int c=0; c<ncomp; ++c)
+    if(mg_t_component[c] == mg_selected_component)
+      mg_selected[c] = true;
                                   // If components are renumbered,
                                   // find the first original
                                   // component corresponding to the
@@ -592,12 +602,19 @@ void MGTransferSelect<number>::build_matrices (
          break;
        }
     }
+
   MGTransferComponentBase::build_matrices (dof, mg_dof);
 
+  interface_dofs.resize(mg_dof.get_tria().n_levels());
+  for(unsigned int l=0; l<mg_dof.get_tria().n_levels(); ++l)
+    interface_dofs[l].resize(mg_dof.n_dofs(l));
+  MGTools::extract_inner_interface_dofs(mg_dof, interface_dofs);
+
   std::vector<unsigned int> global_dof_indices (fe.dofs_per_cell);
   std::vector<unsigned int> level_dof_indices  (fe.dofs_per_cell);
   for (int level=dof.get_tria().n_levels()-1; level>=0; --level)
     {
+      copy_to_and_from_indices[level].clear();
       typename MGDoFHandler<dim,spacedim>::active_cell_iterator
        level_cell = mg_dof.begin_active(level);
       const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
@@ -618,14 +635,15 @@ void MGTransferSelect<number>::build_matrices (
          for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
            {
              const unsigned int component
-               = mg_target_component[fe.system_to_component_index(i).first];
-             if (mg_selected[component])
+               = fe.system_to_component_index(i).first;
+             if (mg_selected[component] && 
+                  !interface_dofs[level][level_dof_indices[i]])
                {
                  const unsigned int level_start
-                   = mg_component_start[level][component];
+                   = mg_component_start[level][mg_target_component[component]];
                  copy_to_and_from_indices[level].insert(
                    std::make_pair(global_dof_indices[i]
-                                  - component_start[selected_component],
+                                  - component_start[target_component[selected_component]],
                                   level_dof_indices[i]-level_start));
                }
            }
@@ -643,7 +661,8 @@ void MGTransferSelect<float>::build_matrices<deal_II_dimension>
  const MGDoFHandler<deal_II_dimension> &,
  unsigned int, unsigned int,
  const std::vector<unsigned int>&,
- const std::vector<unsigned int>&);
+ const std::vector<unsigned int>&,
+ const std::vector<std::set<unsigned int> >&);
 
 template
 void MGTransferSelect<double>::build_matrices<deal_II_dimension>
@@ -651,7 +670,8 @@ void MGTransferSelect<double>::build_matrices<deal_II_dimension>
  const MGDoFHandler<deal_II_dimension> &,
  unsigned int, unsigned int,
  const std::vector<unsigned int>&,
- const std::vector<unsigned int>&);
+ const std::vector<unsigned int>&,
+ const std::vector<std::set<unsigned int> >&);
 
 template void
 MGTransferSelect<float>::copy_to_mg (
index 0e676f06c7cb6bb5612c8e05c98be01f62270f1a..fa7bcc0bab35f73e6464cb056a3d6c8b4a2a8dd6 100644 (file)
@@ -123,7 +123,7 @@ MGTransferPrebuilt<VECTOR>::copy_to_mg (
       --level;
       VECTOR& dst_level = dst[level];
 
-      typedef std::vector<std::pair<unsigned int, unsigned int> >::const_iterator IT;
+      typedef std::map<unsigned int, unsigned int>::const_iterator IT;
       for (IT i= copy_indices[level].begin();
           i != copy_indices[level].end();++i)
        dst_level(i->second) = src(i->first);
@@ -145,7 +145,9 @@ MGTransferPrebuilt<VECTOR>::copy_to_mg (
 template <typename VECTOR>
 template <int dim, int spacedim>
 void MGTransferPrebuilt<VECTOR>::build_matrices (
-  const MGDoFHandler<dim,spacedim> &mg_dof)
+  const MGDoFHandler<dim,spacedim> &mg_dof,
+  const std::vector<std::set<unsigned int> >&boundary_indices
+  )
 {
   //start building the matrices here
   const unsigned int n_levels      = mg_dof.get_tria().n_levels();
@@ -274,42 +276,48 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
   // impose boundary conditions
   // but only in the column of
   // the prolongation matrix
-  std::vector<std::set<unsigned int> > boundary_indices (n_levels);
-  typename FunctionMap<dim>::type boundary;
-  ZeroFunction<dim>                    homogeneous_dirichlet_bc (1);
-  boundary[0] = &homogeneous_dirichlet_bc;
-  MGTools::make_boundary_list(mg_dof, boundary, boundary_indices);
-
-  std::vector<unsigned int> constrain_indices;
-  for (int level=n_levels-2; level>=0; --level)
+
+  if(boundary_indices.size() != 0)
+  {
+    std::vector<unsigned int> constrain_indices;
+    for (int level=n_levels-2; level>=0; --level)
+//    for (unsigned int level=0; level<n_levels-1; ++level)
     {
       if (boundary_indices[level].size() == 0)
-       continue;
+        continue;
 
-                               // need to delete all the columns in the
-                               // matrix that are on the boundary. to fix
-                               // this, create an array as long as there are
-                               // matrix columns, and find which columns we
-                               // need to filter away.
+      // need to delete all the columns in the
+      // matrix that are on the boundary. to fix
+      // this, create an array as long as there are
+      // matrix columns, and find which columns we
+      // need to filter away.
       constrain_indices.resize (0);
       constrain_indices.resize (prolongation_matrices[level]->n(), 0);
       std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
-       endd = boundary_indices[level].end();
+        endd = boundary_indices[level].end();
       for (; dof != endd; ++dof)
-       constrain_indices[*dof] = 1;
+        constrain_indices[*dof] = 1;
 
       const unsigned int n_dofs = prolongation_matrices[level]->m();
       for (unsigned int i=0; i<n_dofs; ++i)
-       {
-         SparseMatrix<double>::iterator start_row = prolongation_matrices[level]->begin(i),
-           end_row = prolongation_matrices[level]->end(i);
-         for(; start_row != end_row; ++start_row)
-           {
-             if (constrain_indices[start_row->column()] == 1)
-               start_row->value() = 0;
-           }
-       }
+      {
+        SparseMatrix<double>::iterator start_row = prolongation_matrices[level]->begin(i),
+          end_row = prolongation_matrices[level]->end(i);
+        for(; start_row != end_row; ++start_row)
+        {
+          std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
+            endd = boundary_indices[level].end();
+          for (; dof != endd; ++dof)
+          {
+            const unsigned int column_number = start_row->column();
+            const unsigned int dof_number = *dof;
+            if(dof_number == column_number)
+              prolongation_matrices[level]->set(i,dof_number,0);
+          }
+        }
+      }
     }
+  }
 
                                // to find the indices that describe the
                                // relation between global dofs and local
@@ -321,11 +329,13 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
                                // the std::vector<std::pair<uint,uint> > that
                                // only contains the active dofs on the
                                // levels.
-
-  find_dofs_on_refinement_edges (mg_dof);
+  interface_dofs.resize(mg_dof.get_tria().n_levels());
+  for(unsigned int l=0; l<mg_dof.get_tria().n_levels(); ++l)
+    interface_dofs[l].resize(mg_dof.n_dofs(l));
+  MGTools::extract_inner_interface_dofs (mg_dof, interface_dofs);
 
   copy_indices.resize(n_levels);
-  std::vector<unsigned int> temp_copy_indices;
+//  std::vector<unsigned int> temp_copy_indices;
   std::vector<unsigned int> global_dof_indices (dofs_per_cell);
   std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
   for (int level=mg_dof.get_tria().n_levels()-1; level>=0; --level)
@@ -336,8 +346,8 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
       const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
        level_end  = mg_dof.end_active(level);
 
-      temp_copy_indices.resize (0);
-      temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int);
+//      temp_copy_indices.resize (0);
+//      temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int);
 
                                       // Compute coarse level right hand side
                                       // by restricting from fine level.
@@ -352,8 +362,14 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
          level_cell->get_mg_dof_indices (level_dof_indices);
 
          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           if(!dofs_on_refinement_edge[level][level_dof_indices[i]])
-             temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+          {
+              if(!interface_dofs[level][level_dof_indices[i]])
+             copy_indices[level].insert(
+               std::make_pair(global_dof_indices[i], level_dof_indices[i]));
+//         temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+//         if(!interface_dofs[level][level_dof_indices[i]])
+//           temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
+          }
        }
 
                                // now all the active dofs got a valid entry,
@@ -362,17 +378,17 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
                                // copy_indices object. Then, insert the pairs
                                // of global index and level index into
                                // copy_indices.
-      const unsigned int n_active_dofs =
-       std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(),
-                      std::bind2nd(std::not_equal_to<unsigned int>(),
-                                   numbers::invalid_unsigned_int));
-      copy_indices[level].resize (n_active_dofs);
-      unsigned int counter = 0;
-      for (unsigned int i=0; i<temp_copy_indices.size(); ++i)
-       if (temp_copy_indices[i] != numbers::invalid_unsigned_int)
-         copy_indices[level][counter++] =
-           std::make_pair<unsigned int,unsigned int> (temp_copy_indices[i], i);
-      Assert (counter == n_active_dofs, ExcInternalError());
+//      const unsigned int n_active_dofs =
+//     std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(),
+//                    std::bind2nd(std::not_equal_to<unsigned int>(),
+//                                 numbers::invalid_unsigned_int));
+//      copy_indices[level].resize (n_active_dofs);
+//      unsigned int counter = 0;
+//      for (unsigned int i=0; i<temp_copy_indices.size(); ++i)
+//     if (temp_copy_indices[i] != numbers::invalid_unsigned_int)
+//       copy_indices[level][counter++] =
+//         std::make_pair<unsigned int,unsigned int> (temp_copy_indices[i], i);
+//      Assert (counter == n_active_dofs, ExcInternalError());
     }
 }
 
@@ -381,19 +397,23 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
 
 template
 void MGTransferPrebuilt<Vector<float> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
 
 template
 void MGTransferPrebuilt<Vector<double> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
 
 template
 void MGTransferPrebuilt<BlockVector<float> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
 
 template
 void MGTransferPrebuilt<BlockVector<double> >::build_matrices<deal_II_dimension>
-(const MGDoFHandler<deal_II_dimension> &mg_dof);
+(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const std::vector<std::set<unsigned int> >&boundary_indices);
 
 template void
 MGTransferPrebuilt<Vector<float> >::copy_to_mg (

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.