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

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

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

revision 5118 by caltinay, Mon Aug 18 11:28:55 2014 UTC revision 5119 by caltinay, Wed Aug 20 04:10:01 2014 UTC
# Line 67  Brick::Brick(int n0, int n1, int n2, dou Line 67  Brick::Brick(int n0, int n1, int n2, dou
67               escript::SubWorld_ptr w) :               escript::SubWorld_ptr w) :
68      RipleyDomain(3, w)      RipleyDomain(3, w)
69  {  {
70      if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)      if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)
71              * static_cast<long>(n2 + 1) > numeric_limits<int>::max())              * static_cast<long>(n2 + 1) > numeric_limits<int>::max())
72          throw RipleyException("The number of elements has overflowed, this "          throw RipleyException("The number of elements has overflowed, this "
73                  "limit may be raised in future releases.");                  "limit may be raised in future releases.");
# Line 98  Brick::Brick(int n0, int n1, int n2, dou Line 98  Brick::Brick(int n0, int n1, int n2, dou
98              if (ranks % d[i] != 0) {              if (ranks % d[i] != 0) {
99                  throw RipleyException("Invalid number of spatial subdivisions");                  throw RipleyException("Invalid number of spatial subdivisions");
100              }              }
101              //remove              //remove
102              ranks /= d[i];              ranks /= d[i];
103          }          }
104          factorise(factors, ranks);          factorise(factors, ranks);
105          if (factors.size() != 0) {          if (factors.size() != 0) {
106              warn = true;              warn = true;
107          }          }
108      }      }
109      while (factors.size() > 0) {      while (factors.size() > 0) {
110          int i = indexOfMax(epr[0],epr[1],epr[2]);          int i = indexOfMax(epr[0],epr[1],epr[2]);
111          int f = factors.back();          int f = factors.back();
# Line 203  Brick::Brick(int n0, int n1, int n2, dou Line 203  Brick::Brick(int n0, int n1, int n2, dou
203          m_offset[2]--;          m_offset[2]--;
204    
205      populateSampleIds();      populateSampleIds();
206      createPattern();  
       
207      for (simap_t::const_iterator i = tagnamestonums.begin();      for (simap_t::const_iterator i = tagnamestonums.begin();
208              i != tagnamestonums.end(); i++) {              i != tagnamestonums.end(); i++) {
209          setTagMap(i->first, i->second);          setTagMap(i->first, i->second);
# Line 212  Brick::Brick(int n0, int n1, int n2, dou Line 211  Brick::Brick(int n0, int n1, int n2, dou
211      addPoints(points, tags);      addPoints(points, tags);
212  }  }
213    
   
214  Brick::~Brick()  Brick::~Brick()
215  {  {
216  }  }
# Line 665  void Brick::readBinaryGridZippedImpl(esc Line 663  void Brick::readBinaryGridZippedImpl(esc
663              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
664                               +(idx2+z)*params.numValues[0]*params.numValues[1]);                               +(idx2+z)*params.numValues[0]*params.numValues[1]);
665              memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));              memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
666                
667              for (int x=0; x<num0; x++) {              for (int x=0; x<num0; x++) {
668                  const int baseIndex = first0+x*params.multiplier[0]                  const int baseIndex = first0+x*params.multiplier[0]
669                                       +(first1+y*params.multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
# Line 792  void Brick::dump(const string& fileName) Line 790  void Brick::dump(const string& fileName)
790          fn+=".silo";          fn+=".silo";
791      }      }
792    
793      int driver=DB_HDF5;          int driver=DB_HDF5;
794      string siloPath;      string siloPath;
795      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
796    
# Line 2107  void Brick::dofToNodes(escript::Data& ou Line 2105  void Brick::dofToNodes(escript::Data& ou
2105      }      }
2106  }  }
2107    
2108    //protected
2109    paso::SystemMatrixPattern_ptr Brick::getPasoMatrixPattern(
2110                                                        bool reducedRowOrder,
2111                                                        bool reducedColOrder) const
2112    {
2113        if (m_pattern.get())
2114            return m_pattern;
2115    
2116        // first call to this method -> create the pattern, then return it
2117        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2118        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2119        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2120        const dim_t numDOF = nDOF0*nDOF1*nDOF2;
2121        const dim_t numShared = m_connector->send->numSharedComponents;
2122        const index_t* sendShared = m_connector->send->shared;
2123        const int x = m_mpiInfo->rank%m_NX[0];
2124        const int y = m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2125        const int z = m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
2126    
2127        // these are for the couple blocks
2128        vector<IndexVector> colIndices(numDOF);
2129        vector<IndexVector> rowIndices(numShared);
2130    
2131        for (dim_t i=0; i < m_connector->send->numNeighbors; i++) {
2132            const dim_t start = m_connector->send->offsetInShared[i];
2133            const dim_t end = m_connector->send->offsetInShared[i+1];
2134            // location of neighbour rank relative to this rank
2135            const int xDiff = m_connector->send->neighbor[i]%m_NX[0] - x;
2136            const int yDiff = m_connector->send->neighbor[i]%(m_NX[0]*m_NX[1])/m_NX[0] - y;
2137            const int zDiff = m_connector->send->neighbor[i]/(m_NX[0]*m_NX[1]) - z;
2138            
2139            if (xDiff==0 && yDiff==0) {
2140                // sharing front or back plane
2141                for (dim_t j = start; j < end; j++) {
2142                    const dim_t i0 = (j-start)%nDOF0;
2143                    const dim_t i1 = (j-start)/nDOF0;
2144                    if (i0 > 0) {
2145                        if (i1 > 0)
2146                            doublyLink(colIndices, rowIndices, sendShared[j-1-nDOF0], j);
2147                        doublyLink(colIndices, rowIndices, sendShared[j-1], j);
2148                        if (i1 < nDOF1-1)
2149                            doublyLink(colIndices, rowIndices, sendShared[j-1+nDOF0], j);
2150                    }
2151                    if (i1 > 0)
2152                        doublyLink(colIndices, rowIndices, sendShared[j-nDOF0], j);
2153                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2154                    if (i1 < nDOF1-1)
2155                        doublyLink(colIndices, rowIndices, sendShared[j+nDOF0], j);
2156                    if (i0 < nDOF0-1) {
2157                        if (i1 > 0)
2158                            doublyLink(colIndices, rowIndices, sendShared[j+1-nDOF0], j);
2159                        doublyLink(colIndices, rowIndices, sendShared[j+1], j);
2160                        if (i1 < nDOF1-1)
2161                            doublyLink(colIndices, rowIndices, sendShared[j+1+nDOF0], j);
2162                    }
2163                }
2164            } else if (xDiff==0 && zDiff==0) {
2165                // sharing top or bottom plane
2166                for (dim_t j = start; j < end; j++) {
2167                    const dim_t i0 = (j-start)%nDOF0;
2168                    const dim_t i1 = (j-start)/nDOF0;
2169                    if (i0 > 0) {
2170                        if (i1 > 0)
2171                            doublyLink(colIndices, rowIndices, sendShared[j]-1-nDOF0*nDOF1, j);
2172                        doublyLink(colIndices, rowIndices, sendShared[j]-1, j);
2173                        if (i1 < nDOF2-1)
2174                            doublyLink(colIndices, rowIndices, sendShared[j]-1+nDOF0*nDOF1, j);
2175                    }
2176                    if (i1 > 0)
2177                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2178                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2179                    if (i1 < nDOF2-1)
2180                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2181                    if (i0 < nDOF0-1) {
2182                        if (i1 > 0)
2183                            doublyLink(colIndices, rowIndices, sendShared[j]+1-nDOF0*nDOF1, j);
2184                        doublyLink(colIndices, rowIndices, sendShared[j]+1, j);
2185                        if (i1 < nDOF2-1)
2186                            doublyLink(colIndices, rowIndices, sendShared[j]+1+nDOF0*nDOF1, j);
2187                    }
2188                }
2189            } else if (yDiff==0 && zDiff==0) {
2190                // sharing left or right plane
2191                for (dim_t j = start; j < end; j++) {
2192                    const dim_t i0 = (j-start)%nDOF1;
2193                    const dim_t i1 = (j-start)/nDOF1;
2194                    if (i0 > 0) {
2195                        if (i1 > 0)
2196                            doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0-nDOF0*nDOF1, j);
2197                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0, j);
2198                        if (i1 < nDOF2-1)
2199                            doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0+nDOF0*nDOF1, j);
2200                    }
2201                    if (i1 > 0)
2202                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2203                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2204                    if (i1 < nDOF2-1)
2205                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2206                    if (i0 < nDOF1-1) {
2207                        if (i1 > 0)
2208                            doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0-nDOF0*nDOF1, j);
2209                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0, j);
2210                        if (i1 < nDOF2-1)
2211                            doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0+nDOF0*nDOF1, j);
2212                    }
2213                }
2214            } else if (xDiff == 0) {
2215                // sharing an edge in x direction
2216                for (dim_t j = start; j < end; j++) {
2217                    if (j > start)
2218                        doublyLink(colIndices, rowIndices, sendShared[j]-1, j);
2219                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2220                    if (j < end-1)
2221                        doublyLink(colIndices, rowIndices, sendShared[j]+1, j);
2222                }
2223            } else if (yDiff == 0) {
2224                // sharing an edge in y direction
2225                for (dim_t j = start; j < end; j++) {
2226                    if (j > start)
2227                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0, j);
2228                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2229                    if (j < end-1)
2230                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0, j);
2231                }
2232            } else if (zDiff == 0) {
2233                // sharing an edge in z direction
2234                for (dim_t j = start; j < end; j++) {
2235                    if (j > start)
2236                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2237                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2238                    if (j < end-1)
2239                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2240                }
2241            } else {
2242                // sharing a node
2243                doublyLink(colIndices, rowIndices, sendShared[start], start);
2244            }
2245        }
2246    
2247    #pragma omp parallel for
2248        for (int i = 0; i < numShared; i++) {
2249            sort(rowIndices[i].begin(), rowIndices[i].end());
2250        }
2251    
2252        // create main and couple blocks
2253        paso::Pattern_ptr mainPattern = createPasoPattern(getConnections(), numDOF);
2254        paso::Pattern_ptr colPattern = createPasoPattern(colIndices, numShared);
2255        paso::Pattern_ptr rowPattern = createPasoPattern(rowIndices, numDOF);
2256    
2257        // allocate paso distribution
2258        paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2259                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2260    
2261        // finally create the system matrix pattern
2262        m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2263                distribution, distribution, mainPattern, colPattern, rowPattern,
2264                m_connector, m_connector));
2265    
2266        // useful debug output
2267        /*
2268        cout << "--- colIndices ---" << endl;
2269        for (size_t i=0; i<colIndices.size(); i++) {
2270            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
2271        }
2272        cout << "--- rowIndices ---" << endl;
2273        for (size_t i=0; i<rowIndices.size(); i++) {
2274            cout << "rowIndices[" << i << "].size()=" << rowIndices[i].size() << endl;
2275        }
2276        */
2277        /*
2278        cout << "--- main_pattern ---" << endl;
2279        cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
2280        for (size_t i=0; i<mainPattern->numOutput+1; i++) {
2281            cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
2282        }
2283        for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
2284            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
2285        }
2286        */
2287        /*
2288        cout << "--- colCouple_pattern ---" << endl;
2289        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
2290        for (size_t i=0; i<colPattern->numOutput+1; i++) {
2291            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
2292        }
2293        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
2294            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
2295        }
2296        */
2297        /*
2298        cout << "--- rowCouple_pattern ---" << endl;
2299        cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
2300        for (size_t i=0; i<rowPattern->numOutput+1; i++) {
2301            cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
2302        }
2303        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
2304            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2305        }
2306        */
2307    
2308        return m_pattern;
2309    }
2310    
2311  //private  //private
2312  void Brick::populateSampleIds()  void Brick::populateSampleIds()
2313  {  {
# Line 2125  void Brick::populateSampleIds() Line 2326  void Brick::populateSampleIds()
2326          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2327      }      }
2328      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2329        
2330      try {      try {
2331          m_nodeId.resize(getNumNodes());          m_nodeId.resize(getNumNodes());
2332          m_dofId.resize(numDOF);          m_dofId.resize(numDOF);
# Line 2133  void Brick::populateSampleIds() Line 2334  void Brick::populateSampleIds()
2334      } catch (const length_error& le) {      } catch (const length_error& le) {
2335          throw RipleyException("The system does not have sufficient memory for a domain of this size.");          throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2336      }      }
2337        
2338      // populate face element counts      // populate face element counts
2339      //left      //left
2340      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2352  void Brick::populateSampleIds() Line 2553  void Brick::populateSampleIds()
2553      setTagMap("front", FRONT);      setTagMap("front", FRONT);
2554      setTagMap("back", BACK);      setTagMap("back", BACK);
2555      updateTagsInUse(FaceElements);      updateTagsInUse(FaceElements);
2556    
2557        populateDofMap();
2558  }  }
2559    
2560  //private  //private
2561  vector<IndexVector> Brick::getConnections() const  vector<IndexVector> Brick::getConnections() const
2562  {  {
2563        // returns a vector v of size numDOF where v[i] is a vector with indices
2564        // of DOFs connected to i (up to 27 in 3D)
2565      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2566      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2567      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2568      const dim_t M = nDOF0*nDOF1*nDOF2;      const dim_t M = nDOF0*nDOF1*nDOF2;
2569      vector<IndexVector> indices(M);      vector<IndexVector> indices(M);
2570        
2571  #pragma omp parallel for  #pragma omp parallel for
2572      for (index_t i=0; i < M; i++) {      for (index_t i=0; i < M; i++) {
2573          const index_t x = i % nDOF0;          const index_t x = i % nDOF0;
# Line 2385  vector<IndexVector> Brick::getConnection Line 2590  vector<IndexVector> Brick::getConnection
2590  }  }
2591    
2592  //private  //private
2593  void Brick::createPattern()  void Brick::populateDofMap()
2594  {  {
2595      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2596      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
# Line 2406  void Brick::createPattern() Line 2611  void Brick::createPattern()
2611          }          }
2612      }      }
2613    
     // build list of shared components and neighbours by looping through  
     // all potential neighbouring ranks and checking if positions are  
     // within bounds  
