/[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 3752 by caltinay, Tue Jan 3 08:06:36 2012 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 63  Rectangle::Rectangle(int n0, int n1, dou Line 63  Rectangle::Rectangle(int n0, int n1, dou
63    
64      // bottom-left node is at (offset0,offset1) in global mesh      // bottom-left node is at (offset0,offset1) in global mesh
65      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
66      if (m_mpiInfo->rank%m_NX>0)      if (m_offset0 > 0)
67          m_offset0--;          m_offset0--;
68      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);      m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
69      if (m_mpiInfo->rank/m_NX>0)      if (m_offset1 > 0)
70          m_offset1--;          m_offset1--;
71    
72      populateSampleIds();      populateSampleIds();
73        createPattern();
74  }  }
75    
76  Rectangle::~Rectangle()  Rectangle::~Rectangle()
# Line 265  bool Rectangle::ownSample(int fsCode, in Line 266  bool Rectangle::ownSample(int fsCode, in
266  {  {
267  #ifdef ESYS_MPI  #ifdef ESYS_MPI
268      if (fsCode == Nodes) {      if (fsCode == Nodes) {
269          const index_t left = (m_offset0==0 ? 0 : 1);          return (m_dofMap[id] < getNumDOF());
         const index_t bottom = (m_offset1==0 ? 0 : 1);  
         const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);  
         const index_t top = (m_mpiInfo->rank/m_NX==m_NY-1 ? m_N1 : m_N1-1);  
         const index_t x=id%m_N0;  
         const index_t y=id/m_N0;  
         return (x>=left && x<right && y>=bottom && y<top);  
270      } else {      } else {
271          stringstream msg;          stringstream msg;
272          msg << "ownSample() not implemented for "          msg << "ownSample() not implemented for "
# Line 755  Paso_SystemMatrixPattern* Rectangle::get Line 750  Paso_SystemMatrixPattern* Rectangle::get
750      if (reducedRowOrder || reducedColOrder)      if (reducedRowOrder || reducedColOrder)
751          throw RipleyException("getPattern() not implemented for reduced order");          throw RipleyException("getPattern() not implemented for reduced order");
752    
753      // connector      return m_pattern;
     RankVector neighbour;  
     IndexVector offsetInShared(1,0);  
     IndexVector sendShared, recvShared;  
     const IndexVector faces=getNumFacesPerBoundary();  
     const index_t nDOF0 = (m_gNE0+1)/m_NX;  
     const index_t nDOF1 = (m_gNE1+1)/m_NY;  
     const int numDOF=nDOF0*nDOF1;  
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     vector<IndexVector> colIndices(numDOF); // for the couple blocks  
     int numShared=0;  
   
     m_dofMap.assign(getNumNodes(), 0);  
 #pragma omp parallel for  
     for (index_t i=bottom; i<m_N1; i++) {  
         for (index_t j=left; j<m_N0; j++) {  
             m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;  
         }  
     }  
   
     // corner node from bottom-left  
     if (faces[0] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX-1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(0);  
         recvShared.push_back(numDOF+numShared);  
         colIndices[0].push_back(numShared);  
         m_dofMap[0]=numDOF+numShared;  
         ++numShared;  
     }  
     // bottom edge  
     if (faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX);  
         offsetInShared.push_back(offsetInShared.back()+nDOF0);  
         for (dim_t i=0; i<nDOF0; i++, numShared++) {  
             sendShared.push_back(i);  
             recvShared.push_back(numDOF+numShared);  
             if (i>0)  
                 colIndices[i-1].push_back(numShared);  
             colIndices[i].push_back(numShared);  
             if (i<nDOF0-1)  
                 colIndices[i+1].push_back(numShared);  
             m_dofMap[i+left]=numDOF+numShared;  
         }  
     }  
     // corner node from bottom-right  
     if (faces[1] == 0 && faces[2] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-m_NX+1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(nDOF0-1);  
         recvShared.push_back(numDOF+numShared);  
         colIndices[nDOF0-2].push_back(numShared);  
         colIndices[nDOF0-1].push_back(numShared);  
         m_dofMap[m_N0-1]=numDOF+numShared;  
         ++numShared;  
     }  
     // right edge  
     if (faces[1] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+1);  
         offsetInShared.push_back(offsetInShared.back()+nDOF1);  
         for (dim_t i=0; i<nDOF1; i++, numShared++) {  
             sendShared.push_back((i+1)*(nDOF0)-1);  
             recvShared.push_back(numDOF+numShared);  
             if (i>0)  
                 colIndices[i*(nDOF0)-1].push_back(numShared);  
             colIndices[(i+1)*(nDOF0)-1].push_back(numShared);  
             if (i<nDOF1-1)  
                 colIndices[(i+2)*(nDOF0)-1].push_back(numShared);  
             m_dofMap[(i+bottom+1)*m_N0-1]=numDOF+numShared;  
         }  
     }  
     // corner node from top-right  
     if (faces[1] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX+1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(numDOF-1);  
         recvShared.push_back(numDOF+numShared);  
         colIndices[numDOF-1].push_back(numShared);  
         m_dofMap[m_N0*m_N1-1]=numDOF+numShared;  
         ++numShared;  
     }  
     // top edge  
     if (faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX);  
         offsetInShared.push_back(offsetInShared.back()+nDOF0);  
         for (dim_t i=0; i<nDOF0; i++, numShared++) {  
             sendShared.push_back(numDOF-nDOF0+i);  
             recvShared.push_back(numDOF+numShared);  
             if (i>0)  
                 colIndices[numDOF-nDOF0+i-1].push_back(numShared);  
             colIndices[numDOF-nDOF0+i].push_back(numShared);  
             if (i<nDOF0-1)  
                 colIndices[numDOF-nDOF0+i+1].push_back(numShared);  
             m_dofMap[m_N0*(m_N1-1)+i+left]=numDOF+numShared;  
         }  
     }  
     // corner node from top-left  
     if (faces[0] == 0 && faces[3] == 0) {  
         neighbour.push_back(m_mpiInfo->rank+m_NX-1);  
         offsetInShared.push_back(offsetInShared.back()+1);  
         sendShared.push_back(numDOF-nDOF0);  
         recvShared.push_back(numDOF+numShared);  
         colIndices[numDOF-nDOF0].push_back(numShared);  
         m_dofMap[m_N0*(m_N1-1)]=numDOF+numShared;  
         ++numShared;  
     }  
     // left edge  
     if (faces[0] == 0) {  
         neighbour.push_back(m_mpiInfo->rank-1);  
         offsetInShared.push_back(offsetInShared.back()+nDOF1);  
         for (dim_t i=0; i<nDOF1; i++, numShared++) {  
             sendShared.push_back(i*nDOF0);  
             recvShared.push_back(numDOF+numShared);  
             if (i>0)  
                 colIndices[(i-1)*nDOF0].push_back(numShared);  
             colIndices[i*nDOF0].push_back(numShared);  
             if (i<nDOF1-1)  
                 colIndices[(i+1)*nDOF0].push_back(numShared);  
             m_dofMap[(i+bottom)*m_N0]=numDOF+numShared;  
         }  
     }  
   
     /*  
     cout << "--- rcv_shcomp ---" << endl;  
     cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;  
     for (size_t i=0; i<neighbour.size(); i++) {  
         cout << "neighbor[" << i << "]=" << neighbour[i]  
             << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;  
     }  
     for (size_t i=0; i<recvShared.size(); i++) {  
         cout << "shared[" << i << "]=" << recvShared[i] << endl;  
     }  
     cout << "--- snd_shcomp ---" << endl;  
     for (size_t i=0; i<sendShared.size(); i++) {  
         cout << "shared[" << i << "]=" << sendShared[i] << endl;  
     }  
     */  
   
     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &sendShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(  
             numDOF, neighbour.size(), &neighbour[0], &recvShared[0],  
             &offsetInShared[0], 1, 0, m_mpiInfo);  
     Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);  
     Paso_SharedComponents_free(snd_shcomp);  
     Paso_SharedComponents_free(rcv_shcomp);  
   
     // create patterns  
     dim_t M, N;  
     IndexVector ptr(1,0);  
     IndexVector index;  
   
     // main pattern  
     for (index_t i=0; i<numDOF; i++) {  
         // always add the node itself  
         index.push_back(i);  
         const int num=insertNeighbours(index, i);  
         ptr.push_back(ptr.back()+num+1);  
     }  
     M=N=ptr.size()-1;  
     // paso will manage the memory  
     index_t* indexC = MEMALLOC(index.size(),index_t);  
     index_t* ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
     /*  
     cout << "--- main_pattern ---" << endl;  
     cout << "M=" << M << ", N=" << N << endl;  
     for (size_t i=0; i<ptr.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr[i] << endl;  
     }  
     for (size_t i=0; i<index.size(); i++) {  
         cout << "index[" << i << "]=" << index[i] << endl;  
     }  
     */  
   
     // column & row couple patterns  
     ptr.assign(1, 0);  
     index.clear();  
   
     for (index_t i=0; i<numDOF; i++) {  
         index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());  
         ptr.push_back(ptr.back()+colIndices[i].size());  
     }  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index.size(), index_t);  
     ptrC = MEMALLOC(ptr.size(), index_t);  
     copy(index.begin(), index.end(), indexC);  
     copy(ptr.begin(), ptr.end(), ptrC);  
     M=ptr.size()-1;  
     N=numShared;  
     Paso_Pattern *colCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             M, N, ptrC, indexC);  
   
     /*  
     cout << "--- colCouple_pattern ---" << endl;  
     cout << "M=" << M << ", N=" << N << endl;  
     for (size_t i=0; i<ptr.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr[i] << endl;  
     }  
     for (size_t i=0; i<index.size(); i++) {  
         cout << "index[" << i << "]=" << index[i] << endl;  
     }  
     */  
   
     // now build the row couple pattern  
     IndexVector ptr2(1,0);  
     IndexVector index2;  
     for (dim_t id=0; id<N; id++) {  
         dim_t n=0;  
         for (dim_t i=0; i<M; i++) {  
             for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {  
                 if (index[j]==id) {  
                     index2.push_back(i);  
                     n++;  
                     break;  
                 }  
             }  
         }  
         ptr2.push_back(ptr2.back()+n);  
     }  
   
     // paso will manage the memory  
     indexC = MEMALLOC(index2.size(), index_t);  
     ptrC = MEMALLOC(ptr2.size(), index_t);  
     copy(index2.begin(), index2.end(), indexC);  
     copy(ptr2.begin(), ptr2.end(), ptrC);  
     Paso_Pattern *rowCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,  
             N, M, ptrC, indexC);  
   
     /*  
     cout << "--- rowCouple_pattern ---" << endl;  
     cout << "M=" << N << ", N=" << M << endl;  
     for (size_t i=0; i<ptr2.size(); i++) {  
         cout << "ptr[" << i << "]=" << ptr2[i] << endl;  
     }  
     for (size_t i=0; i<index2.size(); i++) {  
         cout << "index[" << i << "]=" << index2[i] << endl;  
     }  
     */  
   
     // allocate paso distribution  
     Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);  
   
     Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, distribution, distribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             connector, connector);  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Paso_Distribution_free(distribution);  
     return pattern;  
