/[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 3750 by caltinay, Fri Dec 23 01:20:34 2011 UTC revision 3752 by caltinay, Tue Jan 3 08:06:36 2012 UTC
# Line 43  Rectangle::Rectangle(int n0, int n1, dou Line 43  Rectangle::Rectangle(int n0, int n1, dou
43      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
44          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
45    
46      if (n0%m_NX > 0 || n1%m_NY > 0)      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)
47          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");
48    
49      // local number of elements      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
50      m_NE0 = n0/m_NX;          throw RipleyException("Too few elements for the number of ranks");
51      m_NE1 = n1/m_NY;  
52      // local number of nodes (not necessarily owned)      // local number of elements (including overlap)
53        m_NE0 = (m_NX>1 ? (n0+1)/m_NX : n0);
54        if (m_mpiInfo->rank%m_NX>0 && m_mpiInfo->rank%m_NX<m_NX-1)
55            m_NE0++;
56        m_NE1 = (m_NY>1 ? (n1+1)/m_NY : n1);
57        if (m_mpiInfo->rank/m_NX>0 && m_mpiInfo->rank/m_NX<m_NY-1)
58            m_NE1++;
59    
60        // local number of nodes
61      m_N0 = m_NE0+1;      m_N0 = m_NE0+1;
62      m_N1 = m_NE1+1;      m_N1 = m_NE1+1;
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 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset0 = (n0+1)/m_NX*(m_mpiInfo->rank%m_NX);
66      m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);      if (m_mpiInfo->rank%m_NX>0)
67            m_offset0--;
68        m_offset1 = (n1+1)/m_NY*(m_mpiInfo->rank/m_NX);
69        if (m_mpiInfo->rank/m_NX>0)
70            m_offset1--;
71    
72      populateSampleIds();      populateSampleIds();
73  }  }
74    
# Line 250  const int* Rectangle::borrowSampleRefere Line 264  const int* Rectangle::borrowSampleRefere
264  bool Rectangle::ownSample(int fsCode, index_t id) const  bool Rectangle::ownSample(int fsCode, index_t id) const
265  {  {
266  #ifdef ESYS_MPI  #ifdef ESYS_MPI
267      if (fsCode != ReducedNodes) {      if (fsCode == Nodes) {
268          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];          const index_t left = (m_offset0==0 ? 0 : 1);
269          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;          const index_t bottom = (m_offset1==0 ? 0 : 1);
270          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);
271            const index_t top = (m_mpiInfo->rank/m_NX==m_NY-1 ? m_N1 : m_N1-1);
272            const index_t x=id%m_N0;
273            const index_t y=id/m_N0;
274            return (x>=left && x<right && y>=bottom && y<top);
275      } else {      } else {
276          stringstream msg;          stringstream msg;
277          msg << "ownSample() not implemented for "          msg << "ownSample() not implemented for "
# Line 742  Paso_SystemMatrixPattern* Rectangle::get Line 760  Paso_SystemMatrixPattern* Rectangle::get
760      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
761      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
762      const IndexVector faces=getNumFacesPerBoundary();      const IndexVector faces=getNumFacesPerBoundary();
763        const index_t nDOF0 = (m_gNE0+1)/m_NX;
764        const index_t nDOF1 = (m_gNE1+1)/m_NY;
765        const int numDOF=nDOF0*nDOF1;
766      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
767      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
     const int numDOF=getNumDOF();  
768      vector<IndexVector> colIndices(numDOF); // for the couple blocks      vector<IndexVector> colIndices(numDOF); // for the couple blocks
769      int dofCounter=numDOF;      int numShared=0;
770    
771        m_dofMap.assign(getNumNodes(), 0);
772    #pragma omp parallel for
773        for (index_t i=bottom; i<m_N1; i++) {
774            for (index_t j=left; j<m_N0; j++) {
775                m_dofMap[i*m_N0+j]=(i-bottom)*nDOF0+j-left;
776            }
777        }
778    
779      // corner node from bottom-left      // corner node from bottom-left
780      if (faces[0] == 0 && faces[2] == 0) {      if (faces[0] == 0 && faces[2] == 0) {
781          neighbour.push_back(m_mpiInfo->rank-m_NX-1);          neighbour.push_back(m_mpiInfo->rank-m_NX-1);
782          offsetInShared.push_back(offsetInShared.back()+1);          offsetInShared.push_back(offsetInShared.back()+1);
783          sendShared.push_back(-1);          sendShared.push_back(0);
784          recvShared.push_back(dofCounter);          recvShared.push_back(numDOF+numShared);
785          colIndices[0].push_back(dofCounter-numDOF);          colIndices[0].push_back(numShared);
786          ++dofCounter;          m_dofMap[0]=numDOF+numShared;
787            ++numShared;
788      }      }
789      // bottom edge      // bottom edge
790      if (faces[2] == 0) {      if (faces[2] == 0) {
791          neighbour.push_back(m_mpiInfo->rank-m_NX);          neighbour.push_back(m_mpiInfo->rank-m_NX);
792          offsetInShared.push_back(offsetInShared.back()+m_N0-left);          offsetInShared.push_back(offsetInShared.back()+nDOF0);
793          for (dim_t i=0; i<m_N0-left; i++, dofCounter++) {          for (dim_t i=0; i<nDOF0; i++, numShared++) {
794              sendShared.push_back(i);              sendShared.push_back(i);
795              recvShared.push_back(dofCounter);              recvShared.push_back(numDOF+numShared);
796              if (i>0)              if (i>0)
797                  colIndices[i-1].push_back(dofCounter-numDOF);                  colIndices[i-1].push_back(numShared);
798              colIndices[i].push_back(dofCounter-numDOF);              colIndices[i].push_back(numShared);
799              if (i<m_N0-left-1)              if (i<nDOF0-1)
800                  colIndices[i+1].push_back(dofCounter-numDOF);                  colIndices[i+1].push_back(numShared);
801                m_dofMap[i+left]=numDOF+numShared;
802          }          }
803      }      }
804      // corner node from bottom-right      // corner node from bottom-right
805      if (faces[1] == 0 && faces[2] == 0) {      if (faces[1] == 0 && faces[2] == 0) {
806          neighbour.push_back(m_mpiInfo->rank-m_NX+1);          neighbour.push_back(m_mpiInfo->rank-m_NX+1);
807          offsetInShared.push_back(offsetInShared.back()+1);          offsetInShared.push_back(offsetInShared.back()+1);
808          sendShared.push_back(-1);          sendShared.push_back(nDOF0-1);
809          recvShared.push_back(dofCounter);          recvShared.push_back(numDOF+numShared);
810          colIndices[m_N0-left-2].push_back(dofCounter-numDOF);          colIndices[nDOF0-2].push_back(numShared);
811          colIndices[m_N0-left-1].push_back(dofCounter-numDOF);          colIndices[nDOF0-1].push_back(numShared);
812          ++dofCounter;          m_dofMap[m_N0-1]=numDOF+numShared;
813            ++numShared;
814      }      }
815      // right edge      // right edge
816      if (faces[1] == 0) {      if (faces[1] == 0) {
817          neighbour.push_back(m_mpiInfo->rank+1);          neighbour.push_back(m_mpiInfo->rank+1);
818          offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);          offsetInShared.push_back(offsetInShared.back()+nDOF1);
819          for (dim_t i=0; i<m_N1-bottom; i++, dofCounter++) {          for (dim_t i=0; i<nDOF1; i++, numShared++) {
820              sendShared.push_back((i+1)*(m_N0-left)-1);              sendShared.push_back((i+1)*(nDOF0)-1);
821              recvShared.push_back(dofCounter);              recvShared.push_back(numDOF+numShared);
822              if (i>0)              if (i>0)
823                  colIndices[i*(m_N0-left)-1].push_back(dofCounter-numDOF);                  colIndices[i*(nDOF0)-1].push_back(numShared);
824              colIndices[(i+1)*(m_N0-left)-1].push_back(dofCounter-numDOF);              colIndices[(i+1)*(nDOF0)-1].push_back(numShared);
825              if (i<m_N1-bottom-1)              if (i<nDOF1-1)
826                  colIndices[(i+2)*(m_N0-left)-1].push_back(dofCounter-numDOF);                  colIndices[(i+2)*(nDOF0)-1].push_back(numShared);
827                m_dofMap[(i+bottom+1)*m_N0-1]=numDOF+numShared;
828          }          }
829      }      }
830      // corner node from top-right      // corner node from top-right
# Line 800  Paso_SystemMatrixPattern* Rectangle::get Line 832  Paso_SystemMatrixPattern* Rectangle::get
832          neighbour.push_back(m_mpiInfo->rank+m_NX+1);          neighbour.push_back(m_mpiInfo->rank+m_NX+1);
833          offsetInShared.push_back(offsetInShared.back()+1);          offsetInShared.push_back(offsetInShared.back()+1);
834          sendShared.push_back(numDOF-1);          sendShared.push_back(numDOF-1);
835          recvShared.push_back(dofCounter);          recvShared.push_back(numDOF+numShared);
836          colIndices[numDOF-1].push_back(dofCounter-numDOF);          colIndices[numDOF-1].push_back(numShared);
837          ++dofCounter;          m_dofMap[m_N0*m_N1-1]=numDOF+numShared;
838            ++numShared;
839      }      }
840      // top edge      // top edge
841      if (faces[3] == 0) {      if (faces[3] == 0) {
842          neighbour.push_back(m_mpiInfo->rank+m_NX);          neighbour.push_back(m_mpiInfo->rank+m_NX);
843          offsetInShared.push_back(offsetInShared.back()+m_N0-left);          offsetInShared.push_back(offsetInShared.back()+nDOF0);
844          for (dim_t i=0; i<m_N0-left; i++, dofCounter++) {          for (dim_t i=0; i<nDOF0; i++, numShared++) {
845              sendShared.push_back(numDOF-m_N0+left+i);              sendShared.push_back(numDOF-nDOF0+i);
846              recvShared.push_back(dofCounter);              recvShared.push_back(numDOF+numShared);
847              if (i>0)              if (i>0)
848                  colIndices[numDOF-m_N0+left+i-1].push_back(dofCounter-numDOF);                  colIndices[numDOF-nDOF0+i-1].push_back(numShared);
849              colIndices[numDOF-m_N0+left+i].push_back(dofCounter-numDOF);              colIndices[numDOF-nDOF0+i].push_back(numShared);
850              if (i<m_N0-left-1)              if (i<nDOF0-1)
851                  colIndices[numDOF-m_N0+left+i+1].push_back(dofCounter-numDOF);                  colIndices[numDOF-nDOF0+i+1].push_back(numShared);
852                m_dofMap[m_N0*(m_N1-1)+i+left]=numDOF+numShared;
853          }          }
854      }      }
855      // corner node from top-left      // corner node from top-left
856      if (faces[0] == 0 && faces[3] == 0) {      if (faces[0] == 0 && faces[3] == 0) {
857          neighbour.push_back(m_mpiInfo->rank+m_NX-1);          neighbour.push_back(m_mpiInfo->rank+m_NX-1);
858          offsetInShared.push_back(offsetInShared.back()+1);          offsetInShared.push_back(offsetInShared.back()+1);
859          sendShared.push_back(numDOF-m_N0+left);          sendShared.push_back(numDOF-nDOF0);
860          recvShared.push_back(dofCounter);          recvShared.push_back(numDOF+numShared);
861          colIndices[numDOF-m_N0+left].push_back(dofCounter-numDOF);          colIndices[numDOF-nDOF0].push_back(numShared);
862          ++dofCounter;          m_dofMap[m_N0*(m_N1-1)]=numDOF+numShared;
863            ++numShared;
864      }      }
865      // left edge      // left edge
866      if (faces[0] == 0) {      if (faces[0] == 0) {
867          neighbour.push_back(m_mpiInfo->rank-1);          neighbour.push_back(m_mpiInfo->rank-1);
868          offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);          offsetInShared.push_back(offsetInShared.back()+nDOF1);
869          for (dim_t i=0; i<m_N1-bottom; i++, dofCounter++) {          for (dim_t i=0; i<nDOF1; i++, numShared++) {
870              sendShared.push_back(-1);              sendShared.push_back(i*nDOF0);
871              recvShared.push_back(dofCounter);              recvShared.push_back(numDOF+numShared);
872              if (i>0)              if (i>0)
873                  colIndices[(i-1)*(m_N0-left)].push_back(dofCounter-numDOF);                  colIndices[(i-1)*nDOF0].push_back(numShared);
874              colIndices[i*(m_N0-left)].push_back(dofCounter-numDOF);              colIndices[i*nDOF0].push_back(numShared);
875              if (i<m_N1-bottom-1)              if (i<nDOF1-1)
876                  colIndices[(i+1)*(m_N0-left)].push_back(dofCounter-numDOF);                  colIndices[(i+1)*nDOF0].push_back(numShared);
877                m_dofMap[(i+bottom)*m_N0]=numDOF+numShared;
878          }          }
879      }      }
880      /**/  
881        /*
882      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
883      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;      cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
884      for (size_t i=0; i<neighbour.size(); i++) {      for (size_t i=0; i<neighbour.size(); i++) {
# Line 855  Paso_SystemMatrixPattern* Rectangle::get Line 892  Paso_SystemMatrixPattern* Rectangle::get
892      for (size_t i=0; i<sendShared.size(); i++) {      for (size_t i=0; i<sendShared.size(); i++) {
893          cout << "shared[" << i << "]=" << sendShared[i] << endl;          cout << "shared[" << i << "]=" << sendShared[i] << endl;
894      }      }
895      /**/      */
896    
897      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
898              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 914  Paso_SystemMatrixPattern* Rectangle::get Line 951  Paso_SystemMatrixPattern* Rectangle::get
951      copy(index.begin(), index.end(), indexC);      copy(index.begin(), index.end(), indexC);
952      copy(ptr.begin(), ptr.end(), ptrC);      copy(ptr.begin(), ptr.end(), ptrC);
953      M=ptr.size()-1;      M=ptr.size()-1;
954      N=dofCounter-numDOF;      N=numShared;
955      Paso_Pattern *colCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *colCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
956              M, N, ptrC, indexC);              M, N, ptrC, indexC);
957    
958      /**/      /*
959      cout << "--- colCouple_pattern ---" << endl;      cout << "--- colCouple_pattern ---" << endl;
960      cout << "M=" << M << ", N=" << N << endl;      cout << "M=" << M << ", N=" << N << endl;
961      for (size_t i=0; i<ptr.size(); i++) {      for (size_t i=0; i<ptr.size(); i++) {
# Line 927  Paso_SystemMatrixPattern* Rectangle::get Line 964  Paso_SystemMatrixPattern* Rectangle::get
964      for (size_t i=0; i<index.size(); i++) {      for (size_t i=0; i<index.size(); i++) {
965          cout << "index[" << i << "]=" << index[i] << endl;          cout << "index[" << i << "]=" << index[i] << endl;
966      }      }
967      /**/      */
968    
969      // now build the row couple pattern      // now build the row couple pattern
970      IndexVector ptr2(1,0);      IndexVector ptr2(1,0);
# Line 954  Paso_SystemMatrixPattern* Rectangle::get Line 991  Paso_SystemMatrixPattern* Rectangle::get
991      Paso_Pattern *rowCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *rowCouplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
992              N, M, ptrC, indexC);              N, M, ptrC, indexC);
993    
994      /**/      /*
995      cout << "--- rowCouple_pattern ---" << endl;      cout << "--- rowCouple_pattern ---" << endl;
996      cout << "M=" << N << ", N=" << M << endl;      cout << "M=" << N << ", N=" << M << endl;
997      for (size_t i=0; i<ptr2.size(); i++) {      for (size_t i=0; i<ptr2.size(); i++) {
# Line 963  Paso_SystemMatrixPattern* Rectangle::get Line 1000  Paso_SystemMatrixPattern* Rectangle::get
1000      for (size_t i=0; i<index2.size(); i++) {      for (size_t i=0; i<index2.size(); i++) {
1001          cout << "index[" << i << "]=" << index2[i] << endl;          cout << "index[" << i << "]=" << index2[i] << endl;
1002      }      }
1003      /**/      */
1004    
1005      // allocate paso distribution      // allocate paso distribution
1006      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
# Line 1044  pair<double,double> Rectangle::getFirstC Line 1081  pair<double,double> Rectangle::getFirstC
1081  //protected  //protected
1082  dim_t Rectangle::getNumDOF() const  dim_t Rectangle::getNumDOF() const
1083  {  {
1084      return m_nodeDistribution[m_mpiInfo->rank+1]      return (m_gNE0+1)/m_NX*(m_gNE1+1)/m_NY;
         -m_nodeDistribution[m_mpiInfo->rank];  
1085  }  }
1086    
1087  //protected  //protected
# Line 1084  void Rectangle::assembleCoordinates(escr Line 1120  void Rectangle::assembleCoordinates(escr
1120  //private  //private
1121  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
1122  {  {
1123      // identifiers are ordered from left to right, bottom to top on each rank,      // identifiers are ordered from left to right, bottom to top globablly.
     // except for the shared nodes which are owned by the rank below / to the  
     // left of the current rank  
1124    
1125      // build node distribution vector first.      // build node distribution vector first.
     // m_nodeDistribution[i] is the first node id on rank i, that is  
1126      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1127      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1128      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
1129      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
1130          const index_t x=k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y=k/m_NX;  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1;  
         if (y>0)  
             numNodes-=m_N0;  
         if (x>0 && y>0)  
             numNodes++; // subtracted corner twice -> fix that  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
1131      }      }
1132      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
   
     m_dofId.resize(getNumDOF());  
