/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3753 by caltinay, Tue Jan 3 09:01:49 2012 UTC revision 3754 by caltinay, Wed Jan 4 07:07:37 2012 UTC
# Line 776  Paso_SystemMatrixPattern* Rectangle::get Line 776  Paso_SystemMatrixPattern* Rectangle::get
776          }          }
777      }      }
778    
779      // corner node from bottom-left      // build list of shared components and neighbours by looping through
780      if (faces[0] == 0 && faces[2] == 0) {      // all potential neighbouring ranks and checking if positions are
781          neighbour.push_back(m_mpiInfo->rank-m_NX-1);      // within bounds
782          offsetInShared.push_back(offsetInShared.back()+1);      const int x=m_mpiInfo->rank%m_NX;
783          sendShared.push_back(0);      const int y=m_mpiInfo->rank/m_NX;
784          recvShared.push_back(numDOF+numShared);      for (int i1=-1; i1<2; i1++) {
785          colIndices[0].push_back(numShared);          for (int i0=-1; i0<2; i0++) {
786          m_dofMap[0]=numDOF+numShared;              // skip rank itself
787          ++numShared;              if (i0==0 && i1==0)
788      }                  continue;
789      // bottom edge              // location of neighbour rank
790      if (faces[2] == 0) {              const int nx=x+i0;
791          neighbour.push_back(m_mpiInfo->rank-m_NX);              const int ny=y+i1;
792          offsetInShared.push_back(offsetInShared.back()+nDOF0);              if (nx>=0 && ny>=0 && nx<m_NX && ny<m_NY) {
793          for (dim_t i=0; i<nDOF0; i++, numShared++) {                  neighbour.push_back(ny*m_NX+nx);
794              sendShared.push_back(i);                  if (i0==0) {
795              recvShared.push_back(numDOF+numShared);                      // sharing top or bottom edge
796              if (i>0)                      const int firstDOF=(i1==-1 ? 0 : numDOF-nDOF0);
797                  colIndices[i-1].push_back(numShared);                      const int firstNode=(i1==-1 ? left : m_N0*(m_N1-1)+left);
798              colIndices[i].push_back(numShared);                      offsetInShared.push_back(offsetInShared.back()+nDOF0);
799              if (i<nDOF0-1)                      for (dim_t i=0; i<nDOF0; i++, numShared++) {
800                  colIndices[i+1].push_back(numShared);                          sendShared.push_back(firstDOF+i);
801              m_dofMap[i+left]=numDOF+numShared;                          recvShared.push_back(numDOF+numShared);
802          }                          if (i>0)
803      }                              colIndices[firstDOF+i-1].push_back(numShared);
804      // corner node from bottom-right                          colIndices[firstDOF+i].push_back(numShared);
805      if (faces[1] == 0 && faces[2] == 0) {                          if (i<nDOF0-1)
806          neighbour.push_back(m_mpiInfo->rank-m_NX+1);                              colIndices[firstDOF+i+1].push_back(numShared);
807          offsetInShared.push_back(offsetInShared.back()+1);                          m_dofMap[firstNode+i]=numDOF+numShared;
808          sendShared.push_back(nDOF0-1);                      }
809          recvShared.push_back(numDOF+numShared);                  } else if (i1==0) {
810          colIndices[nDOF0-2].push_back(numShared);                      // sharing left or right edge
811          colIndices[nDOF0-1].push_back(numShared);                      const int firstDOF=(i0==-1 ? 0 : nDOF0-1);
812          m_dofMap[m_N0-1]=numDOF+numShared;                      const int firstNode=(i0==-1 ? bottom*m_N0 : (bottom+1)*m_N0-1);
813          ++numShared;                      offsetInShared.push_back(offsetInShared.back()+nDOF1);
814      }                      for (dim_t i=0; i<nDOF1; i++, numShared++) {
815      // right edge                          sendShared.push_back(firstDOF+i*nDOF0);
816      if (faces[1] == 0) {                          recvShared.push_back(numDOF+numShared);
817          neighbour.push_back(m_mpiInfo->rank+1);                          if (i>0)
818          offsetInShared.push_back(offsetInShared.back()+nDOF1);                              colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);
819          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          colIndices[firstDOF+i*nDOF0].push_back(numShared);
820              sendShared.push_back((i+1)*(nDOF0)-1);                          if (i<nDOF1-1)
821              recvShared.push_back(numDOF+numShared);                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);
822              if (i>0)                          m_dofMap[firstNode+i*m_N0]=numDOF+numShared;
823                  colIndices[i*(nDOF0)-1].push_back(numShared);                      }
824              colIndices[(i+1)*(nDOF0)-1].push_back(numShared);                  } else {
825              if (i<nDOF1-1)                      // sharing a node
826                  colIndices[(i+2)*(nDOF0)-1].push_back(numShared);                      const int dof=(i0+1)/2*(nDOF0-1)+(i1+1)/2*(numDOF-nDOF0);
827              m_dofMap[(i+bottom+1)*m_N0-1]=numDOF+numShared;                      const int node=(i0+1)/2*(m_N0-1)+(i1+1)/2*m_N0*(m_N1-1);
828          }                      offsetInShared.push_back(offsetInShared.back()+1);
829      }                      sendShared.push_back(dof);
830      // corner node from top-right                      recvShared.push_back(numDOF+numShared);
831      if (faces[1] == 0 && faces[3] == 0) {                      colIndices[dof].push_back(numShared);
832          neighbour.push_back(m_mpiInfo->rank+m_NX+1);                      m_dofMap[node]=numDOF+numShared;
833          offsetInShared.push_back(offsetInShared.back()+1);                      ++numShared;
834          sendShared.push_back(numDOF-1);                  }
835          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;  
836          }          }
837      }      }
838    
# Line 892  Paso_SystemMatrixPattern* Rectangle::get Line 850  Paso_SystemMatrixPattern* Rectangle::get
850      for (size_t i=0; i<sendShared.size(); i++) {      for (size_t i=0; i<sendShared.size(); i++) {
851          cout << "shared[" << i << "]=" << sendShared[i] << endl;          cout << "shared[" << i << "]=" << sendShared[i] << endl;
852      }      }
853        cout << "--- dofMap ---" << endl;
854        for (size_t i=0; i<m_dofMap.size(); i++) {
855            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
856        }
857        cout << "--- colIndices ---" << endl;
858        for (size_t i=0; i<colIndices.size(); i++) {
859            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
860        }
861      */      */
862    
863      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
# Line 1191  void Rectangle::populateSampleIds() Line 1157  void Rectangle::populateSampleIds()
1157  //private  //private
1158  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const
1159  {  {
1160      const index_t myN0 = (m_gNE0+1)/m_NX;      const index_t nDOF0 = (m_gNE0+1)/m_NX;
1161      const index_t myN1 = (m_gNE1+1)/m_NY;      const index_t nDOF1 = (m_gNE1+1)/m_NY;
1162      const int x=node%myN0;      const int x=node%nDOF0;
1163      const int y=node/myN0;      const int y=node/nDOF0;
1164      int num=0;      int num=0;
1165      if (y>0) {      // loop through potential neighbours and add to index if positions are
1166          if (x>0) {      // within bounds
1167              // bottom-left      for (int i1=-1; i1<2; i1++) {
1168              index.push_back(node-myN0-1);          for (int i0=-1; i0<2; i0++) {
1169              num++;              // skip node itself
1170          }              if (i0==0 && i1==0)
1171          // bottom                  continue;
1172          index.push_back(node-myN0);              // location of neighbour node
1173          num++;              const int nx=x+i0;
1174          if (x<myN0-1) {              const int ny=y+i1;
1175              // bottom-right              if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
1176              index.push_back(node-myN0+1);                  index.push_back(ny*nDOF0+nx);
1177              num++;                  num++;
1178                }
1179          }          }
1180      }      }
1181      if (x<myN0-1) {  
1182          // right      return num;
1183          index.push_back(node+1);  }
1184          num++;  
1185          if (y<myN1-1) {  //protected
1186              // top-right  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const
1187              index.push_back(node+myN0+1);  {
1188              num++;      const dim_t numComp = in.getDataPointSize();
1189        out.requireWrite();
1190    
1191        const index_t left = (m_offset0==0 ? 0 : 1);
1192        const index_t bottom = (m_offset1==0 ? 0 : 1);
1193        const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);
1194        const index_t top = (m_mpiInfo->rank/m_NX==m_NY-1 ? m_N1 : m_N1-1);
1195        index_t n=0;
1196        for (index_t i=bottom; i<top; i++) {
1197            for (index_t j=left; j<right; j++, n++) {
1198                const double* src=in.getSampleDataRO(j+i*m_N0);
1199                copy(src, src+numComp, out.getSampleDataRW(n));
1200          }          }
1201      }      }
1202      if (y<myN1-1) {  }
1203          // top  
1204          index.push_back(node+myN0);  //protected
1205          num++;  void Rectangle::addToSystemMatrix(Paso_SystemMatrix* mat,
1206          if (x>0) {         const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
1207              // top-left         dim_t num_Sol, const vector<double>& array) const
1208              index.push_back(node+myN0-1);  {
1209              num++;      const dim_t numMyCols = mat->pattern->mainPattern->numInput;
1210        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
1211    
1212        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
1213        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
1214        double* mainBlock_val = mat->mainBlock->val;
1215        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
1216        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
1217        double* col_coupleBlock_val = mat->col_coupleBlock->val;
1218        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
1219        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
1220        double* row_coupleBlock_val = mat->row_coupleBlock->val;
1221    
1222        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1223            // down columns of array
1224            const dim_t j_Eq = nodes_Eq[k_Eq];
1225            const dim_t i_row = j_Eq;
1226    //printf("row:%d\n", i_row);
1227            // only look at the matrix rows stored on this processor
1228            if (i_row < numMyRows) {
1229                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1230                    const dim_t i_col = nodes_Sol[k_Sol];
1231                    if (i_col < numMyCols) {
1232                        for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
1233                            if (mainBlock_index[k] == i_col) {
1234                                mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
1235                                break;
1236                            }
1237                        }
1238                    } else {
1239                        for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1240    //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
1241                            if (col_coupleBlock_index[k] == i_col - numMyCols) {
1242                                col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
1243                                break;
1244                            }
1245                        }
1246                    }
1247                }
1248            } else {
1249                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1250                    // across rows of array
1251                    const dim_t i_col = nodes_Sol[k_Sol];
1252                    if (i_col < numMyCols) {
1253                        for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1254                             k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1255                        {
1256    //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
1257                            if (row_coupleBlock_index[k] == i_col) {
1258                                row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
1259                                break;
1260                            }
1261                        }
1262                    }
1263                }
1264          }          }
1265      }      }
     if (x>0) {  
         // left  
         index.push_back(node-1);  
         num++;  
     }  
   
     return num;  