754  }  }
755    
756  void Rectangle::Print_Mesh_Info(const bool full) const  void Rectangle::Print_Mesh_Info(const bool full) const
# Line 1117  void Rectangle::assembleCoordinates(escr Line 853  void Rectangle::assembleCoordinates(escr
853      }      }
854  }  }
855    
856    //protected
857    dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
858    {
859        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
860        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
861        const int x=node%nDOF0;
862        const int y=node/nDOF0;
863        dim_t num=0;
864        // loop through potential neighbours and add to index if positions are
865        // within bounds
866        for (int i1=-1; i1<2; i1++) {
867            for (int i0=-1; i0<2; i0++) {
868                // skip node itself
869                if (i0==0 && i1==0)
870                    continue;
871                // location of neighbour node
872                const int nx=x+i0;
873                const int ny=y+i1;
874                if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
875                    index.push_back(ny*nDOF0+nx);
876                    num++;
877                }
878            }
879        }
880    
881        return num;
882    }
883    
884    //protected
885    void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
886    {
887        const dim_t numComp = in.getDataPointSize();
888        out.requireWrite();
889    
890        const index_t left = (m_offset0==0 ? 0 : 1);
891        const index_t bottom = (m_offset1==0 ? 0 : 1);
892        const dim_t nDOF0 = (m_gNE0+1)/m_NX;
893        const dim_t nDOF1 = (m_gNE1+1)/m_NY;
894    #pragma omp parallel for
895        for (index_t i=0; i<nDOF1; i++) {
896            for (index_t j=0; j<nDOF0; j++) {
897                const index_t n=j+left+(i+bottom)*m_N0;
898                const double* src=in.getSampleDataRO(n);
899                copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));
900            }
901        }
902    }
903    
904    //protected
905    void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const
906    {
907        const dim_t numComp = in.getDataPointSize();
908        Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
909        in.requireWrite();
910        Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));
911    
912        const dim_t numDOF = getNumDOF();
913        out.requireWrite();
914        const double* buffer = Paso_Coupler_finishCollect(coupler);
915    
916    #pragma omp parallel for
917        for (index_t i=0; i<getNumNodes(); i++) {
918            const double* src=(m_dofMap[i]<numDOF ?
919                    in.getSampleDataRO(m_dofMap[i])
920                    : &buffer[(m_dofMap[i]-numDOF)*numComp]);
921            copy(src, src+numComp, out.getSampleDataRW(i));
922        }
923    }
924    
925  //private  //private
926  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
927  {  {
# Line 1131  void Rectangle::populateSampleIds() Line 936  void Rectangle::populateSampleIds()
936      }      }
937      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
938      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
939        m_dofId.resize(numDOF);
940        m_elementId.resize(getNumElements());
941        m_faceId.resize(getNumFaceElements());
942    
943  #pragma omp parallel for  #pragma omp parallel
944      for (dim_t i1=0; i1<m_N1; i1++) {      {
945          for (dim_t i0=0; i0<m_N0; i0++) {          // nodes
946              m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;  #pragma omp for nowait
947            for (dim_t i1=0; i1<m_N1; i1++) {
948                for (dim_t i0=0; i0<m_N0; i0++) {
949                    m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
950                }
951          }          }
     }  
952    
953      m_dofId.resize(numDOF);          // degrees of freedom
954      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];  #pragma omp for nowait
955      for (dim_t i=0; i<numDOF; i++)          for (dim_t k=0; k<numDOF; k++)
956          m_dofId[i] = firstId+i;              m_dofId[k] = m_nodeDistribution[m_mpiInfo->rank]+k;
957    
958            // elements
959    #pragma omp for nowait
960            for (dim_t i1=0; i1<m_NE1; i1++) {
961                for (dim_t i0=0; i0<m_NE0; i0++) {
962                    m_elementId[i0+i1*m_NE0]=(m_offset1+i1)*m_gNE0+m_offset0+i0;
963                }
964            }
965    
966            // face elements
967    #pragma omp for
968            for (dim_t k=0; k<getNumFaceElements(); k++)
969                m_faceId[k]=k;
970        } // end parallel section
971    
972      m_nodeTags.assign(getNumNodes(), 0);      m_nodeTags.assign(getNumNodes(), 0);
973      updateTagsInUse(Nodes);      updateTagsInUse(Nodes);
974    
     // elements  
     m_elementId.resize(getNumElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumElements(); k++) {  
         m_elementId[k]=k;  
     }  
