/[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 3755 by caltinay, Thu Jan 5 06:51:31 2012 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 70  Rectangle::Rectangle(int n0, int n1, dou Line 70  Rectangle::Rectangle(int n0, int n1, dou
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;  
         }  
     }  
   
     // build list of shared components and neighbours by looping through  
     // all potential neighbouring ranks and checking if positions are  
     // within bounds  
     const int x=m_mpiInfo->rank%m_NX;  
     const int y=m_mpiInfo->rank/m_NX;  
     for (int i1=-1; i1<2; i1++) {  
         for (int i0=-1; i0<2; i0++) {  
             // skip rank itself  
             if (i0==0 && i1==0)  
                 continue;  
             // location of neighbour rank  
             const int nx=x+i0;  
             const int ny=y+i1;  
             if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {  
                 neighbour.push_back(ny*m_NX+nx);  
                 if (i0==0) {  
                     // sharing top or bottom edge  
                     const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);  
                     const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);  
                     offsetInShared.push_back(offsetInShared.back()+nDOF0);  
                     for (dim_t i=0; i<nDOF0; i++, numShared++) {  
                         sendShared.push_back(firstDOF+i);  
                         recvShared.push_back(numDOF+numShared);  
                         if (i>0)  
                             colIndices[firstDOF+i-1].push_back(numShared);  
                         colIndices[firstDOF+i].push_back(numShared);  
                         if (i<nDOF0-1)  
                             colIndices[firstDOF+i+1].push_back(numShared);  
                         m_dofMap[firstNode+i]=numDOF+numShared;  
                     }  
                 } else if (i1==0) {  
                     // sharing left or right edge  
                     const int firstDOF=(i0==-1 ? 0 : nDOF0-1);  
                     const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);  
                     offsetInShared.push_back(offsetInShared.back()+nDOF1);  
                     for (dim_t i=0; i<nDOF1; i++, numShared++) {  
                         sendShared.push_back(firstDOF+i*nDOF0);  
                         recvShared.push_back(numDOF+numShared);  
                         if (i>0)  
                             colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);  
                         colIndices[firstDOF+i*nDOF0].push_back(numShared);  
                         if (i<nDOF1-1)  
                             colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);  
                         m_dofMap[firstNode+i*m_N0]=numDOF+numShared;  
                     }  
                 } else {  
                     // sharing a node  
                     const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);  
                     const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);  
                     offsetInShared.push_back(offsetInShared.back()+1);  
                     sendShared.push_back(dof);  
                     recvShared.push_back(numDOF+numShared);  
                     colIndices[dof].push_back(numShared);  
                     m_dofMap[node]=numDOF+numShared;  
                     ++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;  
     }  
     cout << "--- dofMap ---" << endl;  
     for (size_t i=0; i<m_dofMap.size(); i++) {  
         cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;  
     }  
     cout << "--- colIndices ---" << endl;  
     for (size_t i=0; i<colIndices.size(); i++) {  
         cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << 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 1083  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 1158  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 nDOF0 = (m_gNE0+1)/m_NX;      const dim_t nDOF0 = (m_gNE0+1)/m_NX;
1003      const index_t nDOF1 = (m_gNE1+1)/m_NY;      const dim_t nDOF1 = (m_gNE1+1)/m_NY;
1004      const int x=node%nDOF0;      const index_t left = (m_offset0==0 ? 0 : 1);
1005      const int y=node/nDOF0;      const index_t bottom = (m_offset1==0 ? 0 : 1);
1006      int num=0;  
1007      // loop through potential neighbours and add to index if positions are      // populate node->DOF mapping with own degrees of freedom.
1008        // The rest is assigned in the loop further down
1009        m_dofMap.assign(getNumNodes(), 0);
1010    #pragma omp parallel for
1011        for (index_t i=bottom; i<m_N1; i++) {
1012            for (index_t j=left; j<m_N0; j++) {
1013                m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
1014            }
1015        }
1016    
1017        // build list of shared components and neighbours by looping through
1018        // all potential neighbouring ranks and checking if positions are
1019      // within bounds      // within bounds
1020        const dim_t numDOF=nDOF0*nDOF1;
1021        vector<IndexVector> colIndices(numDOF); // for the couple blocks
1022        RankVector neighbour;
1023        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++) {      for (int i1=-1; i1<2; i1++) {
1029          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1030              // skip node itself              // skip this rank
1031              if (i0==0 && i1==0)              if (i0==0 && i1==0)
1032                  continue;                  continue;
1033              // location of neighbour node              // location of neighbour rank
1034              const int nx=x+i0;              const int nx=x+i0;
1035              const int ny=y+i1;              const int ny=y+i1;
1036              if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
1037                  index.push_back(ny*nDOF0+nx);                  neighbour.push_back(ny*m_NX+nx);
1038                  num++;                  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    
1083      return num;      // create connector
1084  }      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1085                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
1086                &offsetInShared[0], 1, 0, m_mpiInfo);
1087        Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
1088                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
1089                &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  //protected      // create main and couple blocks
1095  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const      Paso_Pattern *mainPattern = createMainPattern();
1096  {      Paso_Pattern *colPattern, *rowPattern;
1097      const dim_t numComp = in.getDataPointSize();      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);
     out.requireWrite();  
