/[escript]/branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3698 by caltinay, Tue Nov 29 00:47:29 2011 UTC revision 3699 by caltinay, Thu Dec 1 22:59:42 2011 UTC
# Line 240  bool Rectangle::ownSample(int fsCode, in Line 240  bool Rectangle::ownSample(int fsCode, in
240  {  {
241  #ifdef ESYS_MPI  #ifdef ESYS_MPI
242      if (fsCode == Nodes) {      if (fsCode == Nodes) {
243          const index_t myFirst=getNumNodes()*m_mpiInfo->rank;          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
244          const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
245          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
246      } else      } else
247          throw RipleyException("ownSample() only implemented for Nodes");          throw RipleyException("ownSample() only implemented for Nodes");
# Line 269  void Rectangle::interpolateOnDomain(escr Line 269  void Rectangle::interpolateOnDomain(escr
269          case Nodes:          case Nodes:
270          case DegreesOfFreedom:          case DegreesOfFreedom:
271              switch (target.getFunctionSpace().getTypeCode()) {              switch (target.getFunctionSpace().getTypeCode()) {
272                    case DegreesOfFreedom:
273                        target=in;
274                        break;
275    
276                  case Elements:                  case Elements:
277                  {                  {
278                      const double tmp0_2 = 0.62200846792814621559;                      const double tmp0_2 = 0.62200846792814621559;
# Line 295  void Rectangle::interpolateOnDomain(escr Line 299  void Rectangle::interpolateOnDomain(escr
299                      break;                      break;
300                  }                  }
301    
                 case DegreesOfFreedom:  
                     target=in;  
                     break;  
                       
302                  default:                  default:
303                      throw RipleyException(msg.str());                      throw RipleyException(msg.str());
304              }              }
305                break;
306          default:          default:
307              throw RipleyException(msg.str());              throw RipleyException(msg.str());
308      }      }
# Line 313  Paso_SystemMatrixPattern* Rectangle::get Line 314  Paso_SystemMatrixPattern* Rectangle::get
314      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
315          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
316    
317  /*      // connector
318      // create distribution      RankVector neighbour;
319      IndexVector dist;      IndexVector offsetInShared(1,0);
320      for (index_t i=0; i<m_mpiInfo->size+1; i++)      IndexVector sendShared, recvShared;
321          dist.push_back(i*getNumNodes());      const IndexVector faces=getNumFacesPerBoundary();
322      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      const index_t left = (m_offset0==0 ? 0 : 1);
323              &dist[0], 1, 0);      const index_t bottom = (m_offset1==0 ? 0 : 1);
324        // corner node from bottom-left
325      // connectors      if (faces[0] == 0 && faces[2] == 0) {
326      dim_t numNeighbours = 0;          neighbour.push_back(m_mpiInfo->rank-m_NX-1);
327      RankVector neighbour(numNeighbours);          offsetInShared.push_back(offsetInShared.back()+1);
328      IndexVector offsetInShared(numNeighbours+1);          sendShared.push_back(m_nodeId[m_N0+left]);
329      IndexVector shared(offsetInShared[numNeighbours]);          recvShared.push_back(m_nodeId[0]);
330        }
331        // bottom edge
332        if (faces[2] == 0) {
333            neighbour.push_back(m_mpiInfo->rank-m_NX);
334            offsetInShared.push_back(offsetInShared.back()+m_N0-left);
335            for (dim_t i=left; i<m_N0; i++) {
336                // easy case, we know the neighbour id's
337                sendShared.push_back(m_nodeId[m_N0+i]);
338                recvShared.push_back(m_nodeId[i]);
339            }
340        }
341        // corner node from bottom-right
342        if (faces[1] == 0 && faces[2] == 0) {
343            neighbour.push_back(m_mpiInfo->rank-m_NX+1);
344            const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
345            const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);
346            const index_t first=m_nodeDistribution[neighbour.back()];
347            offsetInShared.push_back(offsetInShared.back()+1);
348            sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);
349            recvShared.push_back(first+N0*(N1-1));
350        }
351        // left edge
352        if (faces[0] == 0) {
353            neighbour.push_back(m_mpiInfo->rank-1);
354            offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
355            for (dim_t i=bottom; i<m_N1; i++) {
356                // easy case, we know the neighbour id's
357                sendShared.push_back(m_nodeId[i*m_N0+1]);
358                recvShared.push_back(m_nodeId[i*m_N0]);
359            }
360        }
361        // right edge
362        if (faces[1] == 0) {
363            neighbour.push_back(m_mpiInfo->rank+1);
364            const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
365            const index_t first=m_nodeDistribution[neighbour.back()];
366            offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
367            for (dim_t i=bottom; i<m_N1; i++) {
368                sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);
369                recvShared.push_back(first+rightN0*(i-bottom));
370            }
371        }
372        // corner node from top-left
373        if (faces[0] == 0 && faces[3] == 0) {
374            neighbour.push_back(m_mpiInfo->rank+m_NX-1);
375            const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
376            const index_t first=m_nodeDistribution[neighbour.back()];
377            offsetInShared.push_back(offsetInShared.back()+1);
378            sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);
379            recvShared.push_back(first+N0-1);
380        }
381        // top edge
382        if (faces[3] == 0) {
383            neighbour.push_back(m_mpiInfo->rank+m_NX);
384            const index_t first=m_nodeDistribution[neighbour.back()];
385            offsetInShared.push_back(offsetInShared.back()+m_N0-left);
386            for (dim_t i=left; i<m_N0; i++) {
387                sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);
388                recvShared.push_back(first+i-left);
389            }
390        }
391        // corner node from top-right
392        if (faces[1] == 0 && faces[3] == 0) {
393            neighbour.push_back(m_mpiInfo->rank+m_NX+1);
394            const index_t first=m_nodeDistribution[neighbour.back()];
395            offsetInShared.push_back(offsetInShared.back()+1);
396            sendShared.push_back(m_nodeId[m_N0*m_N1-1]);
397            recvShared.push_back(first);
398        }
399        const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];
400        cout << "--- rcv_shcomp ---" << endl;
401        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
402        for (size_t i=0; i<neighbour.size(); i++) {
403            cout << "neighbor[" << i << "]=" << neighbour[i]
404                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
405        }
406        for (size_t i=0; i<recvShared.size(); i++) {
407            cout << "shared[" << i << "]=" << recvShared[i] << endl;
408        }
409        cout << "--- snd_shcomp ---" << endl;
410        for (size_t i=0; i<sendShared.size(); i++) {
411            cout << "shared[" << i << "]=" << sendShared[i] << endl;
412        }
413    
414      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
415              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
416              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
417      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
418              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
419              &offsetInShared[0], 1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
420      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
421        Paso_SharedComponents_free(snd_shcomp);
422        Paso_SharedComponents_free(rcv_shcomp);
423    
424      // create patterns      // create patterns
425      dim_t M=0, N=0;      dim_t M, N;
426      int* ptr=NULL;      IndexVector ptr(1,0);
427      index_t* index=NULL;      IndexVector index;
428    
429        // main pattern
430        for (index_t i=0; i<numDOF; i++) {
431            // always add the node itself
432            index.push_back(i);
433            int num=insertNeighbours(index, i);
434            ptr.push_back(ptr.back()+num+1);
435        }
436        M=N=ptr.size()-1;
437        // paso will manage the memory
438        index_t* indexC = MEMALLOC(index.size(),index_t);
439        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
440        copy(index.begin(), index.end(), indexC);
441        copy(ptr.begin(), ptr.end(), ptrC);
442      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
443              M, N, ptr, index);              M, N, ptrC, indexC);
444      Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
445              M, N, ptr, index);      cout << "--- main_pattern ---" << endl;
446      Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      cout << "M=" << M << ", N=" << N << endl;
447              M, N, ptr, index);      for (size_t i=0; i<ptr.size(); i++) {
448            cout << "ptr[" << i << "]=" << ptr[i] << endl;
449        }
450        for (size_t i=0; i<index.size(); i++) {
451            cout << "index[" << i << "]=" << index[i] << endl;
452        }
453    
454        ptr.clear();
455        index.clear();
456    
457        // column & row couple patterns
458        Paso_Pattern *colCouplePattern, *rowCouplePattern;
459        generateCouplePatterns(&colCouplePattern, &rowCouplePattern);
460    
461        // allocate paso distribution
462        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
463                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
464    
465      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
466              MATRIX_FORMAT_DEFAULT, distribution, distribution,              MATRIX_FORMAT_DEFAULT, distribution, distribution,
# Line 354  Paso_SystemMatrixPattern* Rectangle::get Line 470  Paso_SystemMatrixPattern* Rectangle::get
470      Paso_Pattern_free(colCouplePattern);      Paso_Pattern_free(colCouplePattern);
471      Paso_Pattern_free(rowCouplePattern);      Paso_Pattern_free(rowCouplePattern);
472      Paso_Distribution_free(distribution);      Paso_Distribution_free(distribution);
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
473      return pattern;      return pattern;
 */  
     throw RipleyException("getPattern() not implemented");  
