/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

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

revision 4687 by jfenwick, Wed Feb 19 00:03:29 2014 UTC revision 4753 by sshaw, Mon Mar 17 02:39:44 2014 UTC
# Line 19  Line 19 
19  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
20  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
21  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
22    #include <ripley/LameAssembler3D.h>
23    #include <ripley/domainhelpers.h>
24  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 43  using esysUtils::FileWriter; Line 45  using esysUtils::FileWriter;
45    
46  namespace ripley {  namespace ripley {
47    
48    int indexOfMax(int a, int b, int c) {
49        if (a > b) {
50            if (c > a) {
51                return 2;
52            }
53            return 0;
54        } else if (b > c) {
55            return 1;
56        }
57        return 2;
58    }
59    
60  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
61               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
62               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
# Line 56  Brick::Brick(int n0, int n1, int n2, dou Line 70  Brick::Brick(int n0, int n1, int n2, dou
70          d2=1;          d2=1;
71      }      }
72      bool warn=false;      bool warn=false;
73      // if number of subdivisions is non-positive, try to subdivide by the same  
74      // ratio as the number of elements      std::vector<int> factors;
75      if (d0<=0 && d1<=0 && d2<=0) {      int ranks = m_mpiInfo->size;
76          warn=true;      int epr[3] = {n0,n1,n2};
77          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));      int d[3] = {d0,d1,d2};
78          d0=max(1, d0);      if (d0<=0 || d1<=0 || d2<=0) {
79          d1=max(1, (int)(d0*n1/(float)n0));          for (int i = 0; i < 3; i++) {
80          d2=m_mpiInfo->size/(d0*d1);              if (d[i] < 1) {
81          if (d0*d1*d2 != m_mpiInfo->size) {                  d[i] = 1;
82              // ratios not the same so leave "smallest" side undivided and try                  continue;
83              // dividing 2 sides only              }
84              if (n0>=n1) {              epr[i] = -1; // can no longer be max
85                  if (n1>=n2) {              if (ranks % d[i] != 0) {
86                      d0=d1=0;                  throw RipleyException("Invalid number of spatial subdivisions");
87                      d2=1;              }
88                  } else {              //remove
89                      d0=d2=0;              ranks /= d[i];
90                      d1=1;          }
91                  }          factorise(factors, ranks);
92              } else {          if (factors.size() != 0) {
93                  if (n0>=n2) {              warn = true;
94                      d0=d1=0;          }
95                      d2=1;      }
96                  } else {      while (factors.size() > 0) {
97                      d0=1;          int i = indexOfMax(epr[0],epr[1],epr[2]);
98                      d1=d2=0;          int f = factors.back();
99                  }          factors.pop_back();
100              }          d[i] *= f;
101          }          epr[i] /= f;
     }  
     if (d0<=0 && d1<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));  
         d1=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n1) {  
                 d0=0;  
                 d1=1;  
             } else {  
                 d0=1;  
                 d1=0;  
             }  
         }  
     } else if (d0<=0 && d2<=0) {  
         warn=true;  
         d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d0;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n0>n2) {  
                 d0=0;  
                 d2=1;  
             } else {  
                 d0=1;  
                 d2=0;  
             }  
         }  
     } else if (d1<=0 && d2<=0) {  
         warn=true;  
         d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));  
         d2=m_mpiInfo->size/d1;  
         if (d0*d1*d2 != m_mpiInfo->size) {  
             // ratios not the same so subdivide side with more elements only  
             if (n1>n2) {  
                 d1=0;  
                 d2=1;  
             } else {  
                 d1=1;  
                 d2=0;  
             }  
         }  
     }  
     if (d0<=0) {  
         // d1,d2 are preset, determine d0  
         d0=m_mpiInfo->size/(d1*d2);  
     } else if (d1<=0) {  
         // d0,d2 are preset, determine d1  
         d1=m_mpiInfo->size/(d0*d2);  
     } else if (d2<=0) {  
         // d0,d1 are preset, determine d2  
         d2=m_mpiInfo->size/(d0*d1);  
102      }      }
103        d0 = d[0]; d1 = d[1]; d2 = d[2];
104    
105      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
106      // among number of ranks      // among number of ranks
107      if (d0*d1*d2 != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
108          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
109        }
110      if (warn) {      if (warn) {
111          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
112              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;              << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
# Line 412  void Brick::readNcGrid(escript::Data& ou Line 374  void Brick::readNcGrid(escript::Data& ou
374  #endif  #endif
375  }  }
376    
377    #ifdef USE_BOOSTIO
378    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
379                               const ReaderParameters& params) const
380    {
381        // the mapping is not universally correct but should work on our
382        // supported platforms
383        switch (params.dataType) {
384            case DATATYPE_INT32:
385                readBinaryGridZippedImpl<int>(out, filename, params);
386                break;
387            case DATATYPE_FLOAT32:
388                readBinaryGridZippedImpl<float>(out, filename, params);
389                break;
390            case DATATYPE_FLOAT64:
391                readBinaryGridZippedImpl<double>(out, filename, params);
392                break;
393            default:
394                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
395        }
396    }
397    #endif
398    
399  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
400                             const ReaderParameters& params) const                             const ReaderParameters& params) const
401  {  {
# Line 548  void Brick::readBinaryGridImpl(escript:: Line 532  void Brick::readBinaryGridImpl(escript::
532      f.close();      f.close();
533  }  }
534    
535    #ifdef USE_BOOSTIO
536    template<typename ValueType>
537    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
538                                   const ReaderParameters& params) const
539    {
540        // check destination function space
541        int myN0, myN1, myN2;
542        if (out.getFunctionSpace().getTypeCode() == Nodes) {
543            myN0 = m_NN[0];
544            myN1 = m_NN[1];
545            myN2 = m_NN[2];
546        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
547                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
548            myN0 = m_NE[0];
549            myN1 = m_NE[1];
550            myN2 = m_NE[2];
551        } else
552            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
553    
554        if (params.first.size() != 3)
555            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
556    
557        if (params.numValues.size() != 3)
558            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
559    
560        if (params.multiplier.size() != 3)
561            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
562        for (size_t i=0; i<params.multiplier.size(); i++)
563            if (params.multiplier[i]<1)
564                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
565    
566        // check file existence and size
567        ifstream f(filename.c_str(), ifstream::binary);
568        if (f.fail()) {
569            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
570        }
571        f.seekg(0, ios::end);
572        const int numComp = out.getDataPointSize();
573        int filesize = f.tellg();
574        f.seekg(0, ios::beg);
575        std::vector<char> compressed(filesize);
576        f.read((char*)&compressed[0], filesize);
577        f.close();
578        std::vector<char> decompressed = unzip(compressed);
579        filesize = decompressed.size();
580        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
581        if (filesize < reqsize) {
582            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
583        }
584    
585        // check if this rank contributes anything
586        if (params.first[0] >= m_offset[0]+myN0 ||
587                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
588                params.first[1] >= m_offset[1]+myN1 ||
589                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
590                params.first[2] >= m_offset[2]+myN2 ||
591                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
592            return;
593        }
594    
595        // now determine how much this rank has to write
596    
597        // first coordinates in data object to write to
598        const int first0 = max(0, params.first[0]-m_offset[0]);
599        const int first1 = max(0, params.first[1]-m_offset[1]);
600        const int first2 = max(0, params.first[2]-m_offset[2]);
601        // indices to first value in file
602        const int idx0 = max(0, m_offset[0]-params.first[0]);
603        const int idx1 = max(0, m_offset[1]-params.first[1]);
604        const int idx2 = max(0, m_offset[2]-params.first[2]);
605        // number of values to read
606        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
607        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
608        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
609    
610        out.requireWrite();
611        vector<ValueType> values(num0*numComp);
612        const int dpp = out.getNumDataPointsPerSample();
613    
614        for (int z=0; z<num2; z++) {
615            for (int y=0; y<num1; y++) {
616                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
617                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
618                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
619                
620                for (int x=0; x<num0; x++) {
621                    const int baseIndex = first0+x*params.multiplier[0]
622                                         +(first1+y*params.multiplier[1])*myN0
623                                         +(first2+z*params.multiplier[2])*myN0*myN1;
624                    for (int m2=0; m2<params.multiplier[2]; m2++) {
625                        for (int m1=0; m1<params.multiplier[1]; m1++) {
626                            for (int m0=0; m0<params.multiplier[0]; m0++) {
627                                const int dataIndex = baseIndex+m0
628                                               +m1*myN0
629                                               +m2*myN0*myN1;
630                                double* dest = out.getSampleDataRW(dataIndex);
631                                for (int c=0; c<numComp; c++) {
632                                    ValueType val = values[x*numComp+c];
633    
634                                    if (params.byteOrder != BYTEORDER_NATIVE) {
635                                        char* cval = reinterpret_cast<char*>(&val);
636                                        // this will alter val!!
637                                        byte_swap32(cval);
638                                    }
639                                    if (!std::isnan(val)) {
640                                        for (int q=0; q<dpp; q++) {
641                                            *dest++ = static_cast<double>(val);
642                                        }
643                                    }
644                                }
645                            }
646                        }
647                    }
648                }
649            }
650        }
651    }
652    #endif
653    
654  void Brick::writeBinaryGrid(const escript::Data& in, string filename,  void Brick::writeBinaryGrid(const escript::Data& in, string filename,
655                              int byteOrder, int dataType) const                              int byteOrder, int dataType) const
656  {  {
# Line 2231  void Brick::createPattern() Line 2334  void Brick::createPattern()
2334      RankVector neighbour;      RankVector neighbour;
2335      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2336      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2337      int numShared=0;      int numShared=0, expectedShared=0;;
2338      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2339      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];
2340      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
# Line 2245  void Brick::createPattern() Line 2348  void Brick::createPattern()
2348                  const int nx=x+i0;                  const int nx=x+i0;
2349                  const int ny=y+i1;                  const int ny=y+i1;
2350                  const int nz=z+i2;                  const int nz=z+i2;
2351                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2352                        continue;
2353                    }
2354                    if (i0==0 && i1==0)
2355                        expectedShared += nDOF0*nDOF1;
2356                    else if (i0==0 && i2==0)
2357                        expectedShared += nDOF0*nDOF2;
2358                    else if (i1==0 && i2==0)
2359                        expectedShared += nDOF1*nDOF2;
2360                    else if (i0==0)
2361                        expectedShared += nDOF0;
2362                    else if (i1==0)
2363                        expectedShared += nDOF1;
2364                    else if (i2==0)
2365                        expectedShared += nDOF2;
2366                    else
2367                        expectedShared++;
2368                }
2369            }
2370        }
2371        
2372        vector<IndexVector> rowIndices(expectedShared);
2373        
2374        for (int i2=-1; i2<2; i2++) {
2375            for (int i1=-1; i1<2; i1++) {
2376                for (int i0=-1; i0<2; i0++) {
2377                    // skip this rank
2378                    if (i0==0 && i1==0 && i2==0)
2379                        continue;
2380                    // location of neighbour rank
2381                    const int nx=x+i0;
2382                    const int ny=y+i1;
2383                    const int nz=z+i2;
2384                  if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {                  if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {
2385                      neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);                      neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);
2386                      if (i0==0 && i1==0) {                      if (i0==0 && i1==0) {
# Line 2260  void Brick::createPattern() Line 2396  void Brick::createPattern()
2396                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2397                                  if (j>0) {                                  if (j>0) {
2398                                      if (i>0)                                      if (i>0)
2399                                          colIndices[firstDOF+j-1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2400                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2401                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2402                                          colIndices[firstDOF+j-1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2403                                  }                                  }
2404                                  if (i>0)                                  if (i>0)
2405                                      colIndices[firstDOF+j-nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2406                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2407                                  if (i<nDOF1-1)                                  if (i<nDOF1-1)
2408                                      colIndices[firstDOF+j+nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2409                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2410                                      if (i>0)                                      if (i>0)
2411                                          colIndices[firstDOF+j+1-nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2412                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2413                                      if (i<nDOF1-1)                                      if (i<nDOF1-1)
2414                                          colIndices[firstDOF+j+1+nDOF0].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2415                                  }                                  }
2416                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2417                              }                              }
# Line 2294  void Brick::createPattern() Line 2430  void Brick::createPattern()
2430                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2431                                  if (j>0) {                                  if (j>0) {
2432                                      if (i>0)                                      if (i>0)
2433                                          colIndices[firstDOF+j-1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2434                                      colIndices[firstDOF+j-1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2435                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2436                                          colIndices[firstDOF+j-1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2437                                  }                                  }
2438                                  if (i>0)                                  if (i>0)
2439                                      colIndices[firstDOF+j-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2440                                  colIndices[firstDOF+j].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2441                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2442                                      colIndices[firstDOF+j+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2443                                  if (j<nDOF0-1) {                                  if (j<nDOF0-1) {
2444                                      if (i>0)                                      if (i>0)
2445                                          colIndices[firstDOF+j+1-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2446                                      colIndices[firstDOF+j+1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2447                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2448                                          colIndices[firstDOF+j+1+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2449                                  }                                  }
2450                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2451                              }                              }
# Line 2328  void Brick::createPattern() Line 2464  void Brick::createPattern()
2464                                  recvShared.push_back(numDOF+numShared);                                  recvShared.push_back(numDOF+numShared);
2465                                  if (j>0) {                                  if (j>0) {
2466                                      if (i>0)                                      if (i>0)
2467                                          colIndices[firstDOF+(j-1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2468                                      colIndices[firstDOF+(j-1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2469                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2470                                          colIndices[firstDOF+(j-1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2471                                  }                                  }
2472                                  if (i>0)                                  if (i>0)
2473                                      colIndices[firstDOF+j*nDOF0-nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2474                                  colIndices[firstDOF+j*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2475                                  if (i<nDOF2-1)                                  if (i<nDOF2-1)
2476                                      colIndices[firstDOF+j*nDOF0+nDOF0*nDOF1].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2477                                  if (j<nDOF1-1) {                                  if (j<nDOF1-1) {
2478                                      if (i>0)                                      if (i>0)
2479                                          colIndices[firstDOF+(j+1)*nDOF0-nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2480                                      colIndices[firstDOF+(j+1)*nDOF0].push_back(numShared);                                      doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2481                                      if (i<nDOF2-1)                                      if (i<nDOF2-1)
2482                                          colIndices[firstDOF+(j+1)*nDOF0+nDOF0*nDOF1].push_back(numShared);                                          doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2483                                  }                                  }
2484                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2485                              }                              }
# Line 2359  void Brick::createPattern() Line 2495  void Brick::createPattern()
2495                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2496                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2497                              if (i>0)                              if (i>0)
2498                                  colIndices[firstDOF+i-1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2499                              colIndices[firstDOF+i].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2500                              if (i<nDOF0-1)                              if (i<nDOF0-1)
2501                                  colIndices[firstDOF+i+1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2502                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2503                          }                          }
2504                      } else if (i1==0) {                      } else if (i1==0) {
# Line 2377  void Brick::createPattern() Line 2513  void Brick::createPattern()
2513                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2514                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2515                              if (i>0)                              if (i>0)
2516                                  colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2517                              colIndices[firstDOF+i*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2518                              if (i<nDOF1-1)                              if (i<nDOF1-1)
2519                                  colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2520                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2521                          }                          }
2522                      } else if (i2==0) {                      } else if (i2==0) {
# Line 2395  void Brick::createPattern() Line 2531  void Brick::createPattern()
2531                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2532                              recvShared.push_back(numDOF+numShared);                              recvShared.push_back(numDOF+numShared);
2533                              if (i>0)                              if (i>0)
2534                                  colIndices[firstDOF+(i-1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2535                              colIndices[firstDOF+i*nDOF0*nDOF1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2536                              if (i<nDOF2-1)                              if (i<nDOF2-1)
2537                                  colIndices[firstDOF+(i+1)*nDOF0*nDOF1].push_back(numShared);                                  doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2538                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2539                          }                          }
2540                      } else {                      } else {
# Line 2412  void Brick::createPattern() Line 2548  void Brick::createPattern()
2548                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2549                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2550                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2551                          colIndices[dof].push_back(numShared);                          doublyLink(colIndices, rowIndices, dof, numShared);
2552                          m_dofMap[node]=numDOF+numShared;                          m_dofMap[node]=numDOF+numShared;
2553                          ++numShared;                          ++numShared;
2554                      }                      }
# Line 2421  void Brick::createPattern() Line 2557  void Brick::createPattern()
2557          }          }
2558      }      }
2559    
2560    #pragma omp parallel for
2561        for (int i = 0; i < numShared; i++) {
2562            std::sort(rowIndices[i].begin(), rowIndices[i].end());
2563        }
2564    
2565      // create connector      // create connector
2566      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
2567              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 2435  void Brick::createPattern() Line 2576  void Brick::createPattern()
2576      // create main and couple blocks      // create main and couple blocks
2577      Paso_Pattern *mainPattern = createMainPattern();      Paso_Pattern *mainPattern = createMainPattern();
2578      Paso_Pattern *colPattern, *rowPattern;      Paso_Pattern *colPattern, *rowPattern;
2579      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
2580    
2581      // allocate paso distribution      // allocate paso distribution
2582      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
# Line 2867  namespace Line 3008  namespace
3008      {      {
3009          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
3010          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
3011      double total=0;          double total=0;
3012      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
3013      for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z)
3014      {          {
3015          for (int y=-r;y<=r;++y)              for (int y=-r;y<=r;++y)
3016          {              {
3017          for (int x=-r;x<=r;++x)                  for (int x=-r;x<=r;++x)
3018          {                          {        
3019              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));
3020              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)];    
3021          }                  }
3022          }              }
3023      }          }
3024      double invtotal=1/total;          double invtotal=1/total;
3025      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)
3026      {          {
3027          arr[p]*=invtotal;              arr[p]*=invtotal;
3028      }          }
3029      return arr;          return arr;
3030      }      }
3031            
3032      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 2893  namespace Line 3034  namespace
3034      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)
3035      {      {
3036          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3037      size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3038      double result=0;          double result=0;
3039      for (int z=0;z<2*radius+1;++z)          for (int z=0;z<2*radius+1;++z)
3040      {          {
3041          for (int y=0;y<2*radius+1;++y)              for (int y=0;y<2*radius+1;++y)
3042          {                  {    
3043          for (int x=0;x<2*radius+1;++x)                  for (int x=0;x<2*radius+1;++x)
3044          {                  {
3045              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];
3046          }                  }
3047          }              }
3048      }          }
3049      // use this line for "pass-though" (return the centre point value)          // use this line for "pass-though" (return the centre point value)
3050  //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3051      return result;                return result;      
3052      }      }
3053  }  }
3054    
# Line 2921  escript::Data Brick::randomFill(const es Line 3062  escript::Data Brick::randomFill(const es
3062      int numvals=escript::DataTypes::noValues(shape);      int numvals=escript::DataTypes::noValues(shape);
3063      if (len(filter)>0 && (numvals!=1))      if (len(filter)>0 && (numvals!=1))
3064      {      {
3065      throw RipleyException("Ripley only supports filters for scalar data.");          throw RipleyException("Ripley only supports filters for scalar data.");
3066      }      }
3067      escript::Data res=randomFillWorker(shape, seed, filter);      escript::Data res=randomFillWorker(shape, seed, filter);
3068      if (res.getFunctionSpace()!=what)      if (res.getFunctionSpace()!=what)
3069      {      {
3070      escript::Data r=escript::Data(res, what);          escript::Data r=escript::Data(res, what);
3071      return r;          return r;
3072      }      }
3073      return res;      return res;
3074  }  }
# Line 2976  escript::Data Brick::randomFillWorker(co Line 3117  escript::Data Brick::randomFillWorker(co
3117          throw RipleyException("Brick must be 3D.");          throw RipleyException("Brick must be 3D.");
3118      }      }
3119            
3120      unsigned int radius=0;  // these are only used by gaussian      unsigned int radius=0;  // these are only used by gaussian
3121      double sigma=0.5;      double sigma=0.5;
3122            
3123      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
3124            
3125      if (len(filter)==0)      if (len(filter)==0)
3126      {      {
3127      // nothing special required here yet      // nothing special required here yet
3128      }      }
3129      else if (len(filter)==3)      else if (len(filter)==3)
3130      {      {
3131      boost::python::extract<string> ex(filter[0]);          boost::python::extract<string> ex(filter[0]);
3132      if (!ex.check() || (ex()!="gaussian"))          if (!ex.check() || (ex()!="gaussian"))
3133      {          {
3134          throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
3135      }          }
3136      boost::python::extract<unsigned int> ex1(filter[1]);          boost::python::extract<unsigned int> ex1(filter[1]);
3137      if (!ex1.check())          if (!ex1.check())
3138      {          {
3139          throw RipleyException("Radius of gaussian filter must be a positive integer.");              throw RipleyException("Radius of gaussian filter must be a positive integer.");
3140      }          }
3141      radius=ex1();          radius=ex1();
3142      sigma=0.5;          sigma=0.5;
3143      boost::python::extract<double> ex2(filter[2]);          boost::python::extract<double> ex2(filter[2]);
3144      if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2())<=0)
3145      {          {
3146          throw RipleyException("Sigma must be a postive floating point number.");              throw RipleyException("Sigma must be a postive floating point number.");
3147      }                      }            
3148      }      }
3149      else      else
3150      {      {
# Line 3014  escript::Data Brick::randomFillWorker(co Line 3155  escript::Data Brick::randomFillWorker(co
3155    
3156            
3157      size_t internal[3];      size_t internal[3];
3158      internal[0]=m_NE[0]+1;  // number of points in the internal region      internal[0]=m_NE[0]+1;  // number of points in the internal region
3159      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of
3160      internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of      internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of
3161      size_t ext[3];      size_t ext[3];
3162      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=internal[0]+2*radius;    // includes points we need as input
3163      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=internal[1]+2*radius;    // for smoothing
3164      ext[2]=internal[2]+2*radius;    // for smoothing      ext[2]=internal[2]+2*radius;    // for smoothing
3165            
3166      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3167      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
# Line 3045  escript::Data Brick::randomFillWorker(co Line 3186  escript::Data Brick::randomFillWorker(co
3186        
3187      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))      if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3188      {      {
3189      // since the dimensions are equal for all ranks, this exception      // since the dimensions are equal for all ranks, this exception
3190      // will be thrown on all ranks      // will be thrown on all ranks
3191      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");      throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
3192    
3193      }      }
3194      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
# Line 3063  escript::Data Brick::randomFillWorker(co Line 3204  escript::Data Brick::randomFillWorker(co
3204  #ifdef ESYS_MPI      #ifdef ESYS_MPI    
3205      basex=X*m_gNE[0]/m_NX[0];      basex=X*m_gNE[0]/m_NX[0];
3206      basey=Y*m_gNE[1]/m_NX[1];      basey=Y*m_gNE[1]/m_NX[1];
3207      basez=Z*m_gNE[2]/m_NX[2];          basez=Z*m_gNE[2]/m_NX[2];
 #endif      
     if (seed==0)  
     {  
     seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception  
     }  
     esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, seed, basex, basey, basez);  
 */  