2614      const dim_t numDOF=nDOF0*nDOF1*nDOF2;      const dim_t numDOF=nDOF0*nDOF1*nDOF2;
2615      RankVector neighbour;      RankVector neighbour;
2616      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2617      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2618      int numShared=0, expectedShared=0;;      dim_t numShared=0;
2619      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2620      const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2621      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
     for (int i2=-1; i2<2; i2++) {  
         for (int i1=-1; i1<2; i1++) {  
             for (int i0=-1; i0<2; i0++) {  
                 // skip this rank  
                 if (i0==0 && i1==0 && i2==0)  
                     continue;  
                 // location of neighbour rank  
                 const int nx=x+i0;  
                 const int ny=y+i1;  
                 const int nz=z+i2;  
                 if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {  
                     continue;  
                 }  
                 if (i0==0 && i1==0)  
                     expectedShared += nDOF0*nDOF1;  
                 else if (i0==0 && i2==0)  
                     expectedShared += nDOF0*nDOF2;  
                 else if (i1==0 && i2==0)  
                     expectedShared += nDOF1*nDOF2;  
                 else if (i0==0)  
                     expectedShared += nDOF0;  
                 else if (i1==0)  
                     expectedShared += nDOF1;  
                 else if (i2==0)  
                     expectedShared += nDOF2;  
                 else  
                     expectedShared++;  
             }  
         }  
     }  
       
     vector<IndexVector> colIndices(numDOF); // for the couple blocks  
     vector<IndexVector> rowIndices(expectedShared);  