1133      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
1134    
     // the bottom row and left column are not owned by this rank so the  
     // identifiers need to be computed accordingly  
     const index_t left = (m_offset0==0 ? 0 : 1);  
     const index_t bottom = (m_offset1==0 ? 0 : 1);  
     if (left>0) {  
         const int neighbour=m_mpiInfo->rank-1;  
         const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
1135  #pragma omp parallel for  #pragma omp parallel for
1136          for (dim_t i1=bottom; i1<m_N1; i1++) {      for (dim_t i1=0; i1<m_N1; i1++) {
1137              m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]          for (dim_t i0=0; i0<m_N0; i0++) {
1138                  + (i1-bottom+1)*leftN0-1;              m_nodeId[i0+i1*m_N0] = (m_offset1+i1)*(m_gNE0+1)+m_offset0+i0;
1139          }          }
1140      }      }
     if (bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX;  
         const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  
 #pragma omp parallel for  
         for (dim_t i0=left; i0<m_N0; i0++) {  
             m_nodeId[i0]=m_nodeDistribution[neighbour]  
                 + (bottomN1-1)*bottomN0 + i0 - left;  
         }  
     }  
     if (left>0 && bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);  
         m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;  
     }  
1141    
1142      // the rest of the id's are contiguous      m_dofId.resize(numDOF);
1143      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
1144  #pragma omp parallel for      for (dim_t i=0; i<numDOF; i++)
1145      for (dim_t i1=bottom; i1<m_N1; i1++) {          m_dofId[i] = firstId+i;
1146          for (dim_t i0=left; i0<m_N0; i0++) {  
             const index_t idx=i0-left+(i1-bottom)*(m_N0-left);  
             m_nodeId[i0+i1*m_N0] = firstId+idx;  
             m_dofId[idx] = firstId+idx;  
         }  
     }  
