include/deal.II/lac/block_list.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: block_list.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__block_list_h
00013 #define __deal2__block_list_h
00014 
00015 
00016 #include <deal.II/base/subscriptor.h>
00017 #include <deal.II/base/template_constraints.h>
00018 #include <deal.II/lac/sparsity_pattern.h>
00019 #include <deal.II/fe/fe.h>
00020 
00021 #include <vector>
00022 
00023 DEAL_II_NAMESPACE_OPEN
00024 
00047 class BlockList :
00048   public Subscriptor
00049 {
00050   public:
00052     typedef std::vector<unsigned int> block_container;
00054     typedef block_container::const_iterator const_iterator;
00055 
00068     void create_sparsity_pattern(SparsityPattern& sparsity, unsigned int n) const;
00069 
00076     void add(unsigned int block, const std::vector<unsigned int>& indices);
00077 
00086     void add(unsigned int block,
00087              const std::vector<unsigned int>& indices,
00088              const std::vector<bool>& selected_indices,
00089              unsigned int offset = 0);
00090 
00095     void initialize(unsigned int n_blocks);
00096 
00112     template <typename ITERATOR>
00113     void initialize(unsigned int n_blocks,
00114                     const ITERATOR begin,
00115                     const typename identity<ITERATOR>::type end);
00116 
00135     template <typename ITERATOR>
00136     void initialize_mg(unsigned int n_blocks,
00137                        const ITERATOR begin,
00138                        const typename identity<ITERATOR>::type end);
00139 
00171     template <typename ITERATOR>
00172     void initialize(unsigned int n_blocks,
00173                     const ITERATOR begin,
00174                     const typename identity<ITERATOR>::type end,
00175                     const std::vector<bool>& selected_dofs,
00176                     unsigned int offset = 0);
00208     template <typename ITERATOR>
00209     void initialize_mg(unsigned int n_blocks,
00210                        const ITERATOR begin,
00211                        const typename identity<ITERATOR>::type end,
00212                        const std::vector<bool>& selected_dofs,
00213                        unsigned int offset = 0);
00214 
00226     template <int dim, typename ITERATOR>
00227     void initialize_vertex_patches_mg(unsigned int n_blocks,
00228                        const ITERATOR begin,
00229                        const typename identity<ITERATOR>::type end,
00230                        const std::vector<bool>& selected_dofs = std::vector<bool>(),
00231                        unsigned int offset = 0);
00232 
00240     template <int dim, typename ITERATOR>
00241     unsigned int count_vertex_patches(
00242       const ITERATOR begin,
00243       const typename identity<ITERATOR>::type end,
00244       bool same_level_only) const;
00245 
00251     template <int dim, typename ITERATOR>
00252     bool cell_generates_vertex_patch(const ITERATOR cell,
00253                                      bool same_level_only) const;
00254 
00258     unsigned int size() const;
00259 
00263     unsigned int block_size(unsigned int block) const;
00264 
00268     const_iterator begin(unsigned int block) const;
00272     const_iterator end(unsigned int block) const;
00280     unsigned int local_index(unsigned int block, unsigned int index) const;
00281 
00282   private:
00286     std::vector<block_container> index_sets;
00287 };
00288 
00289 
00290 inline
00291 void
00292 BlockList::create_sparsity_pattern(SparsityPattern& sparsity, unsigned int n) const
00293 {
00294   std::vector<unsigned int> sizes(size());
00295   for (unsigned int b=0;b<size();++b)
00296     sizes[b] = block_size(b);
00297 
00298   sparsity.reinit(size(), n, sizes);
00299   for (unsigned int b=0;b<size();++b)
00300     {
00301       for (const_iterator i = begin(b); i != end(b); ++i)
00302         sparsity.add(b,*i);
00303     }
00304   sparsity.compress();
00305 }
00306 
00307 
00308 inline
00309 void
00310 BlockList::add(const unsigned int block, const std::vector<unsigned int>& indices)
00311 {
00312   AssertIndexRange(block, index_sets.size());
00313 
00314   for (unsigned int i=0;i<indices.size();++i)
00315     {
00316       const unsigned int k = indices[i];
00317       if (k==numbers::invalid_unsigned_int)
00318         continue;
00319       if (std::find(index_sets[block].begin(), index_sets[block].end(), k)
00320           == index_sets[block].end())
00321         index_sets[block].push_back(k);
00322     }
00323 }
00324 
00325 
00326 inline
00327 void
00328 BlockList::add(
00329   const unsigned int block,
00330   const std::vector<unsigned int>& indices,
00331   const std::vector<bool>& selected,
00332   unsigned int offset)
00333 {
00334   AssertIndexRange(block, index_sets.size());
00335   AssertDimension(indices.size(), selected.size());
00336 
00337   for (unsigned int i=0;i<indices.size();++i)
00338     {
00339       const unsigned int k = indices[i];
00340       if (k==numbers::invalid_unsigned_int)
00341         continue;
00342       if (selected[i] && std::find(index_sets[block].begin(), index_sets[block].end(), k-offset)
00343           == index_sets[block].end())
00344         index_sets[block].push_back(k-offset);
00345     }
00346 }
00347 
00348 
00349 inline
00350 void
00351 BlockList::initialize(unsigned int n_blocks)
00352 {
00353   index_sets.resize(n_blocks);
00354 }
00355 
00356 
00357 template <typename ITERATOR>
00358 inline
00359 void
00360 BlockList::initialize(unsigned int n_blocks, const ITERATOR begin, const typename identity<ITERATOR>::type end)
00361 {
00362   index_sets.resize(n_blocks);
00363   std::vector<unsigned int> indices;
00364   unsigned int k = 0;
00365   for (ITERATOR cell = begin; cell != end; ++cell, ++k)
00366     {
00367       indices.resize(cell->get_fe().dofs_per_cell);
00368       cell->get_dof_indices(indices);
00369       add(k, indices);
00370     }
00371 }
00372 
00373 
00374 template <typename ITERATOR>
00375 inline
00376 void
00377 BlockList::initialize_mg(unsigned int n_blocks, const ITERATOR begin, const typename identity<ITERATOR>::type end)
00378 {
00379   index_sets.resize(n_blocks);
00380   std::vector<unsigned int> indices;
00381   unsigned int k = 0;
00382   for (ITERATOR cell = begin; cell != end; ++cell, ++k)
00383     {
00384       indices.resize(cell->get_fe().dofs_per_cell);
00385       cell->get_mg_dof_indices(indices);
00386       add(k, indices);
00387     }
00388 }
00389 
00390 
00391 template <typename ITERATOR>
00392 inline
00393 void
00394 BlockList::initialize(
00395   unsigned int n_blocks,
00396   const ITERATOR begin,
00397   const typename identity<ITERATOR>::type end,
00398   const std::vector<bool>& selected_dofs,
00399   unsigned int offset)
00400 {
00401   index_sets.resize(n_blocks);
00402   std::vector<unsigned int> indices;
00403   unsigned int k = 0;
00404   for (ITERATOR cell = begin; cell != end; ++cell, ++k)
00405     {
00406       indices.resize(cell->get_fe().dofs_per_cell);
00407       AssertDimension(selected_dofs.size(), cell->get_fe().dofs_per_cell);
00408 
00409       cell->get_dof_indices(indices);
00410       add(k, indices, selected_dofs);
00411     }
00412 }
00413 
00414 
00415 template <typename ITERATOR>
00416 inline
00417 void
00418 BlockList::initialize_mg(
00419   unsigned int n_blocks,
00420   const ITERATOR begin,
00421   const typename identity<ITERATOR>::type end,
00422   const std::vector<bool>& selected_dofs,
00423   unsigned int offset)
00424 {
00425   index_sets.resize(n_blocks);
00426   std::vector<unsigned int> indices;
00427   unsigned int k = 0;
00428   for (ITERATOR cell = begin; cell != end; ++cell, ++k)
00429     {
00430       indices.resize(cell->get_fe().dofs_per_cell);
00431       AssertDimension(selected_dofs.size(), cell->get_fe().dofs_per_cell);
00432 
00433       cell->get_mg_dof_indices(indices);
00434       add(k, indices, selected_dofs, offset);
00435     }
00436 }
00437 
00438 
00439 
00440 template <int dim, typename ITERATOR>
00441 inline
00442 bool
00443 BlockList::cell_generates_vertex_patch(
00444   const ITERATOR cell,
00445   bool same_level_only) const
00446 {
00447   switch(dim)
00448     {
00449       case 3:
00450             if (cell->at_boundary(4)) break;
00451             if (same_level_only && cell->neighbor(4)->level() != cell->level()) break;
00452             if (cell->neighbor(4)->at_boundary(0)) break;
00453             if (same_level_only && cell->neighbor(4)->neighbor(0)->level() != cell->level()) break;
00454             if (cell->neighbor(4)->at_boundary(2)) break;
00455             if (same_level_only && cell->neighbor(4)->neighbor(2)->level() != cell->level()) break;
00456             if (cell->neighbor(4)->neighbor(0)->at_boundary(2)) break;
00457             if (same_level_only && cell->neighbor(4)->neighbor(0)->neighbor(2)->level() != cell->level()) break;
00458                                              // No break here
00459       case 2:
00460             if (cell->at_boundary(2)) break;
00461             if (same_level_only && cell->neighbor(2)->level() != cell->level()) break;
00462             if (cell->neighbor(2)->at_boundary(0)) break;
00463             if (same_level_only && cell->neighbor(2)->neighbor(0)->level() != cell->level()) break;
00464       case 1:
00465             if (cell->at_boundary(0)) break;
00466             if (same_level_only && cell->neighbor(0)->level() != cell->level()) break;
00467             return true;
00468             break;
00469       default:
00470             Assert(false, ExcNotImplemented());
00471             break;
00472     }
00473   return false;
00474 }
00475 
00476 
00477 template <int dim, typename ITERATOR>
00478 inline
00479 unsigned int
00480 BlockList::count_vertex_patches(
00481   const ITERATOR begin,
00482   const typename identity<ITERATOR>::type end,
00483   bool same_level_only) const
00484 {
00485   unsigned int count = 0;
00486   for (ITERATOR cell = begin; cell != end; ++cell)
00487     if (cell_generates_vertex_patch<dim>(cell, same_level_only))
00488       ++count;
00489   return count;
00490 }
00491 
00492 
00493 template <int dim, typename ITERATOR>
00494 inline
00495 void
00496 BlockList::initialize_vertex_patches_mg(
00497   unsigned int n_blocks,
00498   const ITERATOR begin,
00499   const typename identity<ITERATOR>::type end,
00500   const std::vector<bool>& selected_dofs,
00501   unsigned int offset)
00502 {
00503   Assert(selected_dofs.size() == 0, ExcNotImplemented());
00504   Assert(offset==0, ExcNotImplemented());
00505   const FiniteElement<dim>& fe = begin->get_fe();
00506   Assert(fe.dofs_per_vertex == 0, ExcNotImplemented());
00507 
00508   index_sets.resize(n_blocks);
00509   std::vector<unsigned int> indices;
00510   unsigned int k = 0;
00511   for (ITERATOR cell = begin; cell != end; ++cell)
00512     {
00513       if (cell_generates_vertex_patch<dim>(cell, true))
00514         {
00515           indices.resize(cell->get_fe().dofs_per_cell);
00516 
00517           switch(dim)
00518             {
00519               case 3:
00520                     cell->neighbor(4)->get_mg_dof_indices(indices);
00521                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00522                       {
00523                         indices[fe.face_to_cell_index(i,1)] = numbers::invalid_unsigned_int;
00524                         indices[fe.face_to_cell_index(i,3)] = numbers::invalid_unsigned_int;
00525                         indices[fe.face_to_cell_index(i,4)] = numbers::invalid_unsigned_int;
00526                       }
00527                     add(k, indices);
00528                     cell->neighbor(4)->neighbor(0)->get_mg_dof_indices(indices);
00529                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00530                       {
00531                         indices[fe.face_to_cell_index(i,0)] = numbers::invalid_unsigned_int;
00532                         indices[fe.face_to_cell_index(i,3)] = numbers::invalid_unsigned_int;
00533                         indices[fe.face_to_cell_index(i,4)] = numbers::invalid_unsigned_int;
00534                       }
00535                     add(k, indices);
00536                     cell->neighbor(4)->neighbor(2)->get_mg_dof_indices(indices);
00537                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00538                       {
00539                         indices[fe.face_to_cell_index(i,1)] = numbers::invalid_unsigned_int;
00540                         indices[fe.face_to_cell_index(i,2)] = numbers::invalid_unsigned_int;
00541                         indices[fe.face_to_cell_index(i,4)] = numbers::invalid_unsigned_int;
00542                       }
00543                     add(k, indices);
00544                     cell->neighbor(4)->neighbor(2)->neighbor(0)->get_mg_dof_indices(indices);
00545                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00546                       {
00547                         indices[fe.face_to_cell_index(i,0)] = numbers::invalid_unsigned_int;
00548                         indices[fe.face_to_cell_index(i,2)] = numbers::invalid_unsigned_int;
00549                         indices[fe.face_to_cell_index(i,4)] = numbers::invalid_unsigned_int;
00550                       }
00551                     add(k, indices);
00552               case 2:
00553                     cell->neighbor(2)->get_mg_dof_indices(indices);
00554                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00555                       {
00556                         indices[fe.face_to_cell_index(i,1)] = numbers::invalid_unsigned_int;
00557                         indices[fe.face_to_cell_index(i,2)] = numbers::invalid_unsigned_int;
00558                         if (dim>2)
00559                           indices[fe.face_to_cell_index(i,5)] = numbers::invalid_unsigned_int;
00560                       }
00561                     add(k, indices);
00562                     cell->neighbor(2)->neighbor(0)->get_mg_dof_indices(indices);
00563                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00564                       {
00565                         indices[fe.face_to_cell_index(i,0)] = numbers::invalid_unsigned_int;
00566                         indices[fe.face_to_cell_index(i,2)] = numbers::invalid_unsigned_int;
00567                         if (dim>2)
00568                           indices[fe.face_to_cell_index(i,5)] = numbers::invalid_unsigned_int;
00569                       }
00570                     add(k, indices);
00571                                                      // no break here
00572               case 1:
00573                     cell->get_mg_dof_indices(indices);
00574                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00575                       {
00576                         indices[fe.face_to_cell_index(i,1)] = numbers::invalid_unsigned_int;
00577                         if (dim>1)
00578                           indices[fe.face_to_cell_index(i,3)] = numbers::invalid_unsigned_int;
00579                         if (dim>2)
00580                           indices[fe.face_to_cell_index(i,5)] = numbers::invalid_unsigned_int;
00581                       }
00582                     add(k, indices);
00583                     cell->neighbor(0)->get_mg_dof_indices(indices);
00584                     for (unsigned int i=0;i<fe.dofs_per_face;++i)
00585                       {
00586                         indices[fe.face_to_cell_index(i,0)] = numbers::invalid_unsigned_int;
00587                         if (dim>1)
00588                           indices[fe.face_to_cell_index(i,3)] = numbers::invalid_unsigned_int;
00589                         if (dim>2)
00590                           indices[fe.face_to_cell_index(i,5)] = numbers::invalid_unsigned_int;
00591                       }
00592                     add(k, indices);
00593                     break;
00594               default:
00595                     Assert(false, ExcNotImplemented());
00596                     break;
00597             }
00598           ++k;
00599         }
00600     }
00601 }
00602 
00603 
00604 inline
00605 unsigned int
00606 BlockList::size() const
00607 {
00608   return index_sets.size();
00609 }
00610 
00611 
00612 inline
00613 unsigned int
00614 BlockList::block_size(unsigned int block) const
00615 {
00616   return index_sets[block].size();
00617 }
00618 
00619 
00620 inline
00621 BlockList::const_iterator
00622 BlockList::begin(unsigned int block) const
00623 {
00624   AssertIndexRange(block, index_sets.size());
00625   return index_sets[block].begin();
00626 }
00627 
00628 
00629 inline
00630 BlockList::const_iterator
00631 BlockList::end(unsigned int block) const
00632 {
00633   AssertIndexRange(block, index_sets.size());
00634   return index_sets[block].end();
00635 }
00636 
00637 
00638 inline
00639 unsigned int
00640 BlockList::local_index(unsigned int block, unsigned int index) const
00641 {
00642   AssertIndexRange(block, index_sets.size());
00643   const block_container& b = index_sets[block];
00644   for (unsigned i=0;i<b.size();++i)
00645     if (b[i] == index)
00646       return i;
00647   return numbers::invalid_unsigned_int;
00648 }
00649 
00650 
00651 DEAL_II_NAMESPACE_CLOSE
00652 
00653 #endif
00654 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:39 by doxygen 1.7.3