2622    
2623        // build list of shared components and neighbours by looping through
2624        // all potential neighbouring ranks and checking if positions are
2625        // within bounds
2626      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2627          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2628              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
# Line 2474  void Brick::createPattern() Line 2646  void Brick::createPattern()
2646                              for (dim_t j=0; j<nDOF0; j++, numShared++) {                              for (dim_t j=0; j<nDOF0; j++, numShared++) {
2647                                  sendShared.push_back(firstDOF+j);                                  sendShared.push_back(firstDOF+j);
2648                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
                                 if (j>0) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);  
                                     if (i<nDOF1-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);  
                                 }  
                                 if (i>0)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);  
                                 doublyLink(colIndices, rowIndices, firstDOF+j, numShared);  
                                 if (i<nDOF1-1)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);  
                                 if (j<nDOF0-1) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);  
                                     if (i<nDOF1-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);  
                                 }  
2649                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2650                              }                              }
2651                          }                          }
# Line 2508  void Brick::createPattern() Line 2661  void Brick::createPattern()
2661                              for (dim_t j=0; j<nDOF0; j++, numShared++) {                              for (dim_t j=0; j<nDOF0; j++, numShared++) {
2662                                  sendShared.push_back(firstDOF+j);                                  sendShared.push_back(firstDOF+j);
2663                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
                                 if (j>0) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);  
                                     if (i<nDOF2-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);  
                                 }  
                                 if (i>0)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);  
                                 doublyLink(colIndices, rowIndices, firstDOF+j, numShared);  
                                 if (i<nDOF2-1)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);  
                                 if (j<nDOF0-1) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);  
                                     if (i<nDOF2-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);  
                                 }  