1147      m_nodeTags.assign(getNumNodes(), 0);      m_nodeTags.assign(getNumNodes(), 0);
1148      updateTagsInUse(Nodes);      updateTagsInUse(Nodes);
1149    
# Line 1193  void Rectangle::populateSampleIds() Line 1187  void Rectangle::populateSampleIds()
1187  //private  //private
1188  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const  int Rectangle::insertNeighbours(IndexVector& index, index_t node) const
1189  {  {
1190      const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);      const index_t myN0 = (m_gNE0+1)/m_NX;
1191      const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);      const index_t myN1 = (m_gNE1+1)/m_NY;
1192      const int x=node%myN0;      const int x=node%myN0;
1193      const int y=node/myN0;      const int y=node/myN0;
1194      int num=0;      int num=0;
# Line 1412  void Rectangle::nodesToDOF(escript::Data Line 1406  void Rectangle::nodesToDOF(escript::Data
1406    
1407      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1408      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1409        const index_t right = (m_mpiInfo->rank%m_NX==m_NX-1 ? m_N0 : m_N0-1);
1410        const index_t top = (m_mpiInfo->rank/m_NX==m_NY-1 ? m_N1 : m_N1-1);
1411      index_t n=0;      index_t n=0;
1412      for (index_t i=bottom; i<m_N1; i++) {      for (index_t i=bottom; i<top; i++) {
1413          for (index_t j=left; j<m_N0; j++, n++) {          for (index_t j=left; j<right; j++, n++) {
1414              const double* src=in.getSampleDataRO(j+i*m_N0);              const double* src=in.getSampleDataRO(j+i*m_N0);
1415              copy(src, src+numComp, out.getSampleDataRW(n));              copy(src, src+numComp, out.getSampleDataRW(n));
1416          }          }
# Line 2132  void Rectangle::assemblePDESingle(Paso_S Line 2128  void Rectangle::assemblePDESingle(Paso_S
2128                      /* GENERATOR SNIP_PDE_SINGLE BOTTOM */                      /* GENERATOR SNIP_PDE_SINGLE BOTTOM */
2129    
2130                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)                      // add to matrix (if add_EM_S) and RHS (if add_EM_F)
2131                      const index_t left = (m_offset0==0 ? 0 : 1);                      const index_t firstNode=m_N0*k1+k0;
                     const index_t bottom = (m_offset1==0 ? 0 : 1);  
                     const index_t firstNode=(m_N0-bottom)*k1+k0-left;  
                     const int numDOF=getNumDOF();  
                     int dof=numDOF;  
   
2132                      IndexVector rowIndex;                      IndexVector rowIndex;
2133                      if (m_offset0==0) {                      rowIndex.push_back(m_dofMap[firstNode]);
2134                          if (m_offset1>0 && k1==0) {                      rowIndex.push_back(m_dofMap[firstNode+1]);
2135                              rowIndex.push_back(dof++);                      rowIndex.push_back(m_dofMap[firstNode+m_N0]);
2136                              rowIndex.push_back(dof++);                      rowIndex.push_back(m_dofMap[firstNode+m_N0+1]);
                             rowIndex.push_back(firstNode);  
                             rowIndex.push_back(firstNode+1);  
                         } else {  
                             rowIndex.push_back(firstNode);  
                             rowIndex.push_back(firstNode+1);  
                             rowIndex.push_back(firstNode+m_N0-left);  
                             rowIndex.push_back(firstNode+m_N0-left+1);  
                         }  
                     } else {  
                         if (k0==0) {  
                             rowIndex.push_back(dof++);  
                             if (m_offset1>0 && k1==0)  
                                 rowIndex.push_back(dof++);  
                             else  
                                 rowIndex.push_back(firstNode);  
                             rowIndex.push_back(dof++);  
                             rowIndex.push_back(firstNode+m_N0-left);  
                         } else {  
                             rowIndex.push_back(firstNode);  
                             rowIndex.push_back(firstNode+1);  
                             rowIndex.push_back(firstNode+m_N0-left);  
                             rowIndex.push_back(firstNode+m_N0-left+1);  
                         }  
                     }  
                     Paso_Coupler_startCollect(mat->col_coupler, &EM_F[0]);  
                     double* recvBuffer=Paso_Coupler_finishCollect(mat->col_coupler);  
                     for (int i=0; i<mat->col_coupler->connector->send->numSharedComponents; i++) {  
                         if (mat->col_coupler->connector->send->shared[i]>=0)  
                             EM_F[mat->col_coupler->connector->send->shared[i]]+=recvBuffer[i];  
                     }  
   
                     IndexVector colIndex=rowIndex;  
                     /*  
                     Paso_Coupler* coupler=Paso_Coupler_alloc(mat->col_coupler->connector, 4);  
                     Paso_Coupler_startCollect(coupler, &EM_S[0]);  
                     recvBuffer=Paso_Coupler_finishCollect(coupler);  
                     for (int i=0; i<coupler->connector->recv->numSharedComponents; i++) {  
                         if (coupler->connector->recv->shared[i]>=0) {  
                             for (int j=0; j<4; j++)  
                                 EM_S.push_back(recvBuffer[i+4*j]);  
                             colIndex.push_back(coupler->connector->recv->shared[i]);  
                         }  
                     }  
                     Paso_Coupler_free(coupler);*/  
   
2137                      if (add_EM_F) {                      if (add_EM_F) {
2138                          cout << "-- AddtoRHS -- " << endl;                          //cout << "-- AddtoRHS -- " << endl;
2139                          double *F_p=rhs.getSampleDataRW(0);                          double *F_p=rhs.getSampleDataRW(0);
                         cout << "F:"<<endl;  
2140                          for (index_t i=0; i<4; i++) {                          for (index_t i=0; i<4; i++) {
2141                              if (rowIndex[i]<getNumDOF()) {                              if (rowIndex[i]<getNumDOF()) {
2142                                  F_p[rowIndex[i]]+=EM_F[i];                                  F_p[rowIndex[i]]+=EM_F[i];
2143                                  cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;                                  //cout << "F[" << rowIndex[i] << "]=" << F_p[rowIndex[i]] << endl;
2144                              }                              }
2145                          }                          }
2146                          cout << "---"<<endl;                          //cout << "---"<<endl;
2147                      }                      }
2148                      if (add_EM_S) {                      if (add_EM_S) {
2149                          cout << "-- AddtoSystem -- " << endl;                          //cout << "-- AddtoSystem -- " << endl;
2150                          addToSystemMatrix(mat, rowIndex, 1, colIndex, 1, EM_S);                          addToSystemMatrix(mat, rowIndex, 1, rowIndex, 1, EM_S);
2151                      }                      }
2152    
2153                  } // end k0 loop                  } // end k0 loop
# Line 2232  void Rectangle::addToSystemMatrix(Paso_S Line 2177  void Rectangle::addToSystemMatrix(Paso_S
2177          // down columns of array          // down columns of array
2178          const dim_t j_Eq = nodes_Eq[k_Eq];          const dim_t j_Eq = nodes_Eq[k_Eq];
2179          const dim_t i_row = j_Eq;          const dim_t i_row = j_Eq;
2180  printf("row:%d\n", i_row);  //printf("row:%d\n", i_row);
2181          // only look at the matrix rows stored on this processor          // only look at the matrix rows stored on this processor
2182          if (i_row < numMyRows) {          if (i_row < numMyRows) {
2183              for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {              for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
# Line 2246  printf("row:%d\n", i_row); Line 2191  printf("row:%d\n", i_row);
2191                      }                      }
2192                  } else {                  } else {
2193                      for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {                      for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
2194  cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;  //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
2195                          if (col_coupleBlock_index[k] == i_col - numMyCols) {                          if (col_coupleBlock_index[k] == i_col - numMyCols) {
2196                              col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];                              col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
2197                              break;                              break;
# Line 2262  cout << "col:" << i_col-numMyCols << " c Line 2207  cout << "col:" << i_col-numMyCols << " c
2207                      for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];                      for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
2208                           k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)                           k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
2209                      {                      {
2210  cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;  //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
2211                          if (row_coupleBlock_index[k] == i_col) {                          if (row_coupleBlock_index[k] == i_col) {
2212                              row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];                              row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
2213                              break;                              break;

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

  ViewVC Help
Powered by ViewVC 1.1.26