975      m_elementTags.assign(getNumElements(), 0);      m_elementTags.assign(getNumElements(), 0);
976      updateTagsInUse(Elements);      updateTagsInUse(Elements);
977    
     // face element id's  
     m_faceId.resize(getNumFaceElements());  
 #pragma omp parallel for  
     for (dim_t k=0; k<getNumFaceElements(); k++) {  
         m_faceId[k]=k;  
     }  
   
978      // generate face offset vector and set face tags      // generate face offset vector and set face tags
979      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
980      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;      const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
# Line 1185  void Rectangle::populateSampleIds() Line 997  void Rectangle::populateSampleIds()
997  }  }
998    
999  //private  //private
1000  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  void Rectangle::createPattern()
1001  {  {
1002      const index_t myN0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1003      const index_t myN1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1004      const int x=node%myN0;      const index_t left = (m_offset0==0 ? 0 : 1);
1005      const int y=node/myN0;      const index_t bottom = (m_offset1==0 ? 0 : 1);
1006      int num=0;  
1007      if (y>0) {      // populate node->DOF mapping with own degrees of freedom.
1008          if (x>0) {      // The rest is assigned in the loop further down
1009              // bottom-left      m_dofMap.assign(getNumNodes(), 0);
1010              index.push_back(node-myN0-1);  #pragma omp parallel for
1011              num++;      for (index_t i=bottom; i<m_N1; i++) {
1012          }          for (index_t j=left; j<m_N0; j++) {
1013          // bottom              m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
         index.push_back(node-myN0);  
         num++;  
         if (x<myN0-1) {  
             // bottom-right  
             index.push_back(node-myN0+1);  
             num++;  
1014          }          }
1015      }      }
1016      if (x<myN0-1) {  
1017          // right      // build list of shared components and neighbours by looping through
1018          index.push_back(node+1);      // all potential neighbouring ranks and checking if positions are
1019          num++;      // within bounds
1020          if (y<myN1-1) {      const dim_t numDOF=nDOF0*nDOF1;
1021              // top-right      vector<IndexVector> colIndices(numDOF); // for the couple blocks
1022              index.push_back(node+myN0+1);      RankVector neighbour;
1023              num++;      IndexVector offsetInShared(1,0);
1024        IndexVector sendShared, recvShared;
1025        int numShared=0;
1026        const int x=m_mpiInfo->rank%m_NX;
1027        const int y=m_mpiInfo->rank/m_NX;
1028        for (int i1=-1; i1<2; i1++) {
1029            for (int i0=-1; i0<2; i0++) {
1030                // skip this rank
1031                if (i0==0 && i1==0)
1032                    continue;
1033                // location of neighbour rank
1034                const int nx=x+i0;
1035                const int ny=y+i1;
1036                if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1037                    neighbour.push_back(ny*m_NX+nx);
1038                    if (i0==0) {
1039                        // sharing top or bottom edge
1040                        const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
1041                        const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
1042                        offsetInShared.push_back(offsetInShared.back()+nDOF0);
1043                        for (dim_t i=0; i<nDOF0; i++, numShared++) {
1044                            sendShared.push_back(firstDOF+i);
1045                            recvShared.push_back(numDOF+numShared);
1046                            if (i>0)
1047                                colIndices[firstDOF+i-1].push_back(numShared);
1048                            colIndices[firstDOF+i].push_back(numShared);
1049                            if (i<nDOF0-1)
1050                                colIndices[firstDOF+i+1].push_back(numShared);
1051                            m_dofMap[firstNode+i]=numDOF+numShared;
1052                        }
1053                    } else if (i1==0) {
1054                        // sharing left or right edge
1055                        const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
1056                        const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
1057                        offsetInShared.push_back(offsetInShared.back()+nDOF1);
1058                        for (dim_t i=0; i<nDOF1; i++, numShared++) {
1059                            sendShared.push_back(firstDOF+i*nDOF0);
1060                            recvShared.push_back(numDOF+numShared);
1061                            if (i>0)
1062                                colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
1063                            colIndices[firstDOF+i*nDOF0].push_back(numShared);
1064                            if (i<nDOF1-1)
1065                                colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
1066                            m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
1067                        }
1068                    } else {
1069                        // sharing a node
1070                        const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
1071                        const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
1072                        offsetInShared.push_back(offsetInShared.back()+1);
1073                        sendShared.push_back(dof);
1074                        recvShared.push_back(numDOF+numShared);
1075                        colIndices[dof].push_back(numShared);
1076                        m_dofMap[node]=numDOF+numShared;
1077                        ++numShared;
1078                    }
1079                }
1080          }          }
1081      }      }
1082      if (y<myN1-1) {  
1083          // top      // create connector
1084          index.push_back(node+myN0);      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1085          num++;              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1086          if (x>0) {              &offsetInShared[0], 1, 0, m_mpiInfo);
1087              // top-left      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1088              index.push_back(node+myN0-1);              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1089              num++;              &offsetInShared[0], 1, 0, m_mpiInfo);
1090          }      m_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1091        Paso_SharedComponents_free(snd_shcomp);
1092        Paso_SharedComponents_free(rcv_shcomp);
1093    
1094        // create main and couple blocks
1095        Paso_Pattern *mainPattern = createMainPattern();
1096        Paso_Pattern *colPattern, *rowPattern;
1097        createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
1098    
1099        // allocate paso distribution
1100        Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1101                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
1102    
1103        // finally create the system matrix
1104        m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1105                distribution, distribution, mainPattern, colPattern, rowPattern,
1106                m_connector, m_connector);
1107    
1108        Paso_Distribution_free(distribution);
1109    
1110        // useful debug output
1111        /*
1112        cout << "--- rcv_shcomp ---" << endl;
1113        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1114        for (size_t i=0; i<neighbour.size(); i++) {
1115            cout << "neighbor[" << i << "]=" << neighbour[i]
1116                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
1117      }      }
1118      if (x>0) {      for (size_t i=0; i<recvShared.size(); i++) {
1119          // left          cout << "shared[" << i << "]=" << recvShared[i] << endl;
         index.push_back(node-1);  
         num++;  
1120      }      }
1121        cout << "--- snd_shcomp ---" << endl;
1122        for (size_t i=0; i<sendShared.size(); i++) {
1123            cout << "shared[" << i << "]=" << sendShared[i] << endl;
1124        }
1125        cout << "--- dofMap ---" << endl;
1126        for (size_t i=0; i<m_dofMap.size(); i++) {
1127            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
1128        }
1129        cout << "--- colIndices ---" << endl;
1130        for (size_t i=0; i<colIndices.size(); i++) {
1131            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
1132        }
1133        */
1134    
1135      return num;      /*
1136        cout << "--- main_pattern ---" << endl;
1137        cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1138        for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1139            cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1140        }
1141        for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1142            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1143        }
1144        */
1145    
1146        /*
1147        cout << "--- colCouple_pattern ---" << endl;
1148        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1149        for (size_t i=0; i<colPattern->numOutput+1; i++) {
1150            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1151        }
1152        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1153            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
1154        }
1155        */
1156    
1157        /*
1158        cout << "--- rowCouple_pattern ---" << endl;
1159        cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
1160        for (size_t i=0; i<rowPattern->numOutput+1; i++) {
1161            cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
1162        }
1163        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
1164            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
1165        }
1166        */
1167    
1168        Paso_Pattern_free(mainPattern);
1169        Paso_Pattern_free(colPattern);
1170        Paso_Pattern_free(rowPattern);
1171  }  }
1172    
1173  //protected  //protected
# Line 1399  void Rectangle::interpolateNodesOnFaces( Line 1333  void Rectangle::interpolateNodesOnFaces(
1333  }  }
1334    
1335  //protected  //protected
 void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const  
 {  
     const dim_t numComp = in.getDataPointSize();  
     out.requireWrite();  
   
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);  
     const index_t top = (m_mpiInfo->rank/m_NX==m_NY-1 ? m_N1 : m_N1-1);  
     index_t n=0;  
     for (index_t i=bottom; i<top; i++) {  
         for (index_t j=left; j<right; j++, n++) {  
             const double* src=in.getSampleDataRO(j+i*m_N0);  
             copy(src, src+numComp, out.getSampleDataRW(n));  
         }  
     }  
 }  
   
 //protected  