1098    
1099      const index_t left = (m_offset0==0 ? 0 : 1);      // allocate paso distribution
1100      const index_t bottom = (m_offset1==0 ? 0 : 1);      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
1101      const index_t nDOF0 = (m_gNE0+1)/m_NX;              const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
     const index_t nDOF1 = (m_gNE1+1)/m_NY;  
 #pragma omp parallel for  
     for (index_t i=0; i<nDOF1; i++) {  
         for (index_t j=0; j<nDOF0; j++) {  
             const index_t n=j+left+(i+bottom)*m_N0;  
             const double* src=in.getSampleDataRO(n);  
             copy(src, src+numComp, out.getSampleDataRW(j+i*nDOF0));  
         }  
     }  
 }  
1102    
1103  //protected      // finally create the system matrix
1104  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const      m_pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
1105  {              distribution, distribution, mainPattern, colPattern, rowPattern,
1106      const dim_t numComp = in.getDataPointSize();              m_connector, m_connector);
     out.requireWrite();  
1107    
1108      //TODO: use coupler to get the rest of the values      Paso_Distribution_free(distribution);
1109    
1110      const index_t left = (m_offset0==0 ? 0 : 1);      // useful debug output
1111      const index_t bottom = (m_offset1==0 ? 0 : 1);      /*
1112      const index_t nDOF0 = (m_gNE0+1)/m_NX;      cout << "--- rcv_shcomp ---" << endl;
1113      const index_t nDOF1 = (m_gNE1+1)/m_NY;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
1114  #pragma omp parallel for      for (size_t i=0; i<neighbour.size(); i++) {
1115      for (index_t i=0; i<nDOF1; i++) {          cout << "neighbor[" << i << "]=" << neighbour[i]
1116          for (index_t j=0; j<nDOF0; j++) {              << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
             const double* src=in.getSampleDataRO(j+i*nDOF0);  
             const index_t n=j+left+(i+bottom)*m_N0;  
             copy(src, src+numComp, out.getSampleDataRW(n));  
         }  
1117      }      }
1118  }      for (size_t i=0; i<recvShared.size(); i++) {
1119            cout << "shared[" << i << "]=" << recvShared[i] << endl;
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  //protected      /*
1136  void Rectangle::addToSystemMatrix(Paso_SystemMatrix* mat,      cout << "--- main_pattern ---" << endl;
1137         const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,      cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
1138         dim_t num_Sol, const vector<double>& array) const      for (size_t i=0; i<mainPattern->numOutput+1; i++) {
1139  {          cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
1140      const dim_t numMyCols = mat->pattern->mainPattern->numInput;      }
1141      const dim_t numMyRows = mat->pattern->mainPattern->numOutput;      for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
1142            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
1143      const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;      }
1144      const index_t* mainBlock_index = mat->mainBlock->pattern->index;      */
1145      double* mainBlock_val = mat->mainBlock->val;  
1146      const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;      /*
1147      const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;      cout << "--- colCouple_pattern ---" << endl;
1148      double* col_coupleBlock_val = mat->col_coupleBlock->val;      cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
1149      const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;      for (size_t i=0; i<colPattern->numOutput+1; i++) {
1150      const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;          cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
1151      double* row_coupleBlock_val = mat->row_coupleBlock->val;      }
1152        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
1153      for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {          cout << "index[" << i << "]=" << colPattern->index[i] << endl;
         // 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;  
                         }  
                     }  
                 }  
             }  
         }  
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

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

  ViewVC Help
Powered by ViewVC 1.1.26