474  }  }
475    
476  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 425  pair<double,double> Rectangle::getFirstC Line 537  pair<double,double> Rectangle::getFirstC
537  //protected  //protected
538  dim_t Rectangle::getNumFaceElements() const  dim_t Rectangle::getNumFaceElements() const
539  {  {
540        const IndexVector faces = getNumFacesPerBoundary();
541      dim_t n=0;      dim_t n=0;
542      //left      for (size_t i=0; i<faces.size(); i++)
543      if (m_offset0==0)          n+=faces[i];
         n+=m_NE1;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         n+=m_NE1;  
     //bottom  
     if (m_offset1==0)  
         n+=m_NE0;  
     //top  
     if (m_mpiInfo->rank/m_NX==m_NY-1)  
         n+=m_NE0;  
   
544      return n;      return n;
545  }  }
546    
# Line 518  void Rectangle::populateSampleIds() Line 620  void Rectangle::populateSampleIds()
620      }      }
621      if (left>0 && bottom>0) {      if (left>0 && bottom>0) {
622          const int neighbour=m_mpiInfo->rank-m_NX-1;          const int neighbour=m_mpiInfo->rank-m_NX-1;
623          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
624          const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);          const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
625          m_nodeId[0]=m_nodeDistribution[neighbour]+bottomN1*bottomN0-1;          m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;
626      }      }
627    
628      // the rest of the id's are contiguous      // the rest of the id's are contiguous
# Line 547  void Rectangle::populateSampleIds() Line 649  void Rectangle::populateSampleIds()
649      }      }
650  }  }
651    
652    //private
653    int Rectangle::insertNeighbours(IndexVector& index, index_t node) const
654    {
655        const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
656        const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
657        const int x=node%myN0;
658        const int y=node/myN0;
659        int num=0;
660        if (y>0) {
661            if (x>0) {
662                // bottom-left
663                index.push_back(node-myN0-1);
664                num++;
665            }
666            // bottom
667            index.push_back(node-myN0);
668            num++;
669            if (x<myN0-1) {
670                // bottom-right
671                index.push_back(node-myN0+1);
672                num++;
673            }
674        }
675        if (x<myN0-1) {
676            // right
677            index.push_back(node+1);
678            num++;
679            if (y<myN1-1) {
680                // top-right
681                index.push_back(node+myN0+1);
682                num++;
683            }
684        }
685        if (y<myN1-1) {
686            // top
687            index.push_back(node+myN0);
688            num++;
689            if (x>0) {
690                // top-left
691                index.push_back(node+myN0-1);
692                num++;
693            }
694        }
695        if (x>0) {
696            // left
697            index.push_back(node-1);
698            num++;
699        }
700    
701        return num;
702    }
703    
704    //private
705    void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
706    {
707        IndexVector ptr(1,0);
708        IndexVector index;
709        const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
710        const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
711        const IndexVector faces=getNumFacesPerBoundary();
712    
713        // bottom edge
714        dim_t n=0;
715        if (faces[0] == 0) {
716            index.push_back(2*(myN0+myN1+1));
717            index.push_back(2*(myN0+myN1+1)+1);
718            n+=2;
719            if (faces[2] == 0) {
720                index.push_back(0);
721                index.push_back(1);
722                index.push_back(2);
723                n+=3;
724            }
725        } else if (faces[2] == 0) {
726            index.push_back(1);
727            index.push_back(2);
728            n+=2;
729        }
730        // n=neighbours of bottom-left corner node
731        ptr.push_back(ptr.back()+n);
732        n=0;
733        if (faces[2] == 0) {
734            for (dim_t i=1; i<myN0-1; i++) {
735                index.push_back(i);
736                index.push_back(i+1);
737                index.push_back(i+2);
738                ptr.push_back(ptr.back()+3);
739            }
740            index.push_back(myN0-1);
741            index.push_back(myN0);
742            n+=2;
743            if (faces[1] == 0) {
744                index.push_back(myN0+1);
745                index.push_back(myN0+2);
746                index.push_back(myN0+3);
747                n+=3;
748            }
749        } else {
750            for (dim_t i=1; i<myN0-1; i++) {
751                ptr.push_back(ptr.back());
752            }
753            if (faces[1] == 0) {
754                index.push_back(myN0+2);
755                index.push_back(myN0+3);
756                n+=2;
757            }
758        }
759        // n=neighbours of bottom-right corner node
760        ptr.push_back(ptr.back()+n);
761    
762        // 2nd row to 2nd last row
763        for (dim_t i=1; i<myN1-1; i++) {
764            // left edge
765            if (faces[0] == 0) {
766                index.push_back(2*(myN0+myN1+2)-i);
767                index.push_back(2*(myN0+myN1+2)-i-1);
768                index.push_back(2*(myN0+myN1+2)-i-2);
769                ptr.push_back(ptr.back()+3);
770            } else {
771                ptr.push_back(ptr.back());
772            }
773            for (dim_t j=1; j<myN0-1; j++) {
774                ptr.push_back(ptr.back());
775            }
776            // right edge
777            if (faces[1] == 0) {
778                index.push_back(myN0+i+1);
779                index.push_back(myN0+i+2);
780                index.push_back(myN0+i+3);
781                ptr.push_back(ptr.back()+3);
782            } else {
783                ptr.push_back(ptr.back());
784            }
785        }
786    
787        // top edge
788        n=0;
789        if (faces[0] == 0) {
790            index.push_back(2*myN0+myN1+5);
791            index.push_back(2*myN0+myN1+4);
792            n+=2;
793            if (faces[3] == 0) {
794                index.push_back(2*myN0+myN1+3);
795                index.push_back(2*myN0+myN1+2);
796                index.push_back(2*myN0+myN1+1);
797                n+=3;
798            }
799        } else if (faces[3] == 0) {
800            index.push_back(2*myN0+myN1+2);
801            index.push_back(2*myN0+myN1+1);
802            n+=2;
803        }
804        // n=neighbours of top-left corner node
805        ptr.push_back(ptr.back()+n);
806        n=0;
807        if (faces[3] == 0) {
808            for (dim_t i=1; i<myN0-1; i++) {
809                index.push_back(2*myN0+myN1+i+1);
810                index.push_back(2*myN0+myN1+i);
811                index.push_back(2*myN0+myN1+i-1);
812                ptr.push_back(ptr.back()+3);
813            }
814            index.push_back(myN0+myN1+4);
815            index.push_back(myN0+myN1+3);
816            n+=2;
817            if (faces[1] == 0) {
818                index.push_back(myN0+myN1+2);
819                index.push_back(myN0+myN1+1);
820                index.push_back(myN0+myN1);
821                n+=3;
822            }
823        } else {
824            for (dim_t i=1; i<myN0-1; i++) {
825                ptr.push_back(ptr.back());
826            }
827            if (faces[1] == 0) {
828                index.push_back(myN0+myN1+1);
829                index.push_back(myN0+myN1);
830                n+=2;
831            }
832        }
833        // n=neighbours of top-right corner node
834        ptr.push_back(ptr.back()+n);
835    
836        dim_t M=ptr.size()-1;
837        map<index_t,index_t> idMap;
838        dim_t N=0;
839        for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {
840            if (idMap.find(*it)==idMap.end()) {
841                idMap[*it]=N;
842                *it=N++;
843            } else {
844                *it=idMap[*it];
845            }
846        }
847    
848        cout << "--- colCouple_pattern ---" << endl;
849        cout << "M=" << M << ", N=" << N << endl;
850        for (size_t i=0; i<ptr.size(); i++) {
851            cout << "ptr[" << i << "]=" << ptr[i] << endl;
852        }
853        for (size_t i=0; i<index.size(); i++) {
854            cout << "index[" << i << "]=" << index[i] << endl;
855        }
856    
857        // now build the row couple pattern
858        IndexVector ptr2(1,0);
859        IndexVector index2;
860        for (dim_t id=0; id<N; id++) {
861            n=0;
862            for (dim_t i=0; i<M; i++) {
863                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
864                    if (index[j]==id) {
865                        index2.push_back(i);
866                        n++;
867                        break;
868                    }
869                }
870            }
871            ptr2.push_back(ptr2.back()+n);
872        }
873    
874        cout << "--- rowCouple_pattern ---" << endl;
875        cout << "M=" << N << ", N=" << M << endl;
876        for (size_t i=0; i<ptr2.size(); i++) {
877            cout << "ptr[" << i << "]=" << ptr2[i] << endl;
878        }
879        for (size_t i=0; i<index2.size(); i++) {
880            cout << "index[" << i << "]=" << index2[i] << endl;
881        }
882    
883        // paso will manage the memory
884        index_t* indexC = MEMALLOC(index.size(), index_t);
885        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
886        copy(index.begin(), index.end(), indexC);
887        copy(ptr.begin(), ptr.end(), ptrC);
888        *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
889    
890        // paso will manage the memory
891        indexC = MEMALLOC(index2.size(), index_t);
892        ptrC = MEMALLOC(ptr2.size(), index_t);
893        copy(index2.begin(), index2.end(), indexC);
894        copy(ptr2.begin(), ptr2.end(), ptrC);
895        *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);
896    }
897    
898  } // end of namespace ripley  } // end of namespace ripley
899    
Legend:
Removed from v.3698  
changed lines
  Added in v.3699
  ViewVC Help
Powered by ViewVC 1.1.26