1336  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1337          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1338          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
# Line 1426  void Rectangle::assemblePDESingle(Paso_S Line 1341  void Rectangle::assemblePDESingle(Paso_S
1341  {  {
1342      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1343      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
1344      /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */      /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1345      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1346      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1347      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
# Line 1508  void Rectangle::assemblePDESingle(Paso_S Line 1423  void Rectangle::assemblePDESingle(Paso_S
1423      rhs.requireWrite();      rhs.requireWrite();
1424  #pragma omp parallel  #pragma omp parallel
1425      {      {
1426          for (index_t k1_0=0; k1_0<2; k1_0++) { // coloring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1427  #pragma omp for  #pragma omp for
1428              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1429                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE0; ++k0)  {
# Line 1517  void Rectangle::assemblePDESingle(Paso_S Line 1432  void Rectangle::assemblePDESingle(Paso_S
1432                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1433                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1434                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
1435                      /* GENERATOR SNIP_PDE_SINGLE TOP */                      /*** GENERATOR SNIP_PDE_SINGLE TOP */
1436                      ///////////////                      ///////////////
1437                      // process A //                      // process A //
1438                      ///////////////                      ///////////////
# Line 2137  void Rectangle::assemblePDESingle(Paso_S Line 2052  void Rectangle::assemblePDESingle(Paso_S
2052                      if (add_EM_F) {                      if (add_EM_F) {
2053                          //cout << "-- AddtoRHS -- " << endl;                          //cout << "-- AddtoRHS -- " << endl;
2054                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2055                          for (index_t i=0; i<4; i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2056                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2057                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
2058                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
# Line 2156  void Rectangle::assemblePDESingle(Paso_S Line 2071  void Rectangle::assemblePDESingle(Paso_S
2071      } // end of parallel region      } // end of parallel region
2072  }  }
2073    
 void Rectangle::addToSystemMatrix(Paso_SystemMatrix* mat,  
        const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,  
        dim_t num_Sol, const vector<double>& array) const  
 {  
     const dim_t numMyCols = mat->pattern->mainPattern->numInput;  
     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;  
   
     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;  
     const index_t* mainBlock_index = mat->mainBlock->pattern->index;  
     double* mainBlock_val = mat->mainBlock->val;  
     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;  
     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;  
     double* col_coupleBlock_val = mat->col_coupleBlock->val;  
     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;  
     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;  
     double* row_coupleBlock_val = mat->row_coupleBlock->val;  
   
     for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {  
         // down columns of array  
         const dim_t j_Eq = nodes_Eq[k_Eq];  
         const dim_t i_row = j_Eq;  
 //printf("row:%d\n", i_row);  
         // only look at the matrix rows stored on this processor  
         if (i_row < numMyRows) {  
             for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {  
                 const dim_t i_col = nodes_Sol[k_Sol];  
                 if (i_col < numMyCols) {  
                     for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {  
                         if (mainBlock_index[k] == i_col) {  
                             mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 } else {  
                     for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {  
 //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;  
                         if (col_coupleBlock_index[k] == i_col - numMyCols) {  
                             col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 }  
             }  
         } else {  
             for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {  
                 // across rows of array  
                 const dim_t i_col = nodes_Sol[k_Sol];  
                 if (i_col < numMyCols) {  
                     for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];  
                          k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)  
                     {  
 //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;  
                         if (row_coupleBlock_index[k] == i_col) {  
                             row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];  
                             break;  
                         }  
                     }  
                 }  
             }  
         }  
     }  
 }  
   
2074  } // end of namespace ripley  } // end of namespace ripley
2075    

Legend:
Removed from v.3752  
changed lines
  Added in v.3756

  ViewVC Help
Powered by ViewVC 1.1.26