/[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 5084 by caltinay, Sun Jun 29 23:29:51 2014 UTC revision 5122 by caltinay, Thu Aug 21 04:32:39 2014 UTC
# Line 16  Line 16 
16    
17  #include <ripley/Brick.h>  #include <ripley/Brick.h>
18  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
 #include <ripley/WaveAssembler3D.h>  
19  #include <ripley/LameAssembler3D.h>  #include <ripley/LameAssembler3D.h>
20    #include <ripley/WaveAssembler3D.h>
21  #include <ripley/blocktools.h>  #include <ripley/blocktools.h>
22  #include <ripley/domainhelpers.h>  #include <ripley/domainhelpers.h>
23  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
# Line 40  Line 40 
40  #include <iomanip>  #include <iomanip>
41  #include <limits>  #include <limits>
42    
43    namespace bp = boost::python;
44  using namespace std;  using namespace std;
45  using esysUtils::FileWriter;  using esysUtils::FileWriter;
46    using escript::AbstractSystemMatrix;
47    
48  namespace ripley {  namespace ripley {
49    
50  int indexOfMax(int a, int b, int c) {  inline int indexOfMax(dim_t a, dim_t b, dim_t c)
51    {
52      if (a > b) {      if (a > b) {
53          if (c > a) {          if (c > a) {
54              return 2;              return 2;
# Line 58  int indexOfMax(int a, int b, int c) { Line 60  int indexOfMax(int a, int b, int c) {
60      return 2;      return 2;
61  }  }
62    
63  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(dim_t n0, dim_t n1, dim_t n2, double x0, double y0, double z0,
64               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
65               const std::vector<double>& points, const std::vector<int>& tags,               const vector<double>& points, const vector<int>& tags,
66               const simap_t& tagnamestonums,               const simap_t& tagnamestonums,
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) > std::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.");
74    
# Line 82  Brick::Brick(int n0, int n1, int n2, dou Line 84  Brick::Brick(int n0, int n1, int n2, dou
84      }      }
85      bool warn=false;      bool warn=false;
86    
87      std::vector<int> factors;      vector<int> factors;
88      int ranks = m_mpiInfo->size;      int ranks = m_mpiInfo->size;
89      int epr[3] = {n0,n1,n2};      dim_t epr[3] = {n0,n1,n2};
90      int d[3] = {d0,d1,d2};      int d[3] = {d0,d1,d2};
91      if (d0<=0 || d1<=0 || d2<=0) {      if (d0<=0 || d1<=0 || d2<=0) {
92          for (int i = 0; i < 3; i++) {          for (int i = 0; i < 3; i++) {
# Line 96  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 131  Brick::Brick(int n0, int n1, int n2, dou Line 133  Brick::Brick(int n0, int n1, int n2, dou
133      m_dx[2] = l2/n2;      m_dx[2] = l2/n2;
134    
135      if ((n0+1)%d0 > 0) {      if ((n0+1)%d0 > 0) {
136          n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;          n0=(dim_t)round((float)(n0+1)/d0+0.5)*d0-1;
137          l0=m_dx[0]*n0;          l0=m_dx[0]*n0;
138          cout << "Warning: Adjusted number of elements and length. N0="          cout << "Warning: Adjusted number of elements and length. N0="
139              << n0 << ", l0=" << l0 << endl;              << n0 << ", l0=" << l0 << endl;
140      }      }
141      if ((n1+1)%d1 > 0) {      if ((n1+1)%d1 > 0) {
142          n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;          n1=(dim_t)round((float)(n1+1)/d1+0.5)*d1-1;
143          l1=m_dx[1]*n1;          l1=m_dx[1]*n1;
144          cout << "Warning: Adjusted number of elements and length. N1="          cout << "Warning: Adjusted number of elements and length. N1="
145              << n1 << ", l1=" << l1 << endl;              << n1 << ", l1=" << l1 << endl;
146      }      }
147      if ((n2+1)%d2 > 0) {      if ((n2+1)%d2 > 0) {
148          n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;          n2=(dim_t)round((float)(n2+1)/d2+0.5)*d2-1;
149          l2=m_dx[2]*n2;          l2=m_dx[2]*n2;
150          cout << "Warning: Adjusted number of elements and length. N2="          cout << "Warning: Adjusted number of elements and length. N2="
151              << n2 << ", l2=" << l2 << endl;              << n2 << ", l2=" << l2 << endl;
# Line 201  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 (map<string, int>::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);
210      }      }
211      addPoints(tags.size(), &points[0], &tags[0]);      addPoints(points, tags);
212  }  }
213    
   
214  Brick::~Brick()  Brick::~Brick()
215  {  {
216  }  }
# Line 239  void Brick::readNcGrid(escript::Data& ou Line 239  void Brick::readNcGrid(escript::Data& ou
239  {  {
240  #ifdef USE_NETCDF  #ifdef USE_NETCDF
241      // check destination function space      // check destination function space
242      int myN0, myN1, myN2;      dim_t myN0, myN1, myN2;
243      if (out.getFunctionSpace().getTypeCode() == Nodes) {      if (out.getFunctionSpace().getTypeCode() == Nodes) {
244          myN0 = m_NN[0];          myN0 = m_NN[0];
245          myN1 = m_NN[1];          myN1 = m_NN[1];
# Line 304  void Brick::readNcGrid(escript::Data& ou Line 304  void Brick::readNcGrid(escript::Data& ou
304      // now determine how much this rank has to write      // now determine how much this rank has to write
305    
306      // first coordinates in data object to write to      // first coordinates in data object to write to
307      const int first0 = max(0, params.first[0]-m_offset[0]);      const dim_t first0 = max(0, params.first[0]-m_offset[0]);
308      const int first1 = max(0, params.first[1]-m_offset[1]);      const dim_t first1 = max(0, params.first[1]-m_offset[1]);
309      const int first2 = max(0, params.first[2]-m_offset[2]);      const dim_t first2 = max(0, params.first[2]-m_offset[2]);
310      // indices to first value in file (not accounting for reverse yet)      // indices to first value in file (not accounting for reverse yet)
311      int idx0 = max(0, m_offset[0]-params.first[0]);      dim_t idx0 = max(0, m_offset[0]-params.first[0]);
312      int idx1 = max(0, m_offset[1]-params.first[1]);      dim_t idx1 = max(0, m_offset[1]-params.first[1]);
313      int idx2 = max(0, m_offset[2]-params.first[2]);      dim_t idx2 = max(0, m_offset[2]-params.first[2]);
314      // number of values to read      // number of values to read
315      const int num0 = min(params.numValues[0]-idx0, myN0-first0);      const dim_t num0 = min(params.numValues[0]-idx0, myN0-first0);
316      const int num1 = min(params.numValues[1]-idx1, myN1-first1);      const dim_t num1 = min(params.numValues[1]-idx1, myN1-first1);
317      const int num2 = min(params.numValues[2]-idx2, myN2-first2);      const dim_t num2 = min(params.numValues[2]-idx2, myN2-first2);
318    
319      // make sure we read the right block if going backwards through file      // make sure we read the right block if going backwards through file
320      if (params.reverse[0])      if (params.reverse[0])
# Line 341  void Brick::readNcGrid(escript::Data& ou Line 341  void Brick::readNcGrid(escript::Data& ou
341      out.requireWrite();      out.requireWrite();
342    
343      // helpers for reversing      // helpers for reversing
344      const int x0 = (params.reverse[0] ? num0-1 : 0);      const dim_t x0 = (params.reverse[0] ? num0-1 : 0);
345      const int x_mult = (params.reverse[0] ? -1 : 1);      const int x_mult = (params.reverse[0] ? -1 : 1);
346      const int y0 = (params.reverse[1] ? num1-1 : 0);      const dim_t y0 = (params.reverse[1] ? num1-1 : 0);
347      const int y_mult = (params.reverse[1] ? -1 : 1);      const int y_mult = (params.reverse[1] ? -1 : 1);
348      const int z0 = (params.reverse[2] ? num2-1 : 0);      const dim_t z0 = (params.reverse[2] ? num2-1 : 0);
349      const int z_mult = (params.reverse[2] ? -1 : 1);      const int z_mult = (params.reverse[2] ? -1 : 1);
350    
351      for (index_t z=0; z<num2; z++) {      for (index_t z=0; z<num2; z++) {
352          for (index_t y=0; y<num1; y++) {          for (index_t y=0; y<num1; y++) {
353  #pragma omp parallel for  #pragma omp parallel for
354              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
355                  const int baseIndex = first0+x*params.multiplier[0]                  const dim_t baseIndex = first0+x*params.multiplier[0]
356                                       +(first1+y*params.multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
357                                       +(first2+z*params.multiplier[2])*myN0*myN1;                                       +(first2+z*params.multiplier[2])*myN0*myN1;
358                  const int srcIndex=(z0+z_mult*z)*num1*num0                  const dim_t srcIndex=(z0+z_mult*z)*num1*num0
359                                    +(y0+y_mult*y)*num0                                    +(y0+y_mult*y)*num0
360                                    +(x0+x_mult*x);                                    +(x0+x_mult*x);
361                  if (!isnan(values[srcIndex])) {                  if (!isnan(values[srcIndex])) {
362                      for (index_t m2=0; m2<params.multiplier[2]; m2++) {                      for (index_t m2=0; m2<params.multiplier[2]; m2++) {
363                          for (index_t m1=0; m1<params.multiplier[1]; m1++) {                          for (index_t m1=0; m1<params.multiplier[1]; m1++) {
364                              for (index_t m0=0; m0<params.multiplier[0]; m0++) {                              for (index_t m0=0; m0<params.multiplier[0]; m0++) {
365                                  const int dataIndex = baseIndex+m0                                  const dim_t dataIndex = baseIndex+m0
366                                                 +m1*myN0                                                 +m1*myN0
367                                                 +m2*myN0*myN1;                                                 +m2*myN0*myN1;
368                                  double* dest = out.getSampleDataRW(dataIndex);                                  double* dest = out.getSampleDataRW(dataIndex);
# Line 428  void Brick::readBinaryGridImpl(escript:: Line 428  void Brick::readBinaryGridImpl(escript::
428                                 const ReaderParameters& params) const                                 const ReaderParameters& params) const
429  {  {
430      // check destination function space      // check destination function space
431      int myN0, myN1, myN2;      dim_t myN0, myN1, myN2;
432      if (out.getFunctionSpace().getTypeCode() == Nodes) {      if (out.getFunctionSpace().getTypeCode() == Nodes) {
433          myN0 = m_NN[0];          myN0 = m_NN[0];
434          myN1 = m_NN[1];          myN1 = m_NN[1];
# Line 462  void Brick::readBinaryGridImpl(escript:: Line 462  void Brick::readBinaryGridImpl(escript::
462      }      }
463      f.seekg(0, ios::end);      f.seekg(0, ios::end);
464      const int numComp = out.getDataPointSize();      const int numComp = out.getDataPointSize();
465      const int filesize = f.tellg();      const dim_t filesize = f.tellg();
466      const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);      const dim_t reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
467      if (filesize < reqsize) {      if (filesize < reqsize) {
468          f.close();          f.close();
469          throw RipleyException("readBinaryGrid(): not enough data in file");          throw RipleyException("readBinaryGrid(): not enough data in file");
# Line 483  void Brick::readBinaryGridImpl(escript:: Line 483  void Brick::readBinaryGridImpl(escript::
483      // now determine how much this rank has to write      // now determine how much this rank has to write
484    
485      // first coordinates in data object to write to      // first coordinates in data object to write to
486      const int first0 = max(0, params.first[0]-m_offset[0]);      const dim_t first0 = max(0, params.first[0]-m_offset[0]);
487      const int first1 = max(0, params.first[1]-m_offset[1]);      const dim_t first1 = max(0, params.first[1]-m_offset[1]);
488      const int first2 = max(0, params.first[2]-m_offset[2]);      const dim_t first2 = max(0, params.first[2]-m_offset[2]);
489      // indices to first value in file (not accounting for reverse yet)      // indices to first value in file (not accounting for reverse yet)
490      int idx0 = max(0, m_offset[0]-params.first[0]);      dim_t idx0 = max(0, (m_offset[0]/params.multiplier[0])-params.first[0]);
491      int idx1 = max(0, m_offset[1]-params.first[1]);      dim_t idx1 = max(0, (m_offset[1]/params.multiplier[1])-params.first[1]);
492      int idx2 = max(0, m_offset[2]-params.first[2]);      dim_t idx2 = max(0, (m_offset[2]/params.multiplier[2])-params.first[2]);
493        // if restX > 0 the first value in the respective dimension has been
494        // written restX times already in a previous rank so this rank only
495        // contributes (multiplier-rank) copies of that value
496        const dim_t rest0 = m_offset[0]%params.multiplier[0];
497        const dim_t rest1 = m_offset[1]%params.multiplier[1];
498        const dim_t rest2 = m_offset[2]%params.multiplier[2];
499    
500      // number of values to read      // number of values to read
501      const int num0 = min(params.numValues[0]-idx0, myN0-first0);      const dim_t num0 = min(params.numValues[0]-idx0, myN0-first0);
502      const int num1 = min(params.numValues[1]-idx1, myN1-first1);      const dim_t num1 = min(params.numValues[1]-idx1, myN1-first1);
503      const int num2 = min(params.numValues[2]-idx2, myN2-first2);      const dim_t num2 = min(params.numValues[2]-idx2, myN2-first2);
504    
505      // make sure we read the right block if going backwards through file      // make sure we read the right block if going backwards through file
506      if (params.reverse[2])      if (params.reverse[2])
# Line 506  void Brick::readBinaryGridImpl(escript:: Line 513  void Brick::readBinaryGridImpl(escript::
513      vector<ValueType> values(num0*numComp);      vector<ValueType> values(num0*numComp);
514      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
515    
516      for (int z=0; z<num2; z++) {      for (dim_t z=0; z<num2; z++) {
517          for (int y=0; y<num1; y++) {          const dim_t m2limit = (z==0 ? params.multiplier[2]-rest2 : params.multiplier[2]);
518              const int fileofs = numComp*(idx0 +          dim_t dataZbase = first2 + z*params.multiplier[2];
519            if (z>0)
520                dataZbase -= rest2;
521    
522            for (dim_t y=0; y<num1; y++) {
523                const dim_t fileofs = numComp*(idx0 +
524                                  (idx1+y)*params.numValues[0] +                                  (idx1+y)*params.numValues[0] +
525                                  (idx2+z_mult*z)*params.numValues[0]*params.numValues[1]);                                  (idx2+z_mult*z)*params.numValues[0]*params.numValues[1]);
526              f.seekg(fileofs*sizeof(ValueType));              f.seekg(fileofs*sizeof(ValueType));
527              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
528                const dim_t m1limit = (y==0 ? params.multiplier[1]-rest1 : params.multiplier[1]);
529              for (int x=0; x<num0; x++) {              dim_t dataYbase = first1 + y*params.multiplier[1];
530                  const int baseIndex = first0+x*params.multiplier[0]              if (y>0)
531                                       +(first1+y*params.multiplier[1])*myN0                  dataYbase -= rest1;
532                                       +(first2+z*params.multiplier[2])*myN0*myN1;  
533                  for (int m2=0; m2<params.multiplier[2]; m2++) {              for (dim_t x=0; x<num0; x++) {
534                      for (int m1=0; m1<params.multiplier[1]; m1++) {                  const dim_t m0limit = (x==0 ? params.multiplier[0]-rest0 : params.multiplier[0]);
535                          for (int m0=0; m0<params.multiplier[0]; m0++) {                  dim_t dataXbase = first0 + x*params.multiplier[0];
536                              const int dataIndex = baseIndex+m0                  if (x>0)
537                                             +m1*myN0                      dataXbase -= rest0;
538                                             +m2*myN0*myN1;                  // write a block of mult0 x mult1 x mult2 identical values into
539                    // Data object
540                    for (dim_t m2=0; m2 < m2limit; m2++) {
541                        const dim_t dataZ = dataZbase + m2;
542                        if (dataZ >= myN2)
543                            break;
544                        for (dim_t m1=0; m1 < m1limit; m1++) {
545                            const dim_t dataY = dataYbase + m1;
546                            if (dataY >= myN1)
547                                break;
548                            for (dim_t m0=0; m0 < m0limit; m0++) {
549                                const dim_t dataX = dataXbase + m0;
550                                if (dataX >= myN0)
551                                    break;
552                                const dim_t dataIndex = dataX + dataY*myN0 + dataZ*myN0*myN1;
553                              double* dest = out.getSampleDataRW(dataIndex);                              double* dest = out.getSampleDataRW(dataIndex);
554                              for (int c=0; c<numComp; c++) {                              for (int c=0; c<numComp; c++) {
555                                  ValueType val = values[x*numComp+c];                                  ValueType val = values[x*numComp+c];
# Line 559  void Brick::readBinaryGridZippedImpl(esc Line 585  void Brick::readBinaryGridZippedImpl(esc
585                                 const ReaderParameters& params) const                                 const ReaderParameters& params) const
586  {  {
587      // check destination function space      // check destination function space
588      int myN0, myN1, myN2;      dim_t myN0, myN1, myN2;
589      if (out.getFunctionSpace().getTypeCode() == Nodes) {      if (out.getFunctionSpace().getTypeCode() == Nodes) {
590          myN0 = m_NN[0];          myN0 = m_NN[0];
591          myN1 = m_NN[1];          myN1 = m_NN[1];
# Line 591  void Brick::readBinaryGridZippedImpl(esc Line 617  void Brick::readBinaryGridZippedImpl(esc
617      }      }
618      f.seekg(0, ios::end);      f.seekg(0, ios::end);
619      const int numComp = out.getDataPointSize();      const int numComp = out.getDataPointSize();
620      int filesize = f.tellg();      dim_t filesize = f.tellg();
621      f.seekg(0, ios::beg);      f.seekg(0, ios::beg);
622      std::vector<char> compressed(filesize);      vector<char> compressed(filesize);
623      f.read((char*)&compressed[0], filesize);      f.read((char*)&compressed[0], filesize);
624      f.close();      f.close();
625      std::vector<char> decompressed = unzip(compressed);      vector<char> decompressed = unzip(compressed);
626      filesize = decompressed.size();      filesize = decompressed.size();
627      const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);      const dim_t reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
628      if (filesize < reqsize) {      if (filesize < reqsize) {
629          throw RipleyException("readBinaryGridFromZipped(): not enough data in file");          throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
630      }      }
# Line 616  void Brick::readBinaryGridZippedImpl(esc Line 642  void Brick::readBinaryGridZippedImpl(esc
642      // now determine how much this rank has to write      // now determine how much this rank has to write
643    
644      // first coordinates in data object to write to      // first coordinates in data object to write to
645      const int first0 = max(0, params.first[0]-m_offset[0]);      const dim_t first0 = max(0, params.first[0]-m_offset[0]);
646      const int first1 = max(0, params.first[1]-m_offset[1]);      const dim_t first1 = max(0, params.first[1]-m_offset[1]);
647      const int first2 = max(0, params.first[2]-m_offset[2]);      const dim_t first2 = max(0, params.first[2]-m_offset[2]);
648      // indices to first value in file      // indices to first value in file
649      const int idx0 = max(0, m_offset[0]-params.first[0]);      const dim_t idx0 = max(0, m_offset[0]-params.first[0]);
650      const int idx1 = max(0, m_offset[1]-params.first[1]);      const dim_t idx1 = max(0, m_offset[1]-params.first[1]);
651      const int idx2 = max(0, m_offset[2]-params.first[2]);      const dim_t idx2 = max(0, m_offset[2]-params.first[2]);
652      // number of values to read      // number of values to read
653      const int num0 = min(params.numValues[0]-idx0, myN0-first0);      const dim_t num0 = min(params.numValues[0]-idx0, myN0-first0);
654      const int num1 = min(params.numValues[1]-idx1, myN1-first1);      const dim_t num1 = min(params.numValues[1]-idx1, myN1-first1);
655      const int num2 = min(params.numValues[2]-idx2, myN2-first2);      const dim_t num2 = min(params.numValues[2]-idx2, myN2-first2);
656    
657      out.requireWrite();      out.requireWrite();
658      vector<ValueType> values(num0*numComp);      vector<ValueType> values(num0*numComp);
659      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
660    
661      for (int z=0; z<num2; z++) {      for (dim_t z=0; z<num2; z++) {
662          for (int y=0; y<num1; y++) {          for (dim_t y=0; y<num1; y++) {
663              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]              const dim_t 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 (dim_t x=0; x<num0; x++) {
668                  const int baseIndex = first0+x*params.multiplier[0]                  const dim_t baseIndex = first0+x*params.multiplier[0]
669                                       +(first1+y*params.multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
670                                       +(first2+z*params.multiplier[2])*myN0*myN1;                                       +(first2+z*params.multiplier[2])*myN0*myN1;
671                  for (int m2=0; m2<params.multiplier[2]; m2++) {                  for (dim_t m2=0; m2<params.multiplier[2]; m2++) {
672                      for (int m1=0; m1<params.multiplier[1]; m1++) {                      for (dim_t m1=0; m1<params.multiplier[1]; m1++) {
673                          for (int m0=0; m0<params.multiplier[0]; m0++) {                          for (dim_t m0=0; m0<params.multiplier[0]; m0++) {
674                              const int dataIndex = baseIndex+m0                              const dim_t dataIndex = baseIndex+m0
675                                             +m1*myN0                                             +m1*myN0
676                                             +m2*myN0*myN1;                                             +m2*myN0*myN1;
677                              double* dest = out.getSampleDataRW(dataIndex);                              double* dest = out.getSampleDataRW(dataIndex);
# Line 697  void Brick::writeBinaryGridImpl(const es Line 723  void Brick::writeBinaryGridImpl(const es
723                                  const string& filename, int byteOrder) const                                  const string& filename, int byteOrder) const
724  {  {
725      // check function space and determine number of points      // check function space and determine number of points
726      int myN0, myN1, myN2;      dim_t myN0, myN1, myN2;
727      int totalN0, totalN1, totalN2;      dim_t totalN0, totalN1, totalN2;
728      if (in.getFunctionSpace().getTypeCode() == Nodes) {      if (in.getFunctionSpace().getTypeCode() == Nodes) {
729          myN0 = m_NN[0];          myN0 = m_NN[0];
730          myN1 = m_NN[1];          myN1 = m_NN[1];
# Line 719  void Brick::writeBinaryGridImpl(const es Line 745  void Brick::writeBinaryGridImpl(const es
745    
746      const int numComp = in.getDataPointSize();      const int numComp = in.getDataPointSize();
747      const int dpp = in.getNumDataPointsPerSample();      const int dpp = in.getNumDataPointsPerSample();
748      const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1*totalN2;      const dim_t fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1*totalN2;
749    
750      if (numComp > 1 || dpp > 1)      if (numComp > 1 || dpp > 1)
751          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
# Line 731  void Brick::writeBinaryGridImpl(const es Line 757  void Brick::writeBinaryGridImpl(const es
757    
758      for (index_t z=0; z<myN2; z++) {      for (index_t z=0; z<myN2; z++) {
759          for (index_t y=0; y<myN1; y++) {          for (index_t y=0; y<myN1; y++) {
760              const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0              const dim_t fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0
761                                  +(m_offset[2]+z)*totalN0*totalN1)*sizeof(ValueType);                                  +(m_offset[2]+z)*totalN0*totalN1)*sizeof(ValueType);
762              ostringstream oss;              ostringstream oss;
763    
# Line 756  void Brick::writeBinaryGridImpl(const es Line 782  void Brick::writeBinaryGridImpl(const es
782      fw.close();      fw.close();
783  }  }
784    
785    void Brick::write(const std::string& filename) const
786    {
787        throw RipleyException("write: not supported");
788    }
789    
790  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
791  {  {
792  #if USE_SILO  #if USE_SILO
# Line 764  void Brick::dump(const string& fileName) Line 795  void Brick::dump(const string& fileName)
795          fn+=".silo";          fn+=".silo";
796      }      }
797    
798      int driver=DB_HDF5;          int driver=DB_HDF5;
799      string siloPath;      string siloPath;
800      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
801    
# Line 820  void Brick::dump(const string& fileName) Line 851  void Brick::dump(const string& fileName)
851      boost::scoped_ptr<double> y(new double[m_NN[1]]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
852      boost::scoped_ptr<double> z(new double[m_NN[2]]);      boost::scoped_ptr<double> z(new double[m_NN[2]]);
853      double* coords[3] = { x.get(), y.get(), z.get() };      double* coords[3] = { x.get(), y.get(), z.get() };
854      const int NN0 = m_NN[0];      const dim_t NN0 = m_NN[0];
855      const int NN1 = m_NN[1];      const dim_t NN1 = m_NN[1];
856      const int NN2 = m_NN[2];      const dim_t NN2 = m_NN[2];
857    
858  #pragma omp parallel  #pragma omp parallel
859      {      {
# Line 839  void Brick::dump(const string& fileName) Line 870  void Brick::dump(const string& fileName)
870              coords[2][i2]=getLocalCoordinate(i2, 2);              coords[2][i2]=getLocalCoordinate(i2, 2);
871          }          }
872      }      }
873      int* dims = const_cast<int*>(getNumNodesPerDim());      dim_t* dims = const_cast<dim_t*>(getNumNodesPerDim());
874    
875      // write mesh      // write mesh
876      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,
# Line 850  void Brick::dump(const string& fileName) Line 881  void Brick::dump(const string& fileName)
881              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
882    
883      // write element ids      // write element ids
884      dims = const_cast<int*>(getNumElementsPerDim());      dims = const_cast<dim_t*>(getNumElementsPerDim());
885      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
886              dims, 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
887    
# Line 905  void Brick::dump(const string& fileName) Line 936  void Brick::dump(const string& fileName)
936  #endif  #endif
937  }  }
938    
939  const int* Brick::borrowSampleReferenceIDs(int fsType) const  const dim_t* Brick::borrowSampleReferenceIDs(int fsType) const
940  {  {
941      switch (fsType) {      switch (fsType) {
942          case Nodes:          case Nodes:
# Line 986  bool Brick::ownSample(int fsType, index_ Line 1017  bool Brick::ownSample(int fsType, index_
1017    
1018  void Brick::setToNormal(escript::Data& out) const  void Brick::setToNormal(escript::Data& out) const
1019  {  {
1020      const int NE0 = m_NE[0];      const dim_t NE0 = m_NE[0];
1021      const int NE1 = m_NE[1];      const dim_t NE1 = m_NE[1];
1022      const int NE2 = m_NE[2];      const dim_t NE2 = m_NE[2];
1023    
1024      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1025          out.requireWrite();          out.requireWrite();
# Line 1170  void Brick::setToSize(escript::Data& out Line 1201  void Brick::setToSize(escript::Data& out
1201          out.requireWrite();          out.requireWrite();
1202          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
1203          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);
1204          const int NE = getNumElements();          const dim_t NE = getNumElements();
1205  #pragma omp parallel for  #pragma omp parallel for
1206          for (index_t k = 0; k < NE; ++k) {          for (index_t k = 0; k < NE; ++k) {
1207              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 1180  void Brick::setToSize(escript::Data& out Line 1211  void Brick::setToSize(escript::Data& out
1211              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1212          out.requireWrite();          out.requireWrite();
1213          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
1214          const int NE0 = m_NE[0];          const dim_t NE0 = m_NE[0];
1215          const int NE1 = m_NE[1];          const dim_t NE1 = m_NE[1];
1216          const int NE2 = m_NE[2];          const dim_t NE2 = m_NE[2];
1217  #pragma omp parallel  #pragma omp parallel
1218          {          {
1219              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
# Line 1287  void Brick::assembleCoordinates(escript: Line 1318  void Brick::assembleCoordinates(escript:
1318      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
1319          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
1320    
1321      const int NN0 = m_NN[0];      const dim_t NN0 = m_NN[0];
1322      const int NN1 = m_NN[1];      const dim_t NN1 = m_NN[1];
1323      const int NN2 = m_NN[2];      const dim_t NN2 = m_NN[2];
1324      arg.requireWrite();      arg.requireWrite();
1325  #pragma omp parallel for  #pragma omp parallel for
1326      for (dim_t i2 = 0; i2 < NN2; i2++) {      for (dim_t i2 = 0; i2 < NN2; i2++) {
# Line 1315  void Brick::assembleGradient(escript::Da Line 1346  void Brick::assembleGradient(escript::Da
1346      const double C4 = .5;      const double C4 = .5;
1347      const double C5 = .62200846792814621559;      const double C5 = .62200846792814621559;
1348      const double C6 = .78867513459481288225;      const double C6 = .78867513459481288225;
1349      const int NE0 = m_NE[0];      const dim_t NE0 = m_NE[0];
1350      const int NE1 = m_NE[1];      const dim_t NE1 = m_NE[1];
1351      const int NE2 = m_NE[2];      const dim_t NE2 = m_NE[2];
1352    
1353      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
1354          out.requireWrite();          out.requireWrite();
# Line 2096  void Brick::dofToNodes(escript::Data& ou Line 2127  void Brick::dofToNodes(escript::Data& ou
2127      }      }
2128  }  }
2129    
2130    //protected
2131    paso::SystemMatrixPattern_ptr Brick::getPasoMatrixPattern(
2132                                                        bool reducedRowOrder,
2133                                                        bool reducedColOrder) const
2134    {
2135        if (m_pattern.get())
2136            return m_pattern;
2137    
2138        // first call to this method -> create the pattern, then return it
2139        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2140        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2141        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2142        const dim_t numDOF = nDOF0*nDOF1*nDOF2;
2143        const dim_t numShared = m_connector->send->numSharedComponents;
2144        const index_t* sendShared = m_connector->send->shared;
2145        const int x = m_mpiInfo->rank%m_NX[0];
2146        const int y = m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2147        const int z = m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
2148    
2149        // these are for the couple blocks
2150        vector<IndexVector> colIndices(numDOF);
2151        vector<IndexVector> rowIndices(numShared);
2152    
2153        for (dim_t i=0; i < m_connector->send->numNeighbors; i++) {
2154            const dim_t start = m_connector->send->offsetInShared[i];
2155            const dim_t end = m_connector->send->offsetInShared[i+1];
2156            // location of neighbour rank relative to this rank
2157            const int xDiff = m_connector->send->neighbor[i]%m_NX[0] - x;
2158            const int yDiff = m_connector->send->neighbor[i]%(m_NX[0]*m_NX[1])/m_NX[0] - y;
2159            const int zDiff = m_connector->send->neighbor[i]/(m_NX[0]*m_NX[1]) - z;
2160            
2161            if (xDiff==0 && yDiff==0) {
2162                // sharing front or back plane
2163                for (dim_t j = start; j < end; j++) {
2164                    const dim_t i0 = (j-start)%nDOF0;
2165                    const dim_t i1 = (j-start)/nDOF0;
2166                    if (i0 > 0) {
2167                        if (i1 > 0)
2168                            doublyLink(colIndices, rowIndices, sendShared[j-1-nDOF0], j);
2169                        doublyLink(colIndices, rowIndices, sendShared[j-1], j);
2170                        if (i1 < nDOF1-1)
2171                            doublyLink(colIndices, rowIndices, sendShared[j-1+nDOF0], j);
2172                    }
2173                    if (i1 > 0)
2174                        doublyLink(colIndices, rowIndices, sendShared[j-nDOF0], j);
2175                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2176                    if (i1 < nDOF1-1)
2177                        doublyLink(colIndices, rowIndices, sendShared[j+nDOF0], j);
2178                    if (i0 < nDOF0-1) {
2179                        if (i1 > 0)
2180                            doublyLink(colIndices, rowIndices, sendShared[j+1-nDOF0], j);
2181                        doublyLink(colIndices, rowIndices, sendShared[j+1], j);
2182                        if (i1 < nDOF1-1)
2183                            doublyLink(colIndices, rowIndices, sendShared[j+1+nDOF0], j);
2184                    }
2185                }
2186            } else if (xDiff==0 && zDiff==0) {
2187                // sharing top or bottom plane
2188                for (dim_t j = start; j < end; j++) {
2189                    const dim_t i0 = (j-start)%nDOF0;
2190                    const dim_t i1 = (j-start)/nDOF0;
2191                    if (i0 > 0) {
2192                        if (i1 > 0)
2193                            doublyLink(colIndices, rowIndices, sendShared[j]-1-nDOF0*nDOF1, j);
2194                        doublyLink(colIndices, rowIndices, sendShared[j]-1, j);
2195                        if (i1 < nDOF2-1)
2196                            doublyLink(colIndices, rowIndices, sendShared[j]-1+nDOF0*nDOF1, j);
2197                    }
2198                    if (i1 > 0)
2199                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2200                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2201                    if (i1 < nDOF2-1)
2202                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2203                    if (i0 < nDOF0-1) {
2204                        if (i1 > 0)
2205                            doublyLink(colIndices, rowIndices, sendShared[j]+1-nDOF0*nDOF1, j);
2206                        doublyLink(colIndices, rowIndices, sendShared[j]+1, j);
2207                        if (i1 < nDOF2-1)
2208                            doublyLink(colIndices, rowIndices, sendShared[j]+1+nDOF0*nDOF1, j);
2209                    }
2210                }
2211            } else if (yDiff==0 && zDiff==0) {
2212                // sharing left or right plane
2213                for (dim_t j = start; j < end; j++) {
2214                    const dim_t i0 = (j-start)%nDOF1;
2215                    const dim_t i1 = (j-start)/nDOF1;
2216                    if (i0 > 0) {
2217                        if (i1 > 0)
2218                            doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0-nDOF0*nDOF1, j);
2219                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0, j);
2220                        if (i1 < nDOF2-1)
2221                            doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0+nDOF0*nDOF1, j);
2222                    }
2223                    if (i1 > 0)
2224                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2225                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2226                    if (i1 < nDOF2-1)
2227                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2228                    if (i0 < nDOF1-1) {
2229                        if (i1 > 0)
2230                            doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0-nDOF0*nDOF1, j);
2231                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0, j);
2232                        if (i1 < nDOF2-1)
2233                            doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0+nDOF0*nDOF1, j);
2234                    }
2235                }
2236            } else if (xDiff == 0) {
2237                // sharing an edge in x direction
2238                for (dim_t j = start; j < end; j++) {
2239                    if (j > start)
2240                        doublyLink(colIndices, rowIndices, sendShared[j]-1, j);
2241                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2242                    if (j < end-1)
2243                        doublyLink(colIndices, rowIndices, sendShared[j]+1, j);
2244                }
2245            } else if (yDiff == 0) {
2246                // sharing an edge in y direction
2247                for (dim_t j = start; j < end; j++) {
2248                    if (j > start)
2249                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0, j);
2250                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2251                    if (j < end-1)
2252                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0, j);
2253                }
2254            } else if (zDiff == 0) {
2255                // sharing an edge in z direction
2256                for (dim_t j = start; j < end; j++) {
2257                    if (j > start)
2258                        doublyLink(colIndices, rowIndices, sendShared[j]-nDOF0*nDOF1, j);
2259                    doublyLink(colIndices, rowIndices, sendShared[j], j);
2260                    if (j < end-1)
2261                        doublyLink(colIndices, rowIndices, sendShared[j]+nDOF0*nDOF1, j);
2262                }
2263            } else {
2264                // sharing a node
2265                doublyLink(colIndices, rowIndices, sendShared[start], start);
2266            }
2267        }
2268    
2269    #pragma omp parallel for
2270        for (dim_t i = 0; i < numShared; i++) {
2271            sort(rowIndices[i].begin(), rowIndices[i].end());
2272        }
2273    
2274        // create main and couple blocks
2275        paso::Pattern_ptr mainPattern = createPasoPattern(getConnections(), numDOF);
2276        paso::Pattern_ptr colPattern = createPasoPattern(colIndices, numShared);
2277        paso::Pattern_ptr rowPattern = createPasoPattern(rowIndices, numDOF);
2278    
2279        // allocate paso distribution
2280        paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2281                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2282    
2283        // finally create the system matrix pattern
2284        m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2285                distribution, distribution, mainPattern, colPattern, rowPattern,
2286                m_connector, m_connector));
2287    
2288        // useful debug output
2289        /*
2290        cout << "--- colIndices ---" << endl;
2291        for (size_t i=0; i<colIndices.size(); i++) {
2292            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
2293        }
2294        cout << "--- rowIndices ---" << endl;
2295        for (size_t i=0; i<rowIndices.size(); i++) {
2296            cout << "rowIndices[" << i << "].size()=" << rowIndices[i].size() << endl;
2297        }
2298        */
2299        /*
2300        cout << "--- main_pattern ---" << endl;
2301        cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
2302        for (size_t i=0; i<mainPattern->numOutput+1; i++) {
2303            cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
2304        }
2305        for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
2306            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
2307        }
2308        */
2309        /*
2310        cout << "--- colCouple_pattern ---" << endl;
2311        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
2312        for (size_t i=0; i<colPattern->numOutput+1; i++) {
2313            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
2314        }
2315        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
2316            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
2317        }
2318        */
2319        /*
2320        cout << "--- rowCouple_pattern ---" << endl;
2321        cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
2322        for (size_t i=0; i<rowPattern->numOutput+1; i++) {
2323            cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
2324        }
2325        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
2326            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2327        }
2328        */
2329    
2330        return m_pattern;
2331    }
2332    
2333  //private  //private
2334  void Brick::populateSampleIds()  void Brick::populateSampleIds()
2335  {  {
# Line 2114  void Brick::populateSampleIds() Line 2348  void Brick::populateSampleIds()
2348          m_nodeDistribution[k]=k*numDOF;          m_nodeDistribution[k]=k*numDOF;
2349      }      }
2350      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2351        
2352      try {      try {
2353          m_nodeId.resize(getNumNodes());          m_nodeId.resize(getNumNodes());
2354          m_dofId.resize(numDOF);          m_dofId.resize(numDOF);
2355          m_elementId.resize(getNumElements());          m_elementId.resize(getNumElements());
2356      } catch (const std::length_error& le) {      } catch (const length_error& le) {
2357          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.");
2358      }      }
2359        
2360      // populate face element counts      // populate face element counts
2361      //left      //left
2362      if (m_offset[0]==0)      if (m_offset[0]==0)
# Line 2341  void Brick::populateSampleIds() Line 2575  void Brick::populateSampleIds()
2575      setTagMap("front", FRONT);      setTagMap("front", FRONT);
2576      setTagMap("back", BACK);      setTagMap("back", BACK);
2577      updateTagsInUse(FaceElements);      updateTagsInUse(FaceElements);
2578    
2579        populateDofMap();
2580    }
2581    
2582    //private
2583    vector<IndexVector> Brick::getConnections() const
2584    {
2585        // returns a vector v of size numDOF where v[i] is a vector with indices
2586        // of DOFs connected to i (up to 27 in 3D)
2587        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2588        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2589        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2590        const dim_t M = nDOF0*nDOF1*nDOF2;
2591        vector<IndexVector> indices(M);
2592    
2593    #pragma omp parallel for
2594        for (index_t i=0; i < M; i++) {
2595            const index_t x = i % nDOF0;
2596            const index_t y = i % (nDOF0*nDOF1)/nDOF0;
2597            const index_t z = i / (nDOF0*nDOF1);
2598            // loop through potential neighbours and add to index if positions are
2599            // within bounds
2600            for (int i2=z-1; i2<z+2; i2++) {
2601                for (int i1=y-1; i1<y+2; i1++) {
2602                    for (int i0=x-1; i0<x+2; i0++) {
2603                        if (i0>=0 && i1>=0 && i2>=0
2604                                && i0<nDOF0 && i1<nDOF1 && i2<nDOF2) {
2605                            indices[i].push_back(i2*nDOF0*nDOF1 + i1*nDOF0 + i0);
2606                        }
2607                    }
2608                }
2609            }
2610        }
2611        return indices;
2612  }  }
2613    
2614  //private  //private
2615  void Brick::createPattern()  void Brick::populateDofMap()
2616  {  {
2617      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2618      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
# Line 2365  void Brick::createPattern() Line 2633  void Brick::createPattern()
2633          }          }
2634      }      }
2635    
     // build list of shared components and neighbours by looping through  
     // all potential neighbouring ranks and checking if positions are  
     // within bounds  
2636      const dim_t numDOF=nDOF0*nDOF1*nDOF2;      const dim_t numDOF=nDOF0*nDOF1*nDOF2;
     vector<IndexVector> colIndices(numDOF); // for the couple blocks  
2637      RankVector neighbour;      RankVector neighbour;
2638      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
2639      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
2640      int numShared=0, expectedShared=0;;      dim_t numShared=0;
2641      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
2642      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];
2643      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
2644      for (int i2=-1; i2<2; i2++) {  
2645          for (int i1=-1; i1<2; i1++) {      // build list of shared components and neighbours by looping through
2646              for (int i0=-1; i0<2; i0++) {      // all potential neighbouring ranks and checking if positions are
2647                  // skip this rank      // within bounds
                 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> rowIndices(expectedShared);  
       
2648      for (int i2=-1; i2<2; i2++) {      for (int i2=-1; i2<2; i2++) {
2649          for (int i1=-1; i1<2; i1++) {          for (int i1=-1; i1<2; i1++) {
2650              for (int i0=-1; i0<2; i0++) {              for (int i0=-1; i0<2; i0++) {
# Line 2426  void Brick::createPattern() Line 2661  void Brick::createPattern()
2661                          // sharing front or back plane                          // sharing front or back plane
2662                          offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF1);                          offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF1);
2663                          for (dim_t i=0; i<nDOF1; i++) {                          for (dim_t i=0; i<nDOF1; i++) {
2664                              const int firstDOF=(i2==-1 ? i*nDOF0                              const dim_t firstDOF=(i2==-1 ? i*nDOF0
2665                                      : i*nDOF0 + nDOF0*nDOF1*(nDOF2-1));                                      : i*nDOF0 + nDOF0*nDOF1*(nDOF2-1));
2666                              const int firstNode=(i2==-1 ? left+(i+bottom)*m_NN[0]                              const dim_t firstNode=(i2==-1 ? left+(i+bottom)*m_NN[0]
2667                                      : left+(i+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1));                                      : left+(i+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1));
2668                              for (dim_t j=0; j<nDOF0; j++, numShared++) {                              for (dim_t j=0; j<nDOF0; j++, numShared++) {
2669                                  sendShared.push_back(firstDOF+j);                                  sendShared.push_back(firstDOF+j);
2670                                  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);  
                                 }  
2671                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2672                              }                              }
2673                          }                          }
# Line 2459  void Brick::createPattern() Line 2675  void Brick::createPattern()
2675                          // sharing top or bottom plane                          // sharing top or bottom plane
2676                          offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF2);
2677                          for (dim_t i=0; i<nDOF2; i++) {                          for (dim_t i=0; i<nDOF2; i++) {
2678                              const int firstDOF=(i1==-1 ? i*nDOF0*nDOF1                              const dim_t firstDOF=(i1==-1 ? i*nDOF0*nDOF1
2679                                      : nDOF0*((i+1)*nDOF1-1));                                      : nDOF0*((i+1)*nDOF1-1));
2680                              const int firstNode=(i1==-1 ?                              const dim_t firstNode=(i1==-1 ?
2681                                      left+(i+front)*m_NN[0]*m_NN[1]                                      left+(i+front)*m_NN[0]*m_NN[1]
2682                                      : left+m_NN[0]*((i+1+front)*m_NN[1]-1));                                      : left+m_NN[0]*((i+1+front)*m_NN[1]-1));
2683                              for (dim_t j=0; j<nDOF0; j++, numShared++) {                              for (dim_t j=0; j<nDOF0; j++, numShared++) {
2684                                  sendShared.push_back(firstDOF+j);                                  sendShared.push_back(firstDOF+j);
2685                                  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);  
                                 }  
2686                                  m_dofMap[firstNode+j]=numDOF+numShared;                                  m_dofMap[firstNode+j]=numDOF+numShared;
2687                              }                              }
2688                          }                          }
# Line 2493  void Brick::createPattern() Line 2690  void Brick::createPattern()
2690                          // sharing left or right plane                          // sharing left or right plane
2691                          offsetInShared.push_back(offsetInShared.back()+nDOF1*nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF1*nDOF2);
2692                          for (dim_t i=0; i<nDOF2; i++) {                          for (dim_t i=0; i<nDOF2; i++) {
2693                              const int firstDOF=(i0==-1 ? i*nDOF0*nDOF1                              const dim_t firstDOF=(i0==-1 ? i*nDOF0*nDOF1
2694                                      : nDOF0*(1+i*nDOF1)-1);                                      : nDOF0*(1+i*nDOF1)-1);
2695                              const int firstNode=(i0==-1 ?                              const dim_t firstNode=(i0==-1 ?
2696                                      (bottom+(i+front)*m_NN[1])*m_NN[0]                                      (bottom+(i+front)*m_NN[1])*m_NN[0]
2697                                      : (bottom+1+(i+front)*m_NN[1])*m_NN[0]-1);                                      : (bottom+1+(i+front)*m_NN[1])*m_NN[0]-1);
2698                              for (dim_t j=0; j<nDOF1; j++, numShared++) {                              for (dim_t j=0; j<nDOF1; j++, numShared++) {
2699                                  sendShared.push_back(firstDOF+j*nDOF0);                                  sendShared.push_back(firstDOF+j*nDOF0);
2700                                  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);  
                                 }  
2701                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;                                  m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2702                              }                              }
2703                          }                          }
2704                      } else if (i0==0) {                      } else if (i0==0) {
2705                          // sharing an edge in x direction                          // sharing an edge in x direction
2706                          offsetInShared.push_back(offsetInShared.back()+nDOF0);                          offsetInShared.push_back(offsetInShared.back()+nDOF0);
2707                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)                          const dim_t firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
2708                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2709                          const int firstNode=left+(i1+1)/2*m_NN[0]*(m_NN[1]-1)                          const dim_t firstNode=left+(i1+1)/2*m_NN[0]*(m_NN[1]-1)
2710                                              +(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);
2711                          for (dim_t i=0; i<nDOF0; i++, numShared++) {                          for (dim_t i=0; i<nDOF0; i++, numShared++) {
2712                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
2713                              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);  
2714                              m_dofMap[firstNode+i]=numDOF+numShared;                              m_dofMap[firstNode+i]=numDOF+numShared;
2715                          }                          }
2716                      } else if (i1==0) {                      } else if (i1==0) {
2717                          // sharing an edge in y direction                          // sharing an edge in y direction
2718                          offsetInShared.push_back(offsetInShared.back()+nDOF1);                          offsetInShared.push_back(offsetInShared.back()+nDOF1);
2719                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const dim_t firstDOF=(i0+1)/2*(nDOF0-1)
2720                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2721                          const int firstNode=bottom*m_NN[0]                          const dim_t firstNode=bottom*m_NN[0]
2722                                              +(i0+1)/2*(m_NN[0]-1)                                              +(i0+1)/2*(m_NN[0]-1)
2723                                              +(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);
2724                          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          for (dim_t i=0; i<nDOF1; i++, numShared++) {
2725                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
2726                              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);  
2727                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2728                          }                          }
2729                      } else if (i2==0) {                      } else if (i2==0) {
2730                          // sharing an edge in z direction                          // sharing an edge in z direction
2731                          offsetInShared.push_back(offsetInShared.back()+nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF2);
2732                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const dim_t firstDOF=(i0+1)/2*(nDOF0-1)
2733                                             +(i1+1)/2*nDOF0*(nDOF1-1);                                             +(i1+1)/2*nDOF0*(nDOF1-1);
2734                          const int firstNode=front*m_NN[0]*m_NN[1]                          const dim_t firstNode=front*m_NN[0]*m_NN[1]
2735                                              +(i0+1)/2*(m_NN[0]-1)                                              +(i0+1)/2*(m_NN[0]-1)
2736                                              +(i1+1)/2*m_NN[0]*(m_NN[1]-1);                                              +(i1+1)/2*m_NN[0]*(m_NN[1]-1);
2737                          for (dim_t i=0; i<nDOF2; i++, numShared++) {                          for (dim_t i=0; i<nDOF2; i++, numShared++) {
2738                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2739                              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);  
2740                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;                              m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2741                          }                          }
2742                      } else {                      } else {
2743                          // sharing a node                          // sharing a node
2744                          const int dof=(i0+1)/2*(nDOF0-1)                          const dim_t dof = (i0+1)/2*(nDOF0-1)
2745                                        +(i1+1)/2*nDOF0*(nDOF1-1)                                         +(i1+1)/2*nDOF0*(nDOF1-1)
2746                                        +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                         +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2747                          const int node=(i0+1)/2*(m_NN[0]-1)                          const dim_t node = (i0+1)/2*(m_NN[0]-1)
2748                                         +(i1+1)/2*m_NN[0]*(m_NN[1]-1)                                          +(i1+1)/2*m_NN[0]*(m_NN[1]-1)
2749                                         +(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);
2750                          offsetInShared.push_back(offsetInShared.back()+1);                          offsetInShared.push_back(offsetInShared.back()+1);
2751                          sendShared.push_back(dof);                          sendShared.push_back(dof);
2752                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
2753                          doublyLink(colIndices, rowIndices, dof, numShared);                          m_dofMap[node] = numDOF+numShared;
                         m_dofMap[node]=numDOF+numShared;  
2754                          ++numShared;                          ++numShared;
2755                      }                      }
2756                  }                  }
# Line 2596  void Brick::createPattern() Line 2758  void Brick::createPattern()
2758          }          }
2759      }      }
2760    
2761  #pragma omp parallel for      // TODO: paso::SharedComponents should take vectors to avoid this
2762      for (int i = 0; i < numShared; i++) {      Esys_MPI_rank* neighPtr = NULL;
2763          std::sort(rowIndices[i].begin(), rowIndices[i].end());      index_t* sendPtr = NULL;
2764        index_t* recvPtr = NULL;
2765        if (neighbour.size() > 0) {
2766            neighPtr = &neighbour[0];
2767            sendPtr = &sendShared[0];
2768            recvPtr = &recvShared[0];
2769      }      }
   
2770      // create connector      // create connector
2771      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(      paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2772              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), neighPtr, sendPtr,
2773              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2774      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(      paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2775              numDOF, neighbour.size(), &neighbour[0], &recvShared[0],              numDOF, neighbour.size(), neighPtr, recvPtr,
2776              &offsetInShared[0], 1, 0, m_mpiInfo));              &offsetInShared[0], 1, 0, m_mpiInfo));
2777      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));      m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2778    
# Line 2629  void Brick::createPattern() Line 2795  void Brick::createPattern()
2795      for (size_t i=0; i<m_dofMap.size(); i++) {      for (size_t i=0; i<m_dofMap.size(); i++) {
2796          cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;          cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
2797      }      }
     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;  
     }  
2798      */      */
2799  }  }
2800    
2801  //private  //private
2802  void Brick::addToMatrixAndRHS(SystemMatrix* S, escript::Data& F,  void Brick::addToMatrixAndRHS(AbstractSystemMatrix* S, escript::Data& F,
2803           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,           const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2804           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const           bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2805  {  {
2806      IndexVector rowIndex;      IndexVector rowIndex(8);
2807      rowIndex.push_back(m_dofMap[firstNode]);      rowIndex[0] = m_dofMap[firstNode];
2808      rowIndex.push_back(m_dofMap[firstNode+1]);      rowIndex[1] = m_dofMap[firstNode+1];
2809      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);      rowIndex[2] = m_dofMap[firstNode+m_NN[0]];
2810      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);      rowIndex[3] = m_dofMap[firstNode+m_NN[0]+1];
2811      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*m_NN[1]]);      rowIndex[4] = m_dofMap[firstNode+m_NN[0]*m_NN[1]];
2812      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*m_NN[1]+1]);      rowIndex[5] = m_dofMap[firstNode+m_NN[0]*m_NN[1]+1];
2813      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)]);      rowIndex[6] = m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)];
2814      rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)+1]);      rowIndex[7] = m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)+1];
2815      if (addF) {      if (addF) {
2816          double *F_p=F.getSampleDataRW(0);          double *F_p=F.getSampleDataRW(0);
2817          for (index_t i=0; i<rowIndex.size(); i++) {          for (index_t i=0; i<rowIndex.size(); i++) {
# Line 2694  void Brick::addToMatrixAndRHS(SystemMatr Line 2823  void Brick::addToMatrixAndRHS(SystemMatr
2823          }          }
2824      }      }
2825      if (addS) {      if (addS) {
2826          S->add(rowIndex, EM_S);          addToSystemMatrix(S, rowIndex, nEq, EM_S);
2827      }      }
2828  }  }
2829    
# Line 3024  namespace Line 3153  namespace
3153      double* get3DGauss(unsigned radius, double sigma)      double* get3DGauss(unsigned radius, double sigma)
3154      {      {
3155          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)];
3156          double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);          double common = pow(M_1_PI * 0.5 / (sigma*sigma), 3./2);
3157          double total=0;          double total=0;
3158          int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
3159          for (int z=-r;z<=r;++z)          for (int z=-r;z<=r;++z) {
3160          {              for (int y=-r;y<=r;++y) {
3161              for (int y=-r;y<=r;++y)                  for (int x=-r;x<=r;++x) {
             {  
                 for (int x=-r;x<=r;++x)  
                 {          
3162                      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));
3163                      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)];
3164                  }                  }
3165              }              }
3166          }          }
3167          double invtotal=1/total;          double invtotal=1/total;
3168          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) {
3169          {              arr[p]*=invtotal;
             arr[p]*=invtotal;  
3170          }          }
3171          return arr;          return arr;
3172      }      }
3173        
3174      // applies conv to source to get a point.      // applies conv to source to get a point.
3175      // (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
3176      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,
3177                          size_t zp, unsigned radius, size_t width, size_t height)
3178      {      {
3179          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;          size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3180          size_t sbase=bx+by*width+bz*width*height;          size_t sbase=bx+by*width+bz*width*height;
3181          double result=0;          double result=0;
3182          for (int z=0;z<2*radius+1;++z)          for (int z=0; z<2*radius+1; ++z) {
3183          {              for (int y=0; y<2*radius+1; ++y) {
3184              for (int y=0;y<2*radius+1;++y)                  for (int x=0; x<2*radius+1; ++x) {
3185              {                          result += conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
                 for (int x=0;x<2*radius+1;++x)  
                 {  
                     result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];  
3186                  }                  }
3187              }              }
3188          }          }
3189          // use this line for "pass-though" (return the centre point value)          // use this line for "pass-though" (return the centre point value)
3190  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];  //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3191          return result;                return result;
3192      }      }
3193  }  }
3194    
3195  /* This is a wrapper for filtered (and non-filtered) randoms  /* This is a wrapper for filtered (and non-filtered) randoms
3196   * For detailed doco see randomFillWorker   * For detailed doco see randomFillWorker
3197  */   */
3198  escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,  escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3199         const escript::FunctionSpace& what,                                  const escript::FunctionSpace& what,
3200         long seed, const boost::python::tuple& filter) const                                  long seed, const bp::tuple& filter) const
3201  {  {
3202      int numvals=escript::DataTypes::noValues(shape);      int numvals=escript::DataTypes::noValues(shape);
3203      if (len(filter)>0 && (numvals!=1))      if (len(filter) > 0 && numvals != 1) {
     {  
3204          throw RipleyException("Ripley only supports filters for scalar data.");          throw RipleyException("Ripley only supports filters for scalar data.");
3205      }      }
3206      escript::Data res=randomFillWorker(shape, seed, filter);      escript::Data res = randomFillWorker(shape, seed, filter);
3207      if (res.getFunctionSpace()!=what)      if (res.getFunctionSpace()!=what) {
3208      {          escript::Data r(res, what);
         escript::Data r=escript::Data(res, what);  
3209          return r;          return r;
3210      }      }
3211      return res;      return res;
# Line 3096  A parameter radius  gives the size of th Line 3217  A parameter radius  gives the size of th
3217  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
3218  in order to complete the stencil.  in order to complete the stencil.
3219    
3220  All local calculation is done on an array called `src`, with  All local calculation is done on an array called `src`, with
3221  dimensions = ext[0] * ext[1] *ext[2].  dimensions = ext[0] * ext[1] *ext[2].
3222  Where ext[i]= internal[i]+2*radius.  Where ext[i]= internal[i]+2*radius.
3223    
3224  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
3225  values themselves but also the external region.  values themselves but also the external region.
3226    
3227  In a hypothetical 1-D case:  In a hypothetical 1-D case:
# Line 3119  need to be known. Line 3240  need to be known.
3240    
3241  pp123(4)56   23(4)567pp  pp123(4)56   23(4)567pp
3242    
3243  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
3244  right hand rank.  right hand rank.
3245    
3246  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.
3247  inset=2*radius+1      inset=2*radius+1
3248  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
3249  that ripley has.  that ripley has.
3250  */  */
3251  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(
3252                            const escript::DataTypes::ShapeType& shape, long seed,
3253                            const bp::tuple& filter) const
3254  {  {
3255      unsigned int radius=0;  // these are only used by gaussian      unsigned int radius=0;  // these are only used by gaussian
3256      double sigma=0.5;      double sigma=0.5;
3257        
3258      unsigned int numvals=escript::DataTypes::noValues(shape);      unsigned int numvals=escript::DataTypes::noValues(shape);
3259        
3260      if (len(filter)==0)      if (len(filter)==0) {
     {  
3261      // nothing special required here yet      // nothing special required here yet
3262      }      } else if (len(filter) == 3) {
3263      else if (len(filter)==3)          bp::extract<string> ex(filter[0]);
3264      {          if (!ex.check() || (ex() != "gaussian")) {
         boost::python::extract<string> ex(filter[0]);  
         if (!ex.check() || (ex()!="gaussian"))  
         {  
3265              throw RipleyException("Unsupported random filter for Brick.");              throw RipleyException("Unsupported random filter for Brick.");
3266          }          }
3267          boost::python::extract<unsigned int> ex1(filter[1]);          bp::extract<unsigned int> ex1(filter[1]);
3268          if (!ex1.check())          if (!ex1.check()) {
         {  
3269              throw RipleyException("Radius of gaussian filter must be a positive integer.");              throw RipleyException("Radius of gaussian filter must be a positive integer.");
3270          }          }
3271          radius=ex1();          radius=ex1();
3272          sigma=0.5;          sigma=0.5;
3273          boost::python::extract<double> ex2(filter[2]);          bp::extract<double> ex2(filter[2]);
3274          if (!ex2.check() || (sigma=ex2())<=0)          if (!ex2.check() || (sigma=ex2()) <= 0) {
         {  
3275              throw RipleyException("Sigma must be a positive floating point number.");              throw RipleyException("Sigma must be a positive floating point number.");
3276          }                      }
3277      }      } else {
     else  
     {  
3278          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
3279      }      }
3280    
# Line 3170  escript::Data Brick::randomFillWorker(co Line 3285  escript::Data Brick::randomFillWorker(co
3285      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
3286      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing      ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3287      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing      ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3288        
3289      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3290      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3291    
3292      if (2*radius>=internal[0]-4)      if (2*radius>=internal[0]-4) {
     {  
3293          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3294      }      }
3295      if (2*radius>=internal[1]-4)      if (2*radius>=internal[1]-4) {
     {  
3296          throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");          throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3297      }      }
3298      if (2*radius>=internal[2]-4)      if (2*radius>=internal[2]-4) {
     {  
3299          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");
3300      }          }
3301      
3302      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3303      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);            esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);
       
 #ifdef ESYS_MPI  
     
     if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))  
     {  
     // since the dimensions are equal for all ranks, this exception  
     // will be thrown on all ranks  
     throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");  
3304    
3305    #ifdef ESYS_MPI
3306        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5)) {
3307            // since the dimensions are equal for all ranks, this exception
3308            // will be thrown on all ranks
3309            throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3310      }      }
3311      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3312      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];
3313      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3314  #endif      #endif
3315    
3316  /*      /*
3317      // if we wanted to test a repeating pattern      // if we wanted to test a repeating pattern
3318      size_t basex=0;      size_t basex=0;
3319      size_t basey=0;      size_t basey=0;
3320      size_t basez=0;      size_t basez=0;
3321  #ifdef ESYS_MPI      #ifdef ESYS_MPI
3322      basex=X*m_gNE[0]/m_NX[0];      basex=X*m_gNE[0]/m_NX[0];
3323      basey=Y*m_gNE[1]/m_NX[1];      basey=Y*m_gNE[1]/m_NX[1];
3324      basez=Z*m_gNE[2]/m_NX[2];      basez=Z*m_gNE[2]/m_NX[2];
3325        cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;
3326  cout << "basex=" << basex << " basey=" << basey << " basez=" << basez << endl;      #endif
       
 #endif      
3327      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);
3328  */  */
3329    
3330  #ifdef ESYS_MPI  #ifdef ESYS_MPI
   
   
   
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      // it's +2 not +1 because a whole element is shared (and hence there is
3333          // there is an overlap of two points both of which need to have "radius" points on either side.      // an overlap of two points both of which need to have "radius" points on
3334            // either side.
3335      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t inset=2*radius+2;
3336      size_t ymidlen=ext[1]-2*inset;    
3337        // how wide is the x-dimension between the two insets
3338        size_t xmidlen=ext[0]-2*inset;
3339        size_t ymidlen=ext[1]-2*inset;
3340      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3341        
3342      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);
3343        
3344      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
3345      MPI_Status stats[50];      MPI_Status stats[50];
3346      short rused=0;      short rused=0;
3347        
3348      messvec incoms;      messvec incoms;
3349      messvec outcoms;      messvec outcoms;
3350        
3351      grid.generateInNeighbours(X, Y, Z ,incoms);      grid.generateInNeighbours(X, Y, Z ,incoms);
3352      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3353        
       