3208            
3209    cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;    
3210            
3211  /*  #endif    
3212  cout << "Pattern:\n";          esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3213  for (int i=0;i<ext[0]*ext[1]*ext[2];)  */
3214  {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }*/  
       
     
3215  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3216    
3217    
3218    
3219      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);
3220      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
3221          // 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.
3222            
3223      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
3224      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;  
3225      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3226            
3227      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);    
3228            
3229      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
3230      MPI_Status stats[50];      MPI_Status stats[50];
3231      short rused=0;      short rused=0;
3232            
# Line 3112  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3239  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3239            
3240      block.copyAllToBuffer(src);      block.copyAllToBuffer(src);
3241            
       
3242      int comserr=0;          int comserr=0;    
3243      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<incoms.size();++i)
3244      {      {
3245      message& m=incoms[i];          message& m=incoms[i];
3246      comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3247      block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3248      }      }
3249    
3250      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3251      {      {
3252      message& m=outcoms[i];          message& m=outcoms[i];
3253      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++));
3254      }          }    
3255            
3256      if (!comserr)      if (!comserr)
# Line 3134  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3260  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3260    
3261      if (comserr)      if (comserr)
3262      {      {
3263      // 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.
3264      // 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.
3265      // 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
3266          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
3267      }      }
3268            
3269      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3270            
       
       
   
