disc » operators

p1operator.hh

Go to the documentation of this file.
00001 // $Id: p1operator.hh 529 2008-12-25 20:51:43Z sander $
00002 
00003 #ifndef DUNE_P1OPERATOR_HH
00004 #define DUNE_P1OPERATOR_HH
00005 
00006 #include<iostream>
00007 #include<vector>
00008 #include<set>
00009 #include<map>
00010 #include<stdio.h>
00011 #include<stdlib.h>
00012 
00013 #include<dune/common/timer.hh>
00014 #include<dune/common/fvector.hh>
00015 #include<dune/common/fmatrix.hh>
00016 #include<dune/common/exceptions.hh>
00017 #include<dune/common/geometrytype.hh>
00018 #include<dune/grid/common/grid.hh>
00019 #include<dune/grid/common/mcmgmapper.hh>
00020 //#include<dune/grid/utility/intersectiongetter.hh>
00021 #include<dune/istl/bvector.hh>
00022 #include<dune/istl/operators.hh>
00023 #include<dune/istl/bcrsmatrix.hh>
00024 #include<dune/disc/functions/p1function.hh> // for parallel extender class
00025 #include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
00026 #include"boundaryconditions.hh"
00027 #include"localstiffness.hh"
00028 
00045 namespace Dune
00046 {
00056   // make a type to sort matrix entries
00057   typedef std::pair<int,int> P1OperatorLink;
00058 
00059   // template meta program for inserting indices
00060   template<int n, int c>
00061   struct P1Operator_meta {
00062         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00063         static void addrowscube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00064                                                          const Refelem& refelem, Matrix& A, std::vector<bool>& visited, 
00065                                                          int hangingnodes, std::set<P1OperatorLink>& links)
00066         {
00067           if (refelem.type(0,0).isCube())
00068                 {
00069                   for (int i=0; i<refelem.size(c); i++) // loop over subentities of codim c of e
00070                         {
00071                           int index = allmapper.template map<c>(e,i);
00072                           if (!visited[index]) 
00073                                 {
00074                                   int corners = refelem.size(i,c,n);
00075                                   for (int j=0; j<corners/2; j++) // uses fact that diagonals are (0,corners-1), (1,corners-2) ...
00076                                         {
00077                                           int alpha = vertexmapper.template map<n>(e,refelem.subEntity(i,c,j,n));
00078                                           int beta = vertexmapper.template map<n>(e,refelem.subEntity(i,c,corners-1-j,n));
00079                                           A.incrementrowsize(alpha);
00080                                           A.incrementrowsize(beta);
00081                                           if (hangingnodes>0) // delete standard links
00082                                                 {
00083                                                   links.erase(P1OperatorLink(alpha,beta));
00084                                                   links.erase(P1OperatorLink(beta,alpha));
00085                                                 }
00086                                         }
00087                                   visited[index] = true;
00088                                 }
00089                         }
00090                 }
00091           if (refelem.type(0,0).isPyramid() && c==1)
00092                 {
00093                   int index = allmapper.template map<c>(e,0);
00094                   if (!visited[index]) 
00095                         {
00096                           int alpha = vertexmapper.template map<n>(e,0);
00097                           int beta = vertexmapper.template map<n>(e,2);
00098                           A.incrementrowsize(alpha);
00099                           A.incrementrowsize(beta);
00100                           if (hangingnodes>0) // delete standard links
00101                                 {
00102                                   links.erase(P1OperatorLink(alpha,beta));
00103                                   links.erase(P1OperatorLink(beta,alpha));
00104                                 }
00105                           alpha = vertexmapper.template map<n>(e,1);
00106                           beta = vertexmapper.template map<n>(e,3);
00107                           A.incrementrowsize(alpha);
00108                           A.incrementrowsize(beta);
00109                           if (hangingnodes>0) // delete standard links
00110                                 {
00111                                   links.erase(P1OperatorLink(alpha,beta));
00112                                   links.erase(P1OperatorLink(beta,alpha));
00113                                 }
00114                           visited[index] = true;
00115                         }
00116                 }         
00117           if (refelem.type(0,0).isPrism() && c==1)
00118                 {
00119                   int index = allmapper.template map<c>(e,1);
00120                   if (!visited[index]) 
00121                         {
00122                           int alpha = vertexmapper.template map<n>(e,0);
00123                           int beta = vertexmapper.template map<n>(e,4);
00124                           A.incrementrowsize(alpha);
00125                           A.incrementrowsize(beta);
00126                           if (hangingnodes>0) // delete standard links
00127                                 {
00128                                   links.erase(P1OperatorLink(alpha,beta));
00129                                   links.erase(P1OperatorLink(beta,alpha));
00130                                 }
00131                           alpha = vertexmapper.template map<n>(e,1);
00132                           beta = vertexmapper.template map<n>(e,3);
00133                           A.incrementrowsize(alpha);
00134                           A.incrementrowsize(beta);
00135                           if (hangingnodes>0) // delete standard links
00136                                 {
00137                                   links.erase(P1OperatorLink(alpha,beta));
00138                                   links.erase(P1OperatorLink(beta,alpha));
00139                                 }
00140                           visited[index] = true;
00141                         }
00142                   index = allmapper.template map<c>(e,2);
00143                   if (!visited[index]) 
00144                         {
00145                           int alpha = vertexmapper.template map<n>(e,1);
00146                           int beta = vertexmapper.template map<n>(e,5);
00147                           A.incrementrowsize(alpha);
00148                           A.incrementrowsize(beta);
00149                           if (hangingnodes>0) // delete standard links
00150                                 {
00151                                   links.erase(P1OperatorLink(alpha,beta));
00152                                   links.erase(P1OperatorLink(beta,alpha));
00153                                 }
00154                           alpha = vertexmapper.template map<n>(e,2);
00155                           beta = vertexmapper.template map<n>(e,4);
00156                           A.incrementrowsize(alpha);
00157                           A.incrementrowsize(beta);
00158                           if (hangingnodes>0) // delete standard links
00159                                 {
00160                                   links.erase(P1OperatorLink(alpha,beta));
00161                                   links.erase(P1OperatorLink(beta,alpha));
00162                                 }
00163                           visited[index] = true;
00164                         }
00165                   index = allmapper.template map<c>(e,3);
00166                   if (!visited[index]) 
00167                         {
00168                           int alpha = vertexmapper.template map<n>(e,0);
00169                           int beta = vertexmapper.template map<n>(e,5);
00170                           A.incrementrowsize(alpha);
00171                           A.incrementrowsize(beta);
00172                           if (hangingnodes>0) // delete standard links
00173                                 {
00174                                   links.erase(P1OperatorLink(alpha,beta));
00175                                   links.erase(P1OperatorLink(beta,alpha));
00176                                 }
00177                           alpha = vertexmapper.template map<n>(e,2);
00178                           beta = vertexmapper.template map<n>(e,3);
00179                           A.incrementrowsize(alpha);
00180                           A.incrementrowsize(beta);
00181                           if (hangingnodes>0) // delete standard links
00182                                 {
00183                                   links.erase(P1OperatorLink(alpha,beta));
00184                                   links.erase(P1OperatorLink(beta,alpha));
00185                                 }
00186                           visited[index] = true;
00187                         }
00188                 }         
00189           P1Operator_meta<n,c-1>::addrowscube(e,vertexmapper,allmapper,refelem,A,visited,hangingnodes,links);
00190           return;
00191         }
00192         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00193         static void addindicescube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00194                                    const Refelem& refelem, Matrix& A, std::vector<bool>& visited)
00195         {
00196           if (refelem.type(0,0).isCube())
00197                 {
00198                   for (int i=0; i<refelem.size(c); i++)
00199                         {
00200                           int index = allmapper.template map<c>(e,i);
00201                           if (!visited[index]) 
00202                                 {
00203                                   int corners = refelem.size(i,c,n);
00204                                   for (int j=0; j<corners/2; j++) // uses fact that diagonals are (0,corners-1), (1,corners-2) ...
00205                                         {
00206                                           int alpha = vertexmapper.template map<n>(e,refelem.subEntity(i,c,j,n));
00207                                           int beta = vertexmapper.template map<n>(e,refelem.subEntity(i,c,corners-1-j,n));
00208                                           A.addindex(alpha,beta);
00209                                           A.addindex(beta,alpha);
00210                                         }
00211                                   visited[index] = true;
00212                                 }
00213                         }
00214                 }
00215           if (refelem.type(0,0).isPyramid() && c==1)
00216                 {
00217                   int index = allmapper.template map<c>(e,0);
00218                   if (!visited[index]) 
00219                         {
00220                           int alpha = vertexmapper.template map<n>(e,0);
00221                           int beta = vertexmapper.template map<n>(e,2);
00222                           A.addindex(alpha,beta);
00223                           A.addindex(beta,alpha);
00224                           alpha = vertexmapper.template map<n>(e,1);
00225                           beta = vertexmapper.template map<n>(e,3);
00226                           A.addindex(alpha,beta);
00227                           A.addindex(beta,alpha);
00228                           visited[index] = true;
00229                         }
00230                 }         
00231           if (refelem.type(0,0).isPrism() && c==1)
00232                 {
00233                   int index = allmapper.template map<c>(e,1);
00234                   if (!visited[index]) 
00235                         {
00236                           int alpha = vertexmapper.template map<n>(e,0);
00237                           int beta = vertexmapper.template map<n>(e,4);
00238                           A.addindex(alpha,beta);
00239                           A.addindex(beta,alpha);
00240                           alpha = vertexmapper.template map<n>(e,1);
00241                           beta = vertexmapper.template map<n>(e,3);
00242                           A.addindex(alpha,beta);
00243                           A.addindex(beta,alpha);
00244                           visited[index] = true;
00245                         }
00246                   index = allmapper.template map<c>(e,2);
00247                   if (!visited[index]) 
00248                         {
00249                           int alpha = vertexmapper.template map<n>(e,1);
00250                           int beta = vertexmapper.template map<n>(e,5);
00251                           A.addindex(alpha,beta);
00252                           A.addindex(beta,alpha);
00253                           alpha = vertexmapper.template map<n>(e,2);
00254                           beta = vertexmapper.template map<n>(e,4);
00255                           A.addindex(alpha,beta);
00256                           A.addindex(beta,alpha);
00257                           visited[index] = true;
00258                         }
00259                   index = allmapper.template map<c>(e,3);
00260                   if (!visited[index]) 
00261                         {
00262                           int alpha = vertexmapper.template map<n>(e,0);
00263                           int beta = vertexmapper.template map<n>(e,5);
00264                           A.addindex(alpha,beta);
00265                           A.addindex(beta,alpha);
00266                           alpha = vertexmapper.template map<n>(e,2);
00267                           beta = vertexmapper.template map<n>(e,3);
00268                           A.addindex(alpha,beta);
00269                           A.addindex(beta,alpha);
00270                           visited[index] = true;
00271                         }
00272                 }         
00273           P1Operator_meta<n,c-1>::addindicescube(e,vertexmapper,allmapper,refelem,A,visited);
00274           return;
00275         }
00276   };
00277   template<int n>
00278   struct P1Operator_meta<n,0> {
00279         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00280         static void addrowscube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00281                                                          const Refelem& refelem, Matrix& A, std::vector<bool>& visited,
00282                                                          int hangingnodes, std::set<P1OperatorLink>& links)
00283         {
00284           if (!refelem.type(0,0).isCube()) return;
00285           int corners = refelem.size(n);
00286           for (int j=0; j<corners/2; j++) // uses fact that diagonals are (0,corners-1), (1,corners-2) ...
00287                 {
00288                   int alpha = vertexmapper.template map<n>(e,refelem.subEntity(0,0,j,n));
00289                   int beta = vertexmapper.template map<n>(e,refelem.subEntity(0,0,corners-1-j,n));
00290                   A.incrementrowsize(alpha);
00291                   A.incrementrowsize(beta);
00292                   if (hangingnodes>0) // delete standard links
00293                         {
00294                           links.erase(P1OperatorLink(alpha,beta));
00295                           links.erase(P1OperatorLink(beta,alpha));
00296                         }
00297                 }
00298           return;
00299         }
00300         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00301         static void addindicescube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00302                                   const Refelem& refelem, Matrix& A, std::vector<bool>& visited)
00303         {
00304           if (!refelem.type(0,0).isCube()) return;
00305           int corners = refelem.size(n);
00306           for (int j=0; j<corners/2; j++) // uses fact that diagonals are (0,corners-1), (1,corners-2) ...
00307                 {
00308                   int alpha = vertexmapper.template map<n>(e,refelem.subEntity(0,0,j,n));
00309                   int beta = vertexmapper.template map<n>(e,refelem.subEntity(0,0,corners-1-j,n));
00310                   A.addindex(alpha,beta);
00311                   A.addindex(beta,alpha);
00312                 }
00313           return;
00314         }
00315   };
00316 
00317   // handles case n=1, c=-1
00318   template<int n>
00319   struct P1Operator_meta<n,-1> {
00320         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00321         static void addrowscube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00322                                                          const Refelem& refelem, Matrix& A, std::vector<bool>& visited,
00323                                                          int hangingnodes, std::set<P1OperatorLink>& links)
00324         {
00325           return;
00326         }
00327         template<class Entity, class VMapper, class AMapper, class Refelem, class Matrix>
00328         static void addindicescube (const Entity& e, const VMapper& vertexmapper, const AMapper& allmapper, 
00329                                   const Refelem& refelem, Matrix& A, std::vector<bool>& visited)
00330         {
00331           return;
00332         }
00333   };
00334 
00335     template<typename Grid, typename RT, bool hasHangingNodes>
00336     struct HangingNodesIdentifier
00337     {};
00338     
00339     template<typename Grid,typename RT>
00340     struct HangingNodesIdentifier<Grid,RT,false>
00341     {
00342       template<typename IndexSet, typename VM>
00343       inline void process(std::vector<unsigned char> S, std::vector<bool>& hanging, std::set<P1OperatorLink>& links, IndexSet& indexset, VM& vertexmapper) const
00344       {}
00345       
00346       inline std::size_t count() const
00347       {
00348         return 0;
00349       }
00350     };
00351 
00352     template<typename Grid,typename RT>
00353     class HangingNodesIdentifier<Grid,RT,true>
00354     {
00355     public:
00356       typedef typename Grid::ctype DT;
00357       enum {n=Grid::dimension};
00358       typedef typename Grid::template Codim<0>::EntityPointer EEntityPointer;
00359 
00360       template<typename GridView, typename VM>
00361       void process(std::vector<unsigned char> S, std::vector<bool>& hanging, std::set<P1OperatorLink>& links, GridView& gridView, VM& vertexmapper)
00362       {
00363         
00364         typedef typename GridView::template Codim<0>::Iterator Iterator;
00365         // LOOP 1 : Prepare hanging node detection
00366         Iterator eendit = gridView.template end<0>();
00367         for (Iterator it = gridView.template begin<0>(); it!=eendit; ++it)
00368           {
00369             Dune::GeometryType gt = it->type();
00370             const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
00371               refelem = ReferenceElements<DT,n>::general(gt);
00372             
00373             // compute S value in vertex
00374             for (int i=0; i<refelem.size(n); i++)
00375               {
00376                 int alpha = vertexmapper.template map<n>(*it,i);
00377                 if (S[alpha]>it->level()) S[alpha] = it->level(); // compute minimum
00378               }
00379           }
00380 
00381         // LOOP 2 : second stage of detecting hanging nodes
00382         for (Iterator it = gridView.template begin<0>(); it!=eendit; ++it)
00383           {
00384             Dune::GeometryType gt = it->type();
00385             const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
00386               refelem = ReferenceElements<DT,n>::general(gt);
00387             
00388             // detect hanging nodes
00389             typename GridView::IntersectionIterator endiit = gridView.iend(*it);
00390             for (typename GridView::IntersectionIterator iit = gridView.ibegin(*it); iit!=endiit; ++iit)
00391               if (iit->neighbor())
00392                 {
00393                   // check if neighbor is on lower level
00394                   const EEntityPointer outside = iit->outside();
00395                   if (it->level()<=outside->level()) continue;
00396                   
00397                   // loop over all vertices of this face
00398                   for (int j=0; j<refelem.size(iit->numberInSelf(),1,n); j++)
00399                     {
00400                       int alpha = vertexmapper.template map<n>(*it,refelem.subEntity(iit->numberInSelf(),1,j,n));
00401                       if (S[alpha]==it->level()) 
00402                         hanging[alpha] = true;
00403                     }
00404                 }
00405           }
00406         
00407         // local to global maps
00408         int l2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
00409         int fl2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
00410         
00411         // LOOP 3 : determine additional links due to hanging nodes
00412         for (Iterator it = gridView.template begin<0>(); it!=eendit; ++it)
00413           {
00414             Dune::GeometryType gt = it->type();
00415             const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
00416               refelem = ReferenceElements<DT,n>::general(gt);
00417             
00418             // build local to global map
00419             bool hasHangingNodes = false; // flag set to true if this element has hanging nodes
00420             for (int i=0; i<refelem.size(n); i++)
00421               {
00422                 l2g[i] = vertexmapper.template map<n>(*it,i);
00423                 if (hanging[l2g[i]]) hasHangingNodes=true;
00424               }
00425             if (!hasHangingNodes) continue;
00426             
00427             // handle father element if hanging nodes were detected
00428             // get father element
00429             const EEntityPointer father = it->father();
00430             
00431             // build local to global map for father
00432             for (int i=0; i<refelem.size(n); i++)
00433               fl2g[i] = vertexmapper.template map<n>(*father,i);
00434             
00435             // a map that inverts l2g
00436             std::map<int,int> g2l;
00437             for (int i=0; i<refelem.size(n); i++)
00438               g2l[l2g[i]] = i;
00439             
00440             // connect all fine nodes to all coarse nodes
00441             for (int i=0; i<refelem.size(n); i++) // nodes in *it
00442               for (int j=0; j<refelem.size(n); j++) // nodes in *father
00443                 if (g2l.find(fl2g[j])==g2l.end())
00444                   {
00445                     links.insert(P1OperatorLink(l2g[i],fl2g[j]));
00446                     links.insert(P1OperatorLink(fl2g[j],l2g[i]));
00447                   }
00448           }
00449         // count hanging nodes
00450         hangingnodes = 0;
00451         for (int i=0; i<vertexmapper.size(); i++) 
00452           if (hanging[i]) hangingnodes++;
00453         
00454       }
00455       inline std::size_t count() const
00456       {
00457         return hangingnodes;
00458       }
00459     private:
00460       std::size_t hangingnodes;
00461     };
00462     
00463 
00464   
00475   template<class G, class RT, class GV, class LC, int m=1>
00476   class P1OperatorBase
00477   {
00478   public:    
00479         // export type used to store the matrix
00480         typedef FieldMatrix<RT,m,m> BlockType; 
00481         typedef BCRSMatrix<BlockType> RepresentationType;
00482 
00483         // mapper: one data element per vertex
00484         template<int dim>
00485         struct P1Layout
00486         {
00487           bool contains (Dune::GeometryType gt)
00488           {
00489               return gt.dim() == 0;
00490           }
00491         }; 
00492 
00493         // mapper: one data element in every entity
00494         template<int dim>
00495         struct AllLayout
00496         {
00497           bool contains (Dune::GeometryType gt)
00498           {
00499                 return true;
00500           }
00501         }; 
00502 
00503         typedef typename G::ctype DT;
00504         enum {n=G::dimension};
00505       typedef typename GV::IndexSet IS;
00506         typedef typename G::template Codim<0>::Entity Entity;
00507         typedef typename GV::template Codim<0>::Iterator Iterator;
00508         typedef typename GV::template Codim<n>::Iterator VIterator;
00509         typedef typename G::template Codim<0>::EntityPointer EEntityPointer;
00510         typedef typename G::Traits::GlobalIdSet IDS;
00511         typedef typename IDS::IdType IdType;
00512         typedef std::set<IdType> GIDSet;
00513         typedef MultipleCodimMultipleGeomTypeMapper<G,IS,P1Layout> VM;
00514         typedef MultipleCodimMultipleGeomTypeMapper<G,IS,AllLayout> AM;
00515 
00516   private:
00517 
00518         // a function to compute the number of nonzeros
00519         // does not work for prisms and pyramids yet ?!
00520         int nnz (const IS& is)
00521         {
00522           int s = 0;
00523           s += is.size(n);   // vertices
00524           s += 2*is.size(n-1); // edges
00525 
00526           for (int c=0; c<n-1; c++)
00527                 {
00528                   s += 2*is.size(GeometryType(GeometryType::cube,G::dimension-c))*(1<<(n-c-1));
00529                 }
00530 
00531           // hanging node correction
00532           s += links.size();
00533 
00534           return s;
00535         }
00536       
00537     
00538         // extra initialization function
00539         // 0) This method is executed before matrix is allocated
00540         // 1) determine hanging nodes as described in the paper
00541         // 2) generate a set with additional links
00542         //    The standard links are deleted later on
00543         bool init (const G& g, const GV& gridView, LC lc, bool extendoverlap)
00544         {
00545           // parallel stuff we need to know
00546           if (extendoverlap && g.overlapSize(0)>0)
00547                 DUNE_THROW(GridError,"P1OperatorBase: extending overlap requires nonoverlapping grid");
00548           extendOverlap = extendoverlap;
00549           extraDOFs = 0;
00550 
00551           // resize the S vector needed for detecting hanging nodes
00552           watch.reset();
00553           // resize hanging and initialize with false
00554           hanging.resize(vertexmapper.size(), false);
00555           // the number of levels never exceeds 100 ...
00556           std::vector<unsigned char> S(vertexmapper.size(), 100);
00557 
00558           // Detect the hanging nodes
00559           typedef HangingNodesIdentifier<G,RT,Dune::Capabilities::template hasHangingNodes<G>::v>
00560             HangingNodesIdentifier;
00561           HangingNodesIdentifier hangingNodes;
00562           hangingNodes.process(S,hanging,links,gridView,vertexmapper);
00563           
00564           // count hanging nodes
00565           hangingnodes = hangingNodes.count();
00566 
00567 //        std::cout << "=== P1OperatorBase hanging node detection + add links " <<  watch.elapsed() << std::endl;
00568 
00569           // compute additional links due to extended overlap
00570           watch.reset();
00571           if (extendOverlap)
00572                 {
00573                   // set of neighbors in global ids for border vertices
00574                   std::map<int,GIDSet> borderlinks;
00575 
00576                   // compute extension
00577                   P1ExtendOverlap<G,GV,VM,LC> extender(lc);
00578                   extender.extend(g,gv,vertexmapper,borderlinks,extraDOFs,gid2index);
00579 
00580                   // put in extra links due to overlap 
00581                   // loop over all neighbors of border vertices
00582                   for (typename std::map<int,GIDSet>::iterator i=borderlinks.begin(); i!=borderlinks.end(); ++i)
00583                         for (typename GIDSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
00584                           links.insert(P1OperatorLink(i->first,gid2index[*j]));
00585                   
00586                   // insert diagonal links for extra DOFs
00587                   for (int i=0; i<extraDOFs; i++)
00588                         links.insert(P1OperatorLink(vertexmapper.size()+i,vertexmapper.size()+i));                
00589                 }
00590 
00591           // Note: links contains now also connections that are standard.
00592           // So below we have throw out these connections again!
00593 //        std::cout << "=== P1OperatorBase parallel extend overlap " <<  watch.elapsed() << std::endl;
00594 
00595           return true;
00596         }
00597 
00598 
00599         // return number of rows/columns
00600         int size () const
00601         {
00602           return vertexmapper.size()+extraDOFs;
00603         }
00604 
00605         struct MatEntry
00606         {
00607           IdType first;
00608           BlockType second;
00609           MatEntry (const IdType& f, const BlockType& s) : first(f),second(s) {}
00610           MatEntry () {}
00611         };
00612 
00613         // A DataHandle class to exchange matrix entries
00614         class MatEntryExchange : 
00615           public CommDataHandleIF<MatEntryExchange,MatEntry> {
00616           typedef typename RepresentationType::RowIterator rowiterator;
00617           typedef typename RepresentationType::ColIterator coliterator;
00618         public:
00620           typedef MatEntry DataType;
00621 
00623           bool contains (int dim, int codim) const
00624           {
00625                 return (codim==dim);
00626           }
00627 
00629           bool fixedsize (int dim, int codim) const
00630           {
00631                 return false;
00632           }
00633 
00638           template<class EntityType>
00639           size_t size (EntityType& e) const
00640           {
00641                 int i=vertexmapper.map(e);
00642                 int n=0;
00643                 for (coliterator j=A[i].begin(); j!=A[i].end(); ++j)
00644                   n++;
00645                 return n;
00646           }
00647 
00649           template<class MessageBuffer, class EntityType>
00650           void gather (MessageBuffer& buff, const EntityType& e) const
00651           {
00652                 int i=vertexmapper.map(e);
00653                 for (coliterator j=A[i].begin(); j!=A[i].end(); ++j)
00654                   {
00655                         typename std::map<int,IdType>::const_iterator it=index2gid.find(j.index());
00656                         if (it==index2gid.end())
00657                           DUNE_THROW(GridError,"MatEntryExchange::gather(): index not in map");
00658                         buff.write(MatEntry(it->second,*j));
00659                   }
00660           }
00661 
00666           template<class MessageBuffer, class EntityType>
00667           void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
00668           {
00669                 int i=vertexmapper.map(e);
00670                 for (size_t k=0; k<n; k++)
00671                   {
00672                         MatEntry m;
00673                         buff.read(m);
00674                         typename std::map<IdType,int>::const_iterator it=gid2index.find(m.first);
00675                         if (it==gid2index.end())
00676                           DUNE_THROW(GridError,"MatEntryExchange::scatter(): gid not in map");
00677                         A[i][it->second] += m.second;
00678                   }
00679           }
00680 
00682           MatEntryExchange (const G& g, const std::map<IdType,int>& g2i, 
00683                                                 const std::map<int,IdType>& i2g,
00684                                                 const VM& vm,
00685                                                 RepresentationType& a)
00686                 : grid(g), gid2index(g2i), index2gid(i2g), vertexmapper(vm), A(a) 
00687           {}
00688  
00689         private:
00690           const G& grid;
00691           const std::map<IdType,int>& gid2index;
00692           const std::map<int,IdType>& index2gid;
00693           const VM& vertexmapper;
00694           RepresentationType& A;
00695         };
00696 
00697   public:
00698 
00699         P1OperatorBase (const G& g, const GV& gridView, LC lcomm, bool extendoverlap=false) 
00700             : grid(g),gv(gridView),lc(lcomm),vertexmapper(g,gridView.indexSet()),allmapper(g,gridView.indexSet()),links(),
00701               initialized(init(g,gridView,lc,extendoverlap)),A(size(),size(),nnz(gridView.indexSet()),RepresentationType::random)
00702                 
00703         {
00704           // be verbose
00705 //        std::cout << g.comm().rank() << ": " << "vector size = " << vertexmapper.size() << " + " << extraDOFs << std::endl;
00706           std::cout << g.comm().rank() << ": " << "making " << size() << "x" 
00707                     << size() << " matrix with " << nnz(gv.indexSet()) << " nonzeros" << std::endl;
00708 //        std::cout << g.comm().rank() << ": " << "allmapper has size " << allmapper.size() << std::endl;
00709 //        std::cout << g.comm().rank() << ": " << "vertexmapper has size " << vertexmapper.size() << std::endl;
00710 //        std::cout << g.comm().rank() << ": " << "hanging nodes=" << hangingnodes << " links=" << links.size() << std::endl;
00711    
00712           // set size of all rows to zero
00713           for (size_t i=0; i<gv.indexSet().size(n); i++)
00714                 A.setrowsize(i,0); 
00715 
00716           // build needs a flag for all entities of all codims
00717           std::vector<bool> visited(allmapper.size());
00718           for (int i=0; i<allmapper.size(); i++) visited[i] = false;
00719 
00720           // LOOP 4 : Compute row sizes
00721           watch.reset();
00722           Iterator eendit = gv.template end<0>();
00723           for (Iterator it = gv.template begin<0>(); it!=eendit; ++it)
00724                 {
00725                   Dune::GeometryType gt = it->type();
00726                   const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
00727                         refelem = ReferenceElements<DT,n>::general(gt);
00728 
00729                   // vertices, c=n
00730                   for (int i=0; i<refelem.size(n); i++)
00731                         {
00732                           int index = allmapper.template map<n>(*it,i);
00733                           int alpha = vertexmapper.template map<n>(*it,i);
00734                           //                      std::cout << "index=" << index << " alpha=" << alpha << std::endl;
00735                           if (!visited[index]) 
00736                                 {
00737                                   A.incrementrowsize(alpha);
00738                                   visited[index] = true;
00739 //                                printf("increment row %04d\n",alpha);
00740                                 }
00741                         }
00742 
00743                   // edges for all element types, c=n-1
00744                   for (int i=0; i<refelem.size(n-1); i++)
00745                         {
00746                           int index = allmapper.template map<n-1>(*it,i);
00747                           int alpha = vertexmapper.template map<n>(*it,refelem.subEntity(i,n-1,0,n));
00748                           int beta = vertexmapper.template map<n>(*it,refelem.subEntity(i,n-1,1,n));
00749                           if (!visited[index]) 
00750                                 {
00751                                   A.incrementrowsize(alpha);
00752                                   A.incrementrowsize(beta);
00753                                   visited[index] = true;
00754                                   if (hangingnodes>0 || extendOverlap) // delete standard links
00755                                         {
00756                                           links.erase(P1OperatorLink(alpha,beta));
00757                                           links.erase(P1OperatorLink(beta,alpha));
00758                                         }
00759                                 }
00760                         }
00761 
00762                   // for codim n-2 to 0 we need a template metaprogram
00763                   if (!gt.isSimplex())
00764                         P1Operator_meta<n,n-2>::addrowscube(*it,vertexmapper,allmapper,refelem,A,visited,hangingnodes+(extendOverlap),links);
00765                 }
00766 
00767           // additional links due to hanging nodes
00768           for (typename std::set<P1OperatorLink>::iterator i=links.begin(); i!=links.end(); ++i)
00769                 A.incrementrowsize(i->first);
00770 
00771           // now the row sizes have been set
00772           A.endrowsizes();
00773 //        std::cout << "=== P1OperatorBase compute row sizes " <<  watch.elapsed() << std::endl;
00774 
00775           // clear the flags for the next round, actually that is not necessary because addindex takes care of this
00776           for (int i=0; i<allmapper.size(); i++) visited[i] = false;
00777 
00778           // LOOP 5 : insert the nonzeros
00779           watch.reset();
00780           for (Iterator it = gv.template begin<0>(); it!=eendit; ++it)
00781                 {
00782                   Dune::GeometryType gt = it->type();
00783                   const typename Dune::ReferenceElementContainer<DT,n>::value_type&
00784                         refelem = ReferenceElements<DT,n>::general(gt);
00785                   //              std::cout << "ELEM " << GeometryName(gt) << std::endl;
00786 
00787                   // vertices, c=n
00788                   for (int i=0; i<refelem.size(n); i++)
00789                         {
00790                           int index = allmapper.template map<n>(*it,i);
00791                           int alpha = vertexmapper.template map<n>(*it,i);
00792                           //                      std::cout << "vertex allindex " << index << std::endl; 
00793                           if (!visited[index]) 
00794                                 {
00795                                   A.addindex(alpha,alpha);
00796                                   visited[index] = true;
00797                                 }
00798                         }
00799 
00800                   // edges for all element types, c=n-1
00801                   for (int i=0; i<refelem.size(n-1); i++)
00802                         {
00803                           int index = allmapper.template map<n-1>(*it,i);
00804                           //                      std::cout << "edge allindex " << index << std::endl; 
00805                           if (!visited[index]) 
00806                                 {
00807                                   int alpha = vertexmapper.template map<n>(*it,refelem.subEntity(i,n-1,0,n));
00808                                   int beta = vertexmapper.template map<n>(*it,refelem.subEntity(i,n-1,1,n));
00809                                   A.addindex(alpha,beta);
00810                                   A.addindex(beta,alpha);
00811                                   visited[index] = true;
00812 //                                printf("adding (%04d,%04d) index=%04d\n",alpha,beta,index);
00813 //                                printf("adding (%04d,%04d) index=%04d\n",beta,alpha,index);
00814                                 }
00815                         }
00816 
00817                   // for codim n-2 to 0 we need a template metaprogram
00818                   if (!gt.isSimplex())
00819                         P1Operator_meta<n,n-2>::addindicescube(*it,vertexmapper,allmapper,refelem,A,visited);
00820                 }
00821 
00822           // additional links due to hanging nodes
00823           for (typename std::set<P1OperatorLink>::iterator i=links.begin(); i!=links.end(); ++i)
00824                 A.addindex(i->first,i->second);
00825 
00826           // now the matrix is ready for use
00827           A.endindices();
00828 //        std::cout << "=== P1OperatorBase index insertion " <<  watch.elapsed() << std::endl;
00829 
00830           // delete additional links
00831           links.clear();
00832 
00833           //      std::cout << grid.comm().rank() << ": " << "matrix initialized" << std::endl;
00834         }
00835 
00837         const RepresentationType& operator* () const
00838         {
00839           return A;
00840         }
00841 
00843         RepresentationType& operator* ()
00844         {
00845           return A;
00846         }
00847 
00849         void sumEntries ()
00850         {
00851           if (!extendOverlap) return;
00852 
00853           // build forward map
00854           std::map<int,IdType> index2gid;
00855           for (typename std::map<IdType,int>::iterator i=gid2index.begin(); i!=gid2index.end(); ++i)
00856                 index2gid[i->second] = i->first;
00857 
00858           // communicate matrix entries
00859           MatEntryExchange datahandle(grid,gid2index,index2gid,vertexmapper,A);
00860           lc.template communicate<MatEntryExchange>(datahandle,
00861                                                                                                   InteriorBorder_InteriorBorder_Interface,
00862                                                                                                   ForwardCommunication);
00863         }
00864 
00865   protected:
00866         Timer watch;
00867         const G& grid;  
00868         const GV gv;
00869         LC lc;
00870         VM vertexmapper;
00871         AM allmapper;
00872         std::vector<bool> hanging;
00873         std::set<P1OperatorLink> links; 
00874         int hangingnodes;
00875         bool extendOverlap;
00876         int extraDOFs;
00877         std::map<IdType,int> gid2index;
00878         bool initialized;
00879         RepresentationType A;
00880   };
00881 
00882 
00883 
00884 
00888   template<class G, class RT, class GV, class LC, int m>
00889   class P1OperatorAssembler : public P1OperatorBase<G,RT,GV,LC,m>
00890   {
00891         typedef typename G::ctype DT;
00892         enum {n=G::dimension};
00893       typedef typename GV::IndexSet IS;
00894         typedef typename G::template Codim<0>::Entity Entity;
00895         typedef typename GV::template Codim<0>::Iterator Iterator;
00896         typedef typename GV::template Codim<n>::Iterator VIterator;
00897         typedef typename G::template Codim<0>::HierarchicIterator HierarchicIterator;
00898         typedef typename G::template Codim<0>::EntityPointer EEntityPointer;
00899         typedef typename P1Function<GV,RT,LC,m>::RepresentationType VectorType;
00900         typedef typename VectorType::block_type VBlockType;
00901         typedef typename P1OperatorBase<G,RT,GV,LC,m>::RepresentationType MatrixType;
00902         typedef typename MatrixType::field_type MFieldType;
00903         typedef typename MatrixType::block_type MBlockType;
00904         typedef typename MatrixType::RowIterator rowiterator;
00905         typedef typename MatrixType::ColIterator coliterator;
00906         typedef typename P1OperatorBase<G,RT,GV,LC,m>::VM VM;
00907     typedef array<BoundaryConditions::Flags,m> BCBlockType;     // componentwise boundary conditions
00908 
00909 
00910         // A DataHandle class to exchange border rows
00911         class AccumulateBCFlags : 
00912       public CommDataHandleIF<AccumulateBCFlags,std::pair<BCBlockType,VBlockType> > {
00913         public:
00915           typedef std::pair<BCBlockType,VBlockType> DataType; // BC flags are stored in a char per vertex
00916 
00918           bool contains (int dim, int codim) const
00919           {
00920                 return (codim==dim);
00921           }
00922 
00924           bool fixedsize (int dim, int codim) const
00925           {
00926                 return true;
00927           }
00928 
00933           template<class EntityType>
00934           size_t size (EntityType& e) const
00935           {
00936                 return 1;
00937           }
00938 
00940           template<class MessageBuffer, class EntityType>
00941           void gather (MessageBuffer& buff, const EntityType& e) const
00942           {
00943                 int alpha = vertexmapper.map(e);
00944                 buff.write(DataType(essential[alpha],f[alpha])); 
00945           }
00946 
00951           template<class MessageBuffer, class EntityType>
00952           void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
00953           {
00954                 DataType x;
00955                 buff.read(x);
00956                 int alpha = vertexmapper.map(e);
00957                 for (int i=0; i<m; i++)
00958                   if (x.first[i]>essential[alpha][i])
00959                         {
00960                           essential[alpha][i] = x.first[i];
00961                           f[alpha][i] = x.second[i];
00962                         }
00963           }
00964 
00966           AccumulateBCFlags (const G& g, const VM& vm, std::vector<BCBlockType>& e, VectorType& _f) 
00967                 : grid(g), essential(e), vertexmapper(vm), f(_f)
00968           {}
00969  
00970         private:
00971           const G& grid;
00972           std::vector<BCBlockType>& essential;
00973           const VM& vertexmapper;
00974           VectorType& f;
00975         };
00976 
00977 
00978 
00979   public:
00980         P1OperatorAssembler (const G& g, const GV& gridView, LC lcomm, bool extendoverlap=false)
00981           : P1OperatorBase<G,RT,GV,LC,m>(g,gridView,lcomm,extendoverlap)
00982         {       }
00983 
00984 
01005         void assemble (LocalStiffness<GV,RT,m>& loc, P1Function<GV,RT,LC,m>& u, P1Function<GV,RT,LC,m>& f)
01006         {
01007           // check size
01008           if ((*u).N()!=this->A.M() || (*f).N()!=this->A.N())
01009                 DUNE_THROW(MathError,"P1OperatorAssembler::assemble(): size mismatch");         
01010 
01011           // clear global stiffness matrix and right hand side
01012           this->watch.reset();
01013           this->A = 0;
01014           *f = 0;
01015 //        std::cout << "=== P1OperatorBase clear matrix " <<  this->watch.elapsed() << std::endl;
01016 
01017           // allocate flag vector to hold flags for essential boundary conditions
01018           std::vector<BCBlockType> essential(this->vertexmapper.size());
01019           for (typename std::vector<BCBlockType>::size_type i=0; i<essential.size(); i++)
01020                 essential[i].assign(BoundaryConditions::neumann);
01021 
01022           // allocate flag vector to note hanging nodes whose row has been assembled
01023           std::vector<unsigned char> treated(this->vertexmapper.size());
01024           for (std::vector<unsigned char>::size_type
01025              i=0; i<treated.size(); i++) treated[i] = false;
01026 
01027           // hanging node stuff
01028           RT alpha[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize][Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
01029           
01030           // local to global id mapping (do not ask vertex mapper repeatedly
01031           int l2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
01032           int fl2g[Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize];
01033 
01034           // run over all leaf elements
01035           Iterator eendit = this->gv.template end<0>();
01036           for (Iterator it = this->gv.template begin<0>(); it!=eendit; ++it)
01037                 {
01038                   // parallelization
01039                   if (it->partitionType()==GhostEntity)
01040                         continue;
01041 
01042                   // get access to shape functions for P1 elements
01043                   Dune::GeometryType gt = it->type();
01044                   const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
01045                         sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1);
01046 
01047                   // get local to global id map
01048                   for (int k=0; k<sfs.size(); k++)
01049                         {
01050                           if (sfs[k].codim()!=n) DUNE_THROW(MathError,"expected codim==dim");
01051                           int alpha = this->vertexmapper.template map<n>(*it,sfs[k].entity());
01052                           l2g[k] = alpha;
01053                         }
01054 
01055                   // build local stiffness matrix for P1 elements
01056                   // inludes rhs and boundary condition information
01057                   loc.template assemble(*it); // assemble local stiffness matrix
01058 
01059                   // assemble constraints in hanging nodes
01060                   // as a by-product determine the interpolation factors WITH RESPECT TO NODES OF FATHER
01061                   bool hasHangingNodes = false;
01062                   for (int k=0; k<sfs.size(); k++) // loop over rows, i.e. test functions
01063                         {
01064                           // process only hanging nodes
01065                           if (!this->hanging[l2g[k]]) continue;
01066                           hasHangingNodes = true;
01067 
01068                           // determine position of hanging node in father
01069                           EEntityPointer father=it->father(); // the father element
01070                           GeometryType gtf = father->type(); // fathers type
01071                           assert(gtf==gt); // in hanging node refinement the element type is preserved
01072                           const FieldVector<DT,n>& cpos=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[k].position();
01073                           FieldVector<DT,n> pos = it->geometryInFather().global(cpos); // map corner to father element
01074 
01075                           // evaluate interpolation factors and local to global mapping for father
01076                           for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1).size(); ++i)
01077                                 {
01078                                   alpha[i][k] = Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1)[i].evaluateFunction(0,pos);
01079                                   fl2g[i] = this->vertexmapper.template map<n>(*father,i);
01080                                 }                                                 
01081 
01082                           // assemble the constraint row once
01083                           if (!treated[l2g[k]])
01084                                 {
01085                                   bool throwflag=false;
01086                                   int cnt=0;
01087                                   for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1).size(); ++i)
01088                                         if (std::abs(alpha[i][k])>1E-4)
01089                                           {
01090                                                 if (this->hanging[fl2g[i]])
01091                                                   {
01092                                                       std::cout << "X i=" << father->geometry().corner(i)
01093                                                                 << " k=" << it->geometry().corner(k)
01094                                                                           << " alpha=" << alpha[i][k]
01095                                                                           << " cnt=" << ++cnt
01096                                                                           << std::endl;
01097                                                         throwflag = true;
01098                                                   }
01099                                                 this->A[l2g[k]][fl2g[i]] = MFieldType(); // Note: Interpolation is done a posteriori
01100                                           }
01101                                   if (throwflag)
01102                                         DUNE_THROW(GridError,"hanging node interpolated from hanging node");
01103                                   for (int comp=0; comp<m; comp++)
01104                                         this->A[l2g[k]][l2g[k]][comp][comp] = static_cast<MFieldType>(1);
01105                                   (*f)[l2g[k]] = 0;
01106                                   (*u)[l2g[k]] = 0;
01107                                   treated[l2g[k]] = true;
01108                                 }
01109                         }
01110 
01111                   // accumulate local matrix into global matrix for non-hanging nodes
01112                   for (int i=0; i<sfs.size(); i++) // loop over rows, i.e. test functions
01113                         {
01114                           // process only non-hanging nodes
01115                           if (this->hanging[l2g[i]]) continue;
01116 
01117                           // accumulate matrix
01118                           for (int j=0; j<sfs.size(); j++)
01119                                 {
01120                                   // process only non-hanging nodes
01121                                   if (this->hanging[l2g[j]]) continue;
01122 
01123                                   // the standard entry
01124                                   this->A[l2g[i]][l2g[j]] += loc.mat(i,j);
01125                                 }
01126                           
01127                           // essential boundary condition and rhs
01128                           for (int comp=0; comp<m; comp++)
01129                                 {
01130                                   if (loc.bc(i)[comp]>essential[l2g[i]][comp])
01131                                         {
01132                                           essential[l2g[i]][comp] = loc.bc(i)[comp];
01133                                           (*f)[l2g[i]][comp] = loc.rhs(i)[comp];
01134                                         }
01135                                   if (essential[l2g[i]][comp]==BoundaryConditions::neumann)
01136                                         (*f)[l2g[i]][comp] += loc.rhs(i)[comp];
01137                                 }
01138                         }
01139 
01140                   // add corrections for hanging nodes
01141                   if (hasHangingNodes)
01142                         {
01143                           // a map that inverts l2g
01144                           std::map<int,int> g2l;
01145                           for (int i=0; i<sfs.size(); i++)
01146                                 g2l[l2g[i]] = i;
01147 
01148                           // a map that inverts fl2g
01149                           std::map<int,int> fg2l;
01150                           for (int i=0; i<sfs.size(); i++)
01151                                 fg2l[fl2g[i]] = i;
01152 
01153                           EEntityPointer father=it->father(); // DEBUG
01154 
01155                           // loop over nodes (rows) in father, assume father has same geometry type
01156                           for (int i=0; i<sfs.size(); ++i)
01157                                 {
01158                                   // corrections to matrix
01159                                   // loop over all nodes (columns) in father, assume father has same geometry type
01160                                   for (int j=0; j<sfs.size(); ++j)
01161                                         {
01162                                           // loop over hanging nodes
01163                                           for (int k=0; k<sfs.size(); k++)
01164                                                 {
01165                                                   // process only hanging nodes
01166                                                   if (!this->hanging[l2g[k]]) continue;
01167 
01168                                                   // first term, with i a coarse vertex
01169                                                   if ( std::abs(alpha[j][k])>1E-4 && g2l.find(fl2g[i])!=g2l.end() )
01170                                                         {
01171                                                           accumulate(this->A[fl2g[i]][fl2g[j]],alpha[j][k],loc.mat(g2l[fl2g[i]],k));
01172                                                         }
01173 
01174                                                   // first term, with i a fine non-hanging vertex
01175                                                   if ( std::abs(alpha[j][k])>1E-4 && fg2l.find(l2g[i])==fg2l.end() && !this->hanging[l2g[i]])
01176                                                         {
01177                                                           accumulate(this->A[l2g[i]][fl2g[j]],alpha[j][k],loc.mat(i,k));
01178                                                         }
01179 
01180                                                   // second term, with j a coarse vertex
01181                                                   if ( std::abs(alpha[i][k])>1E-4 && g2l.find(fl2g[j])!=g2l.end() )
01182                                                         {
01183                                                           accumulate(this->A[fl2g[i]][fl2g[j]],alpha[i][k],loc.mat(k,g2l[fl2g[j]]));
01184                                                         }
01185 
01186                                                   // second term, with j a fine non-hanging vertex
01187                                                   if ( std::abs(alpha[i][k])>1E-4 && fg2l.find(l2g[j])==fg2l.end() && !this->hanging[l2g[j]])
01188                                                         {
01189                                                           accumulate(this->A[fl2g[i]][l2g[j]],alpha[i][k],loc.mat(k,j));
01190                                                         }
01191 
01192                                                   // third term, loop over hanging nodes
01193                                                   for (int l=0; l<sfs.size(); l++)
01194                                                         {
01195                                                           // process only hanging nodes
01196                                                           if (!this->hanging[l2g[l]]) continue;
01197 
01198                                                           if ( std::abs(alpha[i][k])>1E-4 && std::abs(alpha[j][l])>1E-4 )
01199                                                                 {
01200                                                                   accumulate(this->A[fl2g[i]][fl2g[j]],alpha[i][k]*alpha[j][l],loc.mat(k,l));
01201                                                                 }
01202                                                         }
01203                                                 }
01204                                         }
01205 
01206                                   // corrections to rhs
01207                                   for (int comp=0; comp<m; comp++)
01208                                         if (essential[fl2g[i]][comp]==BoundaryConditions::neumann)
01209                                           for (int k=0; k<sfs.size(); k++)
01210                                                 {
01211                                                   // process only hanging nodes
01212                                                   if (!this->hanging[l2g[k]]) continue;
01213                                                   
01214                                                   if ( std::abs(alpha[i][k])>1E-4 && loc.bc(k)[comp]==BoundaryConditions::neumann )
01215                                                         (*f)[fl2g[i]][comp] += alpha[i][k]*loc.rhs(k)[comp];
01216                                                 }
01217                                 }
01218                         }
01219                 }
01220 
01221           // accumulate matrix entries, should do also for systems ...
01222           if (this->extendOverlap)
01223                 this->sumEntries();
01224 
01225           // send around boundary conditions
01226           if (this->extendOverlap)
01227                 {
01228                   AccumulateBCFlags datahandle(this->grid,this->vertexmapper,essential,*f);
01229                   this->lc.template communicate<AccumulateBCFlags>(datahandle,InteriorBorder_InteriorBorder_Interface,ForwardCommunication);
01230                   // on a nonvoerlapping grid we have the correct BC at all interior and border vertices
01231                 }
01232           // what about the overlapping case ?
01233 
01234           // muck up ghost vertices
01235           VIterator vendit = this->gv.template end<n>();
01236           for (VIterator it = this->gv.template begin<n>(); it!=vendit; ++it)
01237                 if (it->partitionType()==GhostEntity)
01238                   {
01239                         // index of this vertex
01240                         int i=this->vertexmapper.map(*it);
01241 
01242                         // muck up row
01243                         coliterator endj=this->A[i].end();
01244                         for (coliterator j=this->A[i].begin(); j!=endj; ++j)
01245                           {
01246                                 (*j) = MFieldType();
01247                                 if (j.index()==i)
01248                                   for (int comp=0; comp<m; comp++)
01249                                         (*j)[comp][comp] = static_cast<MFieldType>(1);
01250                           }
01251                         (*f)[i] = 0;
01252                   }
01253 
01254           // put in essential boundary conditions
01255           rowiterator endi=this->A.end();
01256           for (rowiterator i=this->A.begin(); i!=endi; ++i)
01257                 {
01258                   // muck up extra rows
01259                   if (i.index()>=this->vertexmapper.size())
01260                         {
01261                           coliterator endj=(*i).end();
01262                           for (coliterator j=(*i).begin(); j!=endj; ++j)
01263                                 {
01264                                   (*j) = MFieldType();
01265                                   if (j.index()==i.index())
01266                                         for (int comp=0; comp<m; comp++)
01267                                           (*j)[comp][comp] = static_cast<MFieldType>(1);
01268                                 }
01269                           (*f)[i.index()] = 0;
01270                           continue;
01271                         }
01272 
01273                   // insert dirichlet ans processor boundary conditions
01274                   if (!this->hanging[i.index()])
01275                         {
01276                           for (int icomp=0; icomp<m; icomp++)
01277                                 if (essential[i.index()][icomp]!=BoundaryConditions::neumann)
01278                                   {
01279                                         coliterator endj=(*i).end();
01280                                         for (coliterator j=(*i).begin(); j!=endj; ++j)
01281                                           if (j.index()==i.index())
01282                                                 {
01283                                                   for (int jcomp=0; jcomp<m; jcomp++)
01284                                                         if (icomp==jcomp)
01285                                                           (*j)[icomp][jcomp] = static_cast<MFieldType>(1);
01286                                                         else
01287                                                           (*j)[icomp][jcomp] = MFieldType();                                                      
01288                                                 }
01289                                           else
01290                                                 {
01291                                                   for (int jcomp=0; jcomp<m; jcomp++)
01292                                                         (*j)[icomp][jcomp] = MFieldType();
01293                                                 }
01294                                         (*u)[i.index()][icomp] = (*f)[i.index()][icomp];
01295                                   }
01296                         }
01297                 }
01298         } 
01299 
01301         void interpolateHangingNodes (P1Function<GV,RT,LC,m>& u)
01302         {
01303           // allocate flag vector to note hanging nodes whose row has been assembled
01304           std::vector<unsigned char> treated(this->vertexmapper.size());
01305           for (std::size_t i=0; i<treated.size(); i++) treated[i] = false;
01306 
01307           // run over all leaf elements
01308           Iterator eendit = this->gv.template end<0>();
01309           for (Iterator it = this->gv.template begin<0>(); it!=eendit; ++it)
01310                 {
01311                   // get access to shape functions for P1 elements
01312                   Dune::GeometryType gt = it->type();
01313                   const