3354      block.copyAllToBuffer(src);      block.copyAllToBuffer(src);
3355        
3356      int comserr=0;          int comserr=0;
3357      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0; i < incoms.size(); ++i) {
3358      {          message& m = incoms[i];
3359          message& m=incoms[i];          comserr |= MPI_Irecv(block.getInBuffer(m.destbuffid),
3360          comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));                             block.getBuffSize(m.destbuffid), MPI_DOUBLE,
3361                               m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3362          block.setUsed(m.destbuffid);          block.setUsed(m.destbuffid);
3363      }      }
3364    
3365      for (size_t i=0;i<outcoms.size();++i)      for (size_t i=0; i<outcoms.size(); ++i) {
     {  
3366          message& m=outcoms[i];          message& m=outcoms[i];
3367          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),
3368      }                                   block.getBuffSize(m.srcbuffid), MPI_DOUBLE,
3369                                     m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3370      if (!comserr)      }
3371      {      
3372        if (!comserr) {
3373          comserr=MPI_Waitall(rused, reqs, stats);          comserr=MPI_Waitall(rused, reqs, stats);
3374      }      }
3375    
3376      if (comserr)      if (comserr) {
3377      {          // Yes this is throwing an exception as a result of an MPI error.
3378      // Yes this is throwing an exception as a result of an MPI error.          // and no we don't inform the other ranks that we are doing this.
3379      // and no we don't inform the other ranks that we are doing this.          // however, we have no reason to believe coms work at this point anyway
3380      // however, we have no reason to believe coms work at this point anyway          throw RipleyException("Error in coms for randomFill");
         throw RipleyException("Error in coms for randomFill");        
3381      }      }
3382        
3383      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3384        #endif // ESYS_MPI
3385  #endif      
3386            // the truth of either should imply the truth of the other but let's be safe
3387      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) {
     {  
         
3388          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3389          escript::Data resdat(0, shape, fs , true);          escript::Data resdat(0, shape, fs , true);
3390          // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3391          escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3392        
       
3393          // now we need to copy values over          // now we need to copy values over
3394          for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0; z < internal[2]; ++z) {
3395          {              for (size_t y=0; y < internal[1]; ++y) {
3396              for (size_t y=0;y<(internal[1]);++y)                      for (size_t x=0; x < internal[0]; ++x) {
3397              {                      for (unsigned int i=0; i < numvals; ++i) {
                 for (size_t x=0;x<(internal[0]);++x)  
                 {  
                     for (unsigned int i=0;i<numvals;++i)  
                     {  
3398                          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];
3399                      }                      }
3400                  }                  }
3401              }              }
3402          }            }
3403          delete[] src;          delete[] src;
3404          return resdat;                return resdat;
3405      }      } else { // filter enabled
     else        // filter enabled    
     {  
       
3406          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());          escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3407          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);          escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3408          // don't need to check for exwrite because we just made it          // don't need to check for exwrite because we just made it
3409          escript::DataVector& dv=resdat.getExpandedVectorReference();          escript::DataVector& dv=resdat.getExpandedVectorReference();
3410          double* convolution=get3DGauss(radius, sigma);          double* convolution=get3DGauss(radius, sigma);
3411    
3412          for (size_t z=0;z<(internal[2]);++z)          for (size_t z=0;z<(internal[2]);++z) {
3413          {              for (size_t y=0;y<(internal[1]);++y) {
3414              for (size_t y=0;y<(internal[1]);++y)                      for (size_t x=0;x<(internal[0]);++x) {
             {  
                 for (size_t x=0;x<(internal[0]);++x)  
                 {      
3415                      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]);
3416                
3417                  }                  }
3418              }              }
3419          }          }
3420        
3421          delete[] convolution;          delete[] convolution;
3422          delete[] src;          delete[] src;
3423          return resdat;          return resdat;
       