2664                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2665                              }                              }
2666                          }                          }
# Line 2542  void Brick::createPattern() Line 2676  void Brick::createPattern()
2676                              for (dim_t j=0; j<nDOF1; j++, numShared++) {                              for (dim_t j=0; j<nDOF1; j++, numShared++) {
2677                                  sendShared.push_back(firstDOF+j*nDOF0);                                  sendShared.push_back(firstDOF+j*nDOF0);
2678                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
                                 if (j>0) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);  
                                     if (i<nDOF2-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);  
                                 }  
                                 if (i>0)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);  
                                 doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);  
                                 if (i<nDOF2-1)  
                                     doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);  
                                 if (j<nDOF1-1) {  
                                     if (i>0)  
                                         doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);  
                                     doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);  
                                     if (i<nDOF2-1)  
                                         doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);  
                                 }  
2679                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2680                              }                              }
2681                          }                          }
# Line 2574  void Brick::createPattern() Line 2689  void Brick::createPattern()
2689                          for (dim_t i=0; i<nDOF0; i++, numShared++) {                          for (dim_t i=0; i<nDOF0; i++, numShared++) {
2690                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2691                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
                             if (i>0)  
                                 doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);  
                             doublyLink(colIndices, rowIndices, firstDOF+i, numShared);  
                             if (i<nDOF0-1)  
                                 doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);  