1266  }  }
1267    
1268  //protected  //protected
# Line 1403  void Rectangle::interpolateNodesOnFaces( Line 1428  void Rectangle::interpolateNodesOnFaces(
1428  }  }
1429    
1430  //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  
1431  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,  void Rectangle::assemblePDESingle(Paso_SystemMatrix* mat,
1432          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1433          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
# Line 1430  void Rectangle::assemblePDESingle(Paso_S Line 1436  void Rectangle::assemblePDESingle(Paso_S
1436  {  {
1437      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
1438      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
1439      /* GENERATOR SNIP_PDE_SINGLE_PRE TOP */      /*** GENERATOR SNIP_PDE_SINGLE_PRE TOP */
1440      const double w0 = -0.1555021169820365539*h1/h0;      const double w0 = -0.1555021169820365539*h1/h0;
1441      const double w1 = 0.041666666666666666667;      const double w1 = 0.041666666666666666667;
1442      const double w10 = -0.041666666666666666667*h0/h1;      const double w10 = -0.041666666666666666667*h0/h1;
# Line 1512  void Rectangle::assemblePDESingle(Paso_S Line 1518  void Rectangle::assemblePDESingle(Paso_S
1518      rhs.requireWrite();      rhs.requireWrite();
1519  #pragma omp parallel  #pragma omp parallel
1520      {      {
1521          for (index_t k1_0=0; k1_0<2; k1_0++) { // coloring          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
1522  #pragma omp for  #pragma omp for
1523              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {              for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
1524                  for (index_t k0=0; k0<m_NE0; ++k0)  {                  for (index_t k0=0; k0<m_NE0; ++k0)  {
# Line 1521  void Rectangle::assemblePDESingle(Paso_S Line 1527  void Rectangle::assemblePDESingle(Paso_S
1527                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
1528                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
1529                      const index_t e = k0 + m_NE0*k1;                      const index_t e = k0 + m_NE0*k1;
1530                      /* GENERATOR SNIP_PDE_SINGLE TOP */                      /*** GENERATOR SNIP_PDE_SINGLE TOP */
1531                      ///////////////                      ///////////////
1532                      // process A //                      // process A //
1533                      ///////////////                      ///////////////
# Line 2141  void Rectangle::assemblePDESingle(Paso_S Line 2147  void Rectangle::assemblePDESingle(Paso_S
2147                      if (add_EM_F) {                      if (add_EM_F) {
2148                          //cout << "-- AddtoRHS -- " << endl;                          //cout << "-- AddtoRHS -- " << endl;
2149                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
2150                          for (index_t i=0; i<4; i++) {                          for (index_t i=0; i<rowIndex.size(); i++) {
2151                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2152                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
2153                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
# Line 2160  void Rectangle::assemblePDESingle(Paso_S Line 2166  void Rectangle::assemblePDESingle(Paso_S
2166      } // end of parallel region      } // end of parallel region
2167  }  }
2168    
 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;  
                         }  
                     }  
                 }  
             }  
         }  
     }  
 }  
   
2169  } // end of namespace ripley  } // end of namespace ripley
2170    

Legend:
Removed from v.3753  
changed lines
  Added in v.3754

  ViewVC Help
Powered by ViewVC 1.1.26