3271  #endif      #endif    
3272            
3273  /*          if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
 cout << "Pattern (post transfer):\n";      
 for (int i=0;i<ext[0]*ext[1]*ext[2];)  
 {  
     cout << src[i++] << " ";  
     if (i%ext[0]==0)  
       cout << "\n";  
     if (i%(ext[0]*ext[1])==0)  
       cout << "\n";  
 }   */  
       
     if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe  
3274      {      {
3275                
3276      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3277      escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3278      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3279      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3280            
3281            
3282      // now we need to copy values over          // now we need to copy values over
3283      for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z)
3284      {          {
3285          for (size_t y=0;y<(internal[1]);++y)                  for (size_t y=0;y<(internal[1]);++y)    
3286          {              {
3287          for (size_t x=0;x<(internal[0]);++x)                  for (size_t x=0;x<(internal[0]);++x)
3288          {                  {
3289              for (unsigned int i=0;i<numvals;++i)                      for (unsigned int i=0;i<numvals;++i)
3290              {                      {
3291                  dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];                          dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3292              }                      }
3293          }                  }
3294          }              }
3295      }            }  
3296      delete[] src;          delete[] src;
3297      return resdat;                return resdat;      
3298      }      }
3299      else        // filter enabled        else        // filter enabled  
3300      {      {
3301            
3302      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3303      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3304      // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3305      escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3306      double* convolution=get3DGauss(radius, sigma);          double* convolution=get3DGauss(radius, sigma);
3307    
3308  // cout << "Convolution matrix\n";          for (size_t z=0;z<(internal[2]);++z)
3309  // size_t di=(radius*2+1);          {
3310  // double ctot=0;              for (size_t y=0;y<(internal[1]);++y)    
3311  // for (int i=0;i<di*di*di;++i)              {
3312  // {                  for (size_t x=0;x<(internal[0]);++x)
3313  //     cout << convolution[i] << " ";                  {    
3314  //     ctot+=convolution[i];                      dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3315  //     if ((i+1)%di==0)              
3316  //     {                  }
3317  //  cout << "\n";              }
3318  //     }          }
3319  //     if ((i+1)%(di*di)==0)      
3320  //     {          delete[] convolution;
3321  //  cout << "\n";          delete[] src;
3322  //     }          return resdat;
 // }  
 //  
 // cout << "Sum of matrix is = " << ctot << endl;  
   
     for (size_t z=0;z<(internal[2]);++z)  
     {  
         for (size_t y=0;y<(internal[1]);++y)      
         {  
         for (size_t x=0;x<(internal[0]);++x)  
         {      
             dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);  
               
         }  
         }  
     }  
       
 // cout << "\nResulting matrix:\n";  
 //     for (size_t z=0;z<(internal[2]);++z)  
 //     {  
 //  for (size_t y=0;y<(internal[1]);++y)      
 //  {  
 //      for (size_t x=0;x<(internal[0]);++x)  
 //      {      
 //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";  
 //      }  
 //      cout << endl;  
 //  }  
 //  cout << endl;  
 //     }  
   
       
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3323            
3324      }      }
3325  }  }
# Line 3303  int Brick::findNode(const double *coords Line 3381  int Brick::findNode(const double *coords
3381  void Brick::setAssembler(std::string type, std::map<std::string,  void Brick::setAssembler(std::string type, std::map<std::string,
3382          escript::Data> constants) {          escript::Data> constants) {
3383      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
3384            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3385                throw RipleyException("Domain already using a different custom assembler");
3386            assembler_type = WAVE_ASSEMBLER;
3387          delete assembler;          delete assembler;
3388          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);          assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);
3389        } else if (type.compare("LameAssembler") == 0) {
3390            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
3391                throw RipleyException("Domain already using a different custom assembler");
3392            assembler_type = LAME_ASSEMBLER;
3393            delete assembler;
3394            assembler = new LameAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
3395      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3396          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Brick does not support the"
3397                                  " requested assembler");                                  " requested assembler");
3398      }      }
3399  }  }

Legend:
Removed from v.4687  
changed lines
  Added in v.4753

  ViewVC Help
Powered by ViewVC 1.1.26