2692                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2693                          }                          }
2694                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2592  void Brick::createPattern() Line 2702  void Brick::createPattern()
2702                          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          for (dim_t i=0; i<nDOF1; i++, numShared++) {
2703                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2704                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
                             if (i>0)  
                                 doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);  
                             doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);  
                             if (i<nDOF1-1)  
                                 doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);  
2705                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2706                          }                          }
2707                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2610  void Brick::createPattern() Line 2715  void Brick::createPattern()
2715                          for (dim_t i=0; i<nDOF2; i++, numShared++) {                          for (dim_t i=0; i<nDOF2; i++, numShared++) {
2716                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2717                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
                             if (i>0)  
                                 doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);  
                             doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);  
                             if (i<nDOF2-1)  
                                 doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);  
2718                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2719                          }                          }
2720                      } else {                      } else {
2721                          // sharing a node                          // sharing a node
2722                          const int dof=(i0+1)/2*(nDOF0-1)                          const int dof = (i0+1)/2*(nDOF0-1)
2723                                        +(i1+1)/2*nDOF0*(nDOF1-1)                                         +(i1+1)/2*nDOF0*(nDOF1-1)
2724                                        +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                         +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2725                          const int node=(i0+1)/2*(m_NN[0]-1)                          const int node = (i0+1)/2*(m_NN[0]-1)
2726                                         +(i1+1)/2*m_NN[0]*(m_NN[1]-1)                                          +(i1+1)/2*m_NN[0]*(m_NN[1]-1)
2727                                         +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);                                          +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
2728                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2729                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2730                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2731                          doublyLink(colIndices, rowIndices, dof, numShared);                          m_dofMap[node] = numDOF+numShared;
                         m_dofMap[node]=numDOF+numShared;  
2732                          ++numShared;                          ++numShared;
2733                      }                      }
2734                  }                  }
# Line 2637  void Brick::createPattern() Line 2736  void Brick::createPattern()
2736          }          }
2737      }      }
2738    
 #pragma omp parallel for  
     for (int i = 0; i < numShared; i++) {  
         sort(rowIndices[i].begin(), rowIndices[i].end());  
     }  
   
   