3424      }      }
3425  }  }
3426    
3427  int Brick::findNode(const double *coords) const  dim_t Brick::findNode(const double *coords) const
3428  {  {
3429      const int NOT_MINE = -1;      const dim_t NOT_MINE = -1;
3430      //is the found element even owned by this rank      //is the found element even owned by this rank
3431      // (inside owned or shared elements but will map to an owned element)      // (inside owned or shared elements but will map to an owned element)
3432      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
# Line 3350  int Brick::findNode(const double *coords Line 3442  int Brick::findNode(const double *coords
3442      double x = coords[0] - m_origin[0];      double x = coords[0] - m_origin[0];
3443      double y = coords[1] - m_origin[1];      double y = coords[1] - m_origin[1];
3444      double z = coords[2] - m_origin[2];      double z = coords[2] - m_origin[2];
3445        
3446      //check if the point is even inside the domain      //check if the point is even inside the domain
3447      if (x < 0 || y < 0 || z < 0      if (x < 0 || y < 0 || z < 0
3448              || x > m_length[0] || y > m_length[1] || z > m_length[2])              || x > m_length[0] || y > m_length[1] || z > m_length[2])
3449          return NOT_MINE;          return NOT_MINE;
3450            
3451      // distance in elements      // distance in elements
3452      int ex = (int) floor(x / m_dx[0]);      dim_t ex = (dim_t) floor(x / m_dx[0]);
3453      int ey = (int) floor(y / m_dx[1]);      dim_t ey = (dim_t) floor(y / m_dx[1]);
3454      int ez = (int) floor(z / m_dx[2]);      dim_t ez = (dim_t) floor(z / m_dx[2]);
3455      // set the min distance high enough to be outside the element plus a bit      // set the min distance high enough to be outside the element plus a bit
3456      int closest = NOT_MINE;      dim_t closest = NOT_MINE;
3457      double minDist = 1;      double minDist = 1;
3458      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3459          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
# Line 3383  int Brick::findNode(const double *coords Line 3475  int Brick::findNode(const double *coords
3475          }          }
3476      }      }
3477      if (closest == NOT_MINE) {      if (closest == NOT_MINE) {
3478          throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");          throw RipleyException("Unable to map appropriate dirac point to a "
3479                             "node, implementation problem in Brick::findNode()");
3480      }      }
3481      return closest;      return closest;
3482  }  }
3483    
3484  Assembler_ptr Brick::createAssembler(std::string type, std::map<std::string,  Assembler_ptr Brick::createAssembler(string type, const DataMap& constants) const
         escript::Data> constants) const  
3485  {  {
3486      if (type.compare("DefaultAssembler") == 0) {      if (type.compare("DefaultAssembler") == 0) {
3487          return Assembler_ptr(new DefaultAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));          return Assembler_ptr(new DefaultAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
# Line 3398  Assembler_ptr Brick::createAssembler(std Line 3490  Assembler_ptr Brick::createAssembler(std
3490      } else if (type.compare("LameAssembler") == 0) {      } else if (type.compare("LameAssembler") == 0) {
3491          return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));          return Assembler_ptr(new LameAssembler3D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
3492      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
3493          throw RipleyException("Ripley::Brick does not support the"          throw RipleyException("Ripley::Brick does not support the requested "
3494                                  " requested assembler");                                "assembler");
3495      }      }
3496  }  }
3497    

Legend:
Removed from v.5084  
changed lines
  Added in v.5122

  ViewVC Help
Powered by ViewVC 1.1.26