2739      // TODO: paso::SharedComponents should take vectors to avoid this      // TODO: paso::SharedComponents should take vectors to avoid this
2740      Esys_MPI_rank* neighPtr = NULL;      Esys_MPI_rank* neighPtr = NULL;
2741      index_t* sendPtr = NULL;      index_t* sendPtr = NULL;
# Line 2661  void Brick::createPattern() Line 2754  void Brick::createPattern()
2754              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2755      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2756    
     // create main and couple blocks  
     paso::Pattern_ptr mainPattern = createPasoPattern(getConnections(), numDOF);  
     paso::Pattern_ptr colPattern = createPasoPattern(colIndices, numShared);  
     paso::Pattern_ptr rowPattern = createPasoPattern(rowIndices, numDOF);  
   
     // allocate paso distribution  
     paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,  
             const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));  
   
     // finally create the system matrix pattern  
     m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,  
             distribution, distribution, mainPattern, colPattern, rowPattern,  
             m_connector, m_connector));  
   
2757      // useful debug output      // useful debug output
2758      /*      /*
2759      cout << "--- rcv_shcomp ---" << endl;      cout << "--- rcv_shcomp ---" << endl;
# Line 2694  void Brick::createPattern() Line 2773  void Brick::createPattern()
2773      for (size_t i=0; i<m_dofMap.size(); i++) {      for (size_t i=0; i<m_dofMap.size(); i++) {
2774          cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;          cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
2775      }      }
     cout << "--- colIndices ---" << endl;  
     for (size_t i=0; i<colIndices.size(); i++) {  
         cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;  
     }  
     */  
   
     /*  
     cout << "--- main_pattern ---" << endl;  
     cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;  
     for (size_t i=0; i<mainPattern->numOutput+1; i++) {  
         cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;  
     }  
     for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {  
         cout << "index[" << i << "]=" << mainPattern->index[i] << endl;  
     }  
     */  
   
     /*  
     cout << "--- colCouple_pattern ---" << endl;  
     cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;  
     for (size_t i=0; i<colPattern->numOutput+1; i++) {  
         cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;  
     }  
     for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {  
         cout << "index[" << i << "]=" << colPattern->index[i] << endl;  
     }  
     */  
   
     /*  
     cout << "--- rowCouple_pattern ---" << endl;  
     cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;  
     for (size_t i=0; i<rowPattern->numOutput+1; i++) {  
         cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;  
     }  
     for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {  
         cout << "index[" << i << "]=" << rowPattern->index[i] << endl;  
     }  
2776      */      */
2777  }  }
2778    
# Line 3097  namespace Line 3139  namespace
3139              for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
3140              {              {
3141                  for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3142                  {                          {
3143                      arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));                      arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
3144                      total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];                          total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];
3145                  }                  }
3146              }              }
3147          }          }
3148          double invtotal=1/total;          double invtotal=1/total;
3149          for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)          for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3150          {          {
3151              arr[p]*=invtotal;              arr[p]*=invtotal;
3152          }          }
3153          return arr;          return arr;
3154      }      }
3155        
3156      // applies conv to source to get a point.      // applies conv to source to get a point.
3157      // (xp, yp) are the coords in the source matrix not the destination matrix      // (xp, yp) are the coords in the source matrix not the destination matrix
3158      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)      double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
# Line 3121  namespace Line 3163  namespace
3163          for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3164          {          {
3165              for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3166              {                  {
3167                  for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3168                  {                  {
3169                      result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];                      result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
# Line 3130  namespace Line 3172  namespace
3172          }          }
3173          // use this line for "pass-though" (return the centre point value)          // use this line for "pass-though" (return the centre point value)
3174  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3175          return result;                return result;
3176      }      }
3177  }  }
3178    
3179  /* This is a wrapper for filtered (and non-filtered) randoms  /* This is a wrapper for filtered (and non-filtered) randoms
3180   * For detailed doco see randomFillWorker   * For detailed doco see randomFillWorker
3181  */  */
3182  escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,  escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3183         const escript::FunctionSpace& what,         const escript::FunctionSpace& what,
3184         long seed, const boost::python::tuple& filter) const         long seed, const boost::python::tuple& filter) const
# Line 3161  A parameter radius  gives the size of th Line 3203  A parameter radius  gives the size of th
3203  A point on the left hand edge for example, will still require `radius` extra points to the left  A point on the left hand edge for example, will still require `radius` extra points to the left
3204  in order to complete the stencil.  in order to complete the stencil.
3205    
3206  All local calculation is done on an array called `src`, with  All local calculation is done on an array called `src`, with
3207  dimensions = ext[0] * ext[1] *ext[2].  dimensions = ext[0] * ext[1] *ext[2].
3208  Where ext[i]= internal[i]+2*radius.  Where ext[i]= internal[i]+2*radius.
3209    
3210  Now for MPI there is overlap to deal with. We need to share both the overlapping  Now for MPI there is overlap to deal with. We need to share both the overlapping
3211  values themselves but also the external region.  values themselves but also the external region.
3212    
3213  In a hypothetical 1-D case:  In a hypothetical 1-D case:
# Line 3184  need to be known. Line 3226  need to be known.
3226    
3227  pp123(4)56   23(4)567pp  pp123(4)56   23(4)567pp
3228    
3229  Now in our case, we wout set all the values 23456 on the left rank and send them to the  Now in our case, we wout set all the values 23456 on the left rank and send them to the
3230  right hand rank.  right hand rank.
3231    
3232  So the edges _may_ need to be shared at a distance `inset` from all boundaries.  So the edges _may_ need to be shared at a distance `inset` from all boundaries.
3233  inset=2*radius+1      inset=2*radius+1
3234  This is to ensure that values at distance `radius` from the shared/overlapped element  This is to ensure that values at distance `radius` from the shared/overlapped element
3235  that ripley has.  that ripley has.
3236  */  */
# Line 3196  escript::Data Brick::randomFillWorker(co Line 3238  escript::Data Brick::randomFillWorker(co
3238  {  {
3239      unsigned int radius=0;  // these are only used by gaussian      unsigned int radius=0;  // these are only used by gaussian
3240      double sigma=0.5;      double sigma=0.5;
3241        
3242      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
3243        
3244      if (len(filter)==0)      if (len(filter)==0)
3245      {      {
3246      // nothing special required here yet      // nothing special required here yet
3247      }      }
3248      else if (len(filter)==3)      else if (len(filter)==3)
3249      {      {
3250          boost::python::extract<string> ex(filter[0]);          boost::python::extract<string> ex(filter[0]);
3251          if (!ex.check() || (ex()!="gaussian"))          if (!ex.check() || (ex()!="gaussian"))
3252          {          {
3253              throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
3254          }          }
# Line 3221  escript::Data Brick::randomFillWorker(co Line 3263  escript::Data Brick::randomFillWorker(co
3263          if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
3264          {          {
3265              throw RipleyException("Sigma must be a positive floating point number.");              throw RipleyException("Sigma must be a positive floating point number.");
3266          }                      }
3267      }      }
3268      else      else
3269      {      {
# Line 3235  escript::Data Brick::randomFillWorker(co Line 3277  escript::Data Brick::randomFillWorker(co
3277      ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input      ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3278      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3279      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3280        
3281      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3282      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3283    
3284      if (2*radius>=internal[0]-4)      if (2*radius>=internal[0]-4)
# Line 3250  escript::Data Brick::randomFillWorker(co Line 3292  escript::Data Brick::randomFillWorker(co
3292      if (2*radius>=internal[2]-4)      if (2*radius>=internal[2]-4)
3293      {      {
3294          throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");          throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3295      }          }
3296      
3297      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3298      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);            esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);
3299        
3300  #ifdef ESYS_MPI  #ifdef ESYS_MPI
     
3301      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3302      {      {
3303      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
# Line 3267  escript::Data Brick::randomFillWorker(co Line 3308  escript::Data Brick::randomFillWorker(co
3308      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3309      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3310      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3311  #endif      #endif
3312    
3313  /*      /*
3314      // if we wanted to test a repeating pattern      // if we wanted to test a repeating pattern
3315      size_t basex=0;      size_t basex=0;
3316      size_t basey=0;      size_t basey=0;
3317      size_t basez=0;      size_t basez=0;
3318  #ifdef ESYS_MPI      #ifdef ESYS_MPI
3319      basex=X*m_gNE[0]/m_NX[0];      basex=X*m_gNE[0]/m_NX[0];
3320      basey=Y*m_gNE[1]/m_NX[1];      basey=Y*m_gNE[1]/m_NX[1];
3321      basez=Z*m_gNE[2]/m_NX[2];      basez=Z*m_gNE[2]/m_NX[2];
3322        cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;
3323  cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;      #endif
       
 #endif      
3324      esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);      esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3325  */  */
3326    
# Line 3290  cout << "basex=" << basex << " basey=" < Line 3329  cout << "basex=" << basex << " basey=" <
3329    
3330    
3331      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3332      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3333          // there is an overlap of two points both of which need to have "radius" points on either side.          // there is an overlap of two points both of which need to have "radius" points on either side.
3334        
3335      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3336      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;
3337      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3338        
3339      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);
3340        
3341      MPI_Request reqs[50];       // a non-tight upper bound on how many we need      MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3342      MPI_Status stats[50];      MPI_Status stats[50];
3343      short rused=0;      short rused=0;
3344        
3345      messvec incoms;      messvec incoms;
3346      messvec outcoms;      messvec outcoms;
3347        
3348      grid.generateInNeighbours(X, Y, Z ,incoms);      grid.generateInNeighbours(X, Y, Z ,incoms);
3349      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3350        
       
3351      block.copyAllToBuffer(src);      block.copyAllToBuffer(src);
3352        
3353      int comserr=0;          int comserr=0;
3354      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3355      {      {
3356          message& m=incoms[i];          message& m=incoms[i];
# Line 3324  cout << "basex=" << basex << " basey=" < Line 3362  cout << "basex=" << basex << " basey=" <
3362      {      {
3363          message& m=outcoms[i];          message& m=outcoms[i];
3364          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3365      }          }
3366        
3367      if (!comserr)      if (!comserr)
3368      {          {
3369          comserr=MPI_Waitall(rused, reqs, stats);          comserr=MPI_Waitall(rused, reqs, stats);
3370      }      }
3371    
# Line 3336  cout << "basex=" << basex << " basey=" < Line 3374  cout << "basex=" << basex << " basey=" <
3374      // Yes this is throwing an exception as a result of an MPI error.      // Yes this is throwing an exception as a result of an MPI error.
3375      // and no we don't inform the other ranks that we are doing this.      // and no we don't inform the other ranks that we are doing this.
3376      // however, we have no reason to believe coms work at this point anyway      // however, we have no reason to believe coms work at this point anyway
3377          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");
3378      }      }
3379        
3380      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3381        
3382  #endif      #endif
3383        
3384      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3385      {      {
         
3386          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3387          escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3388          // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3389          escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3390        
       
3391          // now we need to copy values over          // now we need to copy values over
3392          for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3393          {          {
3394              for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)
3395              {              {
3396                  for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3397                  {                  {
# Line 3365  cout << "basex=" << basex << " basey=" < Line 3401  cout << "basex=" << basex << " basey=" <
3401                      }                      }
3402                  }                  }
3403              }              }
3404          }            }
3405          delete[] src;          delete[] src;
3406          return resdat;                return resdat;
3407      }      }
3408      else        // filter enabled        else // filter enabled
3409      {      {
       
3410          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3411          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3412          // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
# Line 3380  cout << "basex=" << basex << " basey=" < Line 3415  cout << "basex=" << basex << " basey=" <
3415    
3416          for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3417          {          {
3418              for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)
3419              {              {
3420                  for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3421                  {                      {
3422                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3423                
3424                  }                  }
3425              }              }
3426          }          }
3427        
3428          delete[] convolution;          delete[] convolution;
3429          delete[] src;          delete[] src;
3430          return resdat;          return resdat;
       
3431      }      }
3432  }  }
3433    
3434    int Brick::findNode(const double *coords) const
3435    {
   
   
   
 int Brick::findNode(const double *coords) const {  
3436      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3437      //is the found element even owned by this rank      //is the found element even owned by this rank
3438      // (inside owned or shared elements but will map to an owned element)      // (inside owned or shared elements but will map to an owned element)
# Line 3419  int Brick::findNode(const double *coords Line 3449  int Brick::findNode(const double *coords
3449      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3450      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3451      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3452        
3453      //check if the point is even inside the domain      //check if the point is even inside the domain
3454      if (x < 0 || y < 0 || z < 0      if (x < 0 || y < 0 || z < 0
3455              || x > m_length[0] || y > m_length[1] || z > m_length[2])              || x > m_length[0] || y > m_length[1] || z > m_length[2])
3456          return NOT_MINE;          return NOT_MINE;
3457            
3458      // distance in elements      // distance in elements
3459      int ex = (int) floor(x / m_dx[0]);      int ex = (int) floor(x / m_dx[0]);
3460      int ey = (int) floor(y / m_dx[1]);      int ey = (int) floor(y / m_dx[1]);

Legend:
Removed from v.5118  
changed lines
  Added in v.5119

  ViewVC Help
Powered by ViewVC 1.1.26