/[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 4988 by caltinay, Wed Jun 4 01:15:10 2014 UTC revision 5136 by caltinay, Tue Sep 9 07:13:55 2014 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 #include <limits>  
   
17  #include <ripley/Brick.h>  #include <ripley/Brick.h>
 #include <esysUtils/esysFileWriter.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>
22  #include <ripley/domainhelpers.h>  #include <ripley/domainhelpers.h>
23    #include <esysUtils/esysFileWriter.h>
24    #include <esysUtils/EsysRandom.h>
25    #include <paso/SystemMatrix.h>
26    
27  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
28    
29  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 36  Line 38 
38  #endif  #endif
39    
40  #include <iomanip>  #include <iomanip>
41    #include <limits>
42    
43  #include "esysUtils/EsysRandom.h"  namespace bp = boost::python;
 #include "blocktools.h"  
   
   
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 TagMap& 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 (TagMap::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 531  void Brick::readBinaryGridImpl(escript:: Line 557  void Brick::readBinaryGridImpl(escript::
557                                  if (params.byteOrder != BYTEORDER_NATIVE) {                                  if (params.byteOrder != BYTEORDER_NATIVE) {
558                                      char* cval = reinterpret_cast<char*>(&val);                                      char* cval = reinterpret_cast<char*>(&val);
559                                      // this will alter val!!                                      // this will alter val!!
560                                      byte_swap32(cval);                                      if (sizeof(ValueType)>4) {
561                                            byte_swap64(cval);
562                                        } else {
563                                            byte_swap32(cval);
564                                        }
565                                  }                                  }
566                                  if (!std::isnan(val)) {                                  if (!isnan(val)) {
567                                      for (int q=0; q<dpp; q++) {                                      for (int q=0; q<dpp; q++) {
568                                          *dest++ = static_cast<double>(val);                                          *dest++ = static_cast<double>(val);
569                                      }                                      }
# Line 555  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 587  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 612  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 653  void Brick::readBinaryGridZippedImpl(esc Line 683  void Brick::readBinaryGridZippedImpl(esc
683                                      // this will alter val!!                                      // this will alter val!!
684                                      byte_swap32(cval);                                      byte_swap32(cval);
685                                  }                                  }
686                                  if (!std::isnan(val)) {                                  if (!isnan(val)) {
687                                      for (int q=0; q<dpp; q++) {                                      for (int q=0; q<dpp; q++) {
688                                          *dest++ = static_cast<double>(val);                                          *dest++ = static_cast<double>(val);
689                                      }                                      }
# Line 693  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 715  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 727  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 738  void Brick::writeBinaryGridImpl(const es Line 768  void Brick::writeBinaryGridImpl(const es
768                      oss.write((char*)&fvalue, sizeof(fvalue));                      oss.write((char*)&fvalue, sizeof(fvalue));
769                  } else {                  } else {
770                      char* value = reinterpret_cast<char*>(&fvalue);                      char* value = reinterpret_cast<char*>(&fvalue);
771                      oss.write(byte_swap32(value), sizeof(fvalue));                      if (sizeof(fvalue)>4) {
772                            byte_swap64(value);
773                        } else {
774                            byte_swap32(value);
775                        }
776                        oss.write(value, sizeof(fvalue));
777                  }                  }
778              }              }
779              fw.writeAt(oss, fileofs);              fw.writeAt(oss, fileofs);
# Line 747  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 755  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 811  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 dim_t NN0 = m_NN[0];
855        const dim_t NN1 = m_NN[1];
856        const dim_t NN2 = m_NN[2];
857    
858  #pragma omp parallel  #pragma omp parallel
859      {      {
860  #pragma omp for  #pragma omp for
861          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {          for (dim_t i0 = 0; i0 < NN0; i0++) {
862              coords[0][i0]=getLocalCoordinate(i0, 0);              coords[0][i0]=getLocalCoordinate(i0, 0);
863          }          }
864  #pragma omp for  #pragma omp for
865          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {          for (dim_t i1 = 0; i1 < NN1; i1++) {
866              coords[1][i1]=getLocalCoordinate(i1, 1);              coords[1][i1]=getLocalCoordinate(i1, 1);
867          }          }
868  #pragma omp for  #pragma omp for
869          for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {          for (dim_t i2 = 0; i2 < NN2; i2++) {
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 837  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 892  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 973  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 dim_t NE0 = m_NE[0];
1021        const dim_t NE1 = m_NE[1];
1022        const dim_t NE2 = m_NE[2];
1023    
1024      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1025          out.requireWrite();          out.requireWrite();
1026  #pragma omp parallel  #pragma omp parallel
1027          {          {
1028              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1029  #pragma omp for nowait  #pragma omp for nowait
1030                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1031                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1032                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1033                          // set vector at four quadrature points                          // set vector at four quadrature points
1034                          *o++ = -1.; *o++ = 0.; *o++ = 0.;                          *o++ = -1.; *o++ = 0.; *o++ = 0.;
# Line 993  void Brick::setToNormal(escript::Data& o Line 1041  void Brick::setToNormal(escript::Data& o
1041    
1042              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1043  #pragma omp for nowait  #pragma omp for nowait
1044                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1045                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1046                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1047                          // set vector at four quadrature points                          // set vector at four quadrature points
1048                          *o++ = 1.; *o++ = 0.; *o++ = 0.;                          *o++ = 1.; *o++ = 0.; *o++ = 0.;
# Line 1007  void Brick::setToNormal(escript::Data& o Line 1055  void Brick::setToNormal(escript::Data& o
1055    
1056              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1057  #pragma omp for nowait  #pragma omp for nowait
1058                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1059                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1060                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1061                          // set vector at four quadrature points                          // set vector at four quadrature points
1062                          *o++ = 0.; *o++ = -1.; *o++ = 0.;                          *o++ = 0.; *o++ = -1.; *o++ = 0.;
# Line 1021  void Brick::setToNormal(escript::Data& o Line 1069  void Brick::setToNormal(escript::Data& o
1069    
1070              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1071  #pragma omp for nowait  #pragma omp for nowait
1072                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1073                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1074                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1075                          // set vector at four quadrature points                          // set vector at four quadrature points
1076                          *o++ = 0.; *o++ = 1.; *o++ = 0.;                          *o++ = 0.; *o++ = 1.; *o++ = 0.;
# Line 1035  void Brick::setToNormal(escript::Data& o Line 1083  void Brick::setToNormal(escript::Data& o
1083    
1084              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1085  #pragma omp for nowait  #pragma omp for nowait
1086                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1087                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1088                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1089                          // set vector at four quadrature points                          // set vector at four quadrature points
1090                          *o++ = 0.; *o++ = 0.; *o++ = -1.;                          *o++ = 0.; *o++ = 0.; *o++ = -1.;
# Line 1049  void Brick::setToNormal(escript::Data& o Line 1097  void Brick::setToNormal(escript::Data& o
1097    
1098              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1099  #pragma omp for nowait  #pragma omp for nowait
1100                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1101                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1102                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1103                          // set vector at four quadrature points                          // set vector at four quadrature points
1104                          *o++ = 0.; *o++ = 0.; *o++ = 1.;                          *o++ = 0.; *o++ = 0.; *o++ = 1.;
# Line 1067  void Brick::setToNormal(escript::Data& o Line 1115  void Brick::setToNormal(escript::Data& o
1115          {          {
1116              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1117  #pragma omp for nowait  #pragma omp for nowait
1118                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1119                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1120                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1121                          *o++ = -1.;                          *o++ = -1.;
1122                          *o++ = 0.;                          *o++ = 0.;
# Line 1079  void Brick::setToNormal(escript::Data& o Line 1127  void Brick::setToNormal(escript::Data& o
1127    
1128              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1129  #pragma omp for nowait  #pragma omp for nowait
1130                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1131                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1132                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1133                          *o++ = 1.;                          *o++ = 1.;
1134                          *o++ = 0.;                          *o++ = 0.;
# Line 1091  void Brick::setToNormal(escript::Data& o Line 1139  void Brick::setToNormal(escript::Data& o
1139    
1140              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1141  #pragma omp for nowait  #pragma omp for nowait
1142                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1143                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1144                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1145                          *o++ = 0.;                          *o++ = 0.;
1146                          *o++ = -1.;                          *o++ = -1.;
# Line 1103  void Brick::setToNormal(escript::Data& o Line 1151  void Brick::setToNormal(escript::Data& o
1151    
1152              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1153  #pragma omp for nowait  #pragma omp for nowait
1154                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1155                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1156                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1157                          *o++ = 0.;                          *o++ = 0.;
1158                          *o++ = 1.;                          *o++ = 1.;
# Line 1115  void Brick::setToNormal(escript::Data& o Line 1163  void Brick::setToNormal(escript::Data& o
1163    
1164              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1165  #pragma omp for nowait  #pragma omp for nowait
1166                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1167                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1168                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1169                          *o++ = 0.;                          *o++ = 0.;
1170                          *o++ = 0.;                          *o++ = 0.;
# Line 1127  void Brick::setToNormal(escript::Data& o Line 1175  void Brick::setToNormal(escript::Data& o
1175    
1176              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1177  #pragma omp for nowait  #pragma omp for nowait
1178                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1179                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1180                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1181                          *o++ = 0.;                          *o++ = 0.;
1182                          *o++ = 0.;                          *o++ = 0.;
# Line 1153  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 dim_t NE = getNumElements();
1205  #pragma omp parallel for  #pragma omp parallel for
1206          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < NE; ++k) {
1207              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
1208              fill(o, o+numQuad, size);              fill(o, o+numQuad, size);
1209          }          }
# Line 1162  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 dim_t NE0 = m_NE[0];
1215            const dim_t NE1 = m_NE[1];
1216            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) {
1220                  const double size=min(m_dx[1],m_dx[2]);                  const double size=min(m_dx[1],m_dx[2]);
1221  #pragma omp for nowait  #pragma omp for nowait
1222                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1223                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1224                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1225                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1226                      }                      }
# Line 1178  void Brick::setToSize(escript::Data& out Line 1230  void Brick::setToSize(escript::Data& out
1230              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1231                  const double size=min(m_dx[1],m_dx[2]);                  const double size=min(m_dx[1],m_dx[2]);
1232  #pragma omp for nowait  #pragma omp for nowait
1233                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1234                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1235                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1236                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1237                      }                      }
# Line 1189  void Brick::setToSize(escript::Data& out Line 1241  void Brick::setToSize(escript::Data& out
1241              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1242                  const double size=min(m_dx[0],m_dx[2]);                  const double size=min(m_dx[0],m_dx[2]);
1243  #pragma omp for nowait  #pragma omp for nowait
1244                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1245                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1246                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1247                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1248                      }                      }
# Line 1200  void Brick::setToSize(escript::Data& out Line 1252  void Brick::setToSize(escript::Data& out
1252              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1253                  const double size=min(m_dx[0],m_dx[2]);                  const double size=min(m_dx[0],m_dx[2]);
1254  #pragma omp for nowait  #pragma omp for nowait
1255                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1256                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1257                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1258                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1259                      }                      }
# Line 1211  void Brick::setToSize(escript::Data& out Line 1263  void Brick::setToSize(escript::Data& out
1263              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1264                  const double size=min(m_dx[0],m_dx[1]);                  const double size=min(m_dx[0],m_dx[1]);
1265  #pragma omp for nowait  #pragma omp for nowait
1266                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1267                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1268                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1269                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1270                      }                      }
# Line 1222  void Brick::setToSize(escript::Data& out Line 1274  void Brick::setToSize(escript::Data& out
1274              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1275                  const double size=min(m_dx[0],m_dx[1]);                  const double size=min(m_dx[0],m_dx[1]);
1276  #pragma omp for nowait  #pragma omp for nowait
1277                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1278                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1279                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1280                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1281                      }                      }
# Line 1266  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 dim_t NN0 = m_NN[0];
1322        const dim_t NN1 = m_NN[1];
1323        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 < m_NN[2]; i2++) {      for (dim_t i2 = 0; i2 < NN2; i2++) {
1327          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {          for (dim_t i1 = 0; i1 < NN1; i1++) {
1328              for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {              for (dim_t i0 = 0; i0 < NN0; i0++) {
1329                  double* point = arg.getSampleDataRW(i0+m_NN[0]*i1+m_NN[0]*m_NN[1]*i2);                  double* point = arg.getSampleDataRW(i0+NN0*i1+NN0*NN1*i2);
1330                  point[0] = getLocalCoordinate(i0, 0);                  point[0] = getLocalCoordinate(i0, 0);
1331                  point[1] = getLocalCoordinate(i1, 1);                  point[1] = getLocalCoordinate(i1, 1);
1332                  point[2] = getLocalCoordinate(i2, 2);                  point[2] = getLocalCoordinate(i2, 2);
# Line 1291  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 dim_t NE0 = m_NE[0];
1350        const dim_t NE1 = m_NE[1];
1351        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 1305  void Brick::assembleGradient(escript::Da Line 1363  void Brick::assembleGradient(escript::Da
1363              vector<double> f_110(numComp);              vector<double> f_110(numComp);
1364              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1365  #pragma omp for  #pragma omp for
1366              for (index_t k2=0; k2 < m_NE[2]; ++k2) {              for (index_t k2=0; k2 < NE2; ++k2) {
1367                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1368                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1369                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1370                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1371                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1316  void Brick::assembleGradient(escript::Da Line 1374  void Brick::assembleGradient(escript::Da
1374                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1375                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1376                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1377                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,NE0,NE1));
1378                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1379                              const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];
1380                              const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];                              const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];
# Line 1372  void Brick::assembleGradient(escript::Da Line 1430  void Brick::assembleGradient(escript::Da
1430              vector<double> f_110(numComp);              vector<double> f_110(numComp);
1431              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1432  #pragma omp for  #pragma omp for
1433              for (index_t k2=0; k2 < m_NE[2]; ++k2) {              for (index_t k2=0; k2 < NE2; ++k2) {
1434                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1435                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1436                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1437                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1438                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1383  void Brick::assembleGradient(escript::Da Line 1441  void Brick::assembleGradient(escript::Da
1441                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1442                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1443                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1444                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,NE0,NE1));
1445                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1446                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1447                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1407  void Brick::assembleGradient(escript::Da Line 1465  void Brick::assembleGradient(escript::Da
1465              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1466              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1467  #pragma omp for nowait  #pragma omp for nowait
1468                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1469                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1470                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1471                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1472                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1417  void Brick::assembleGradient(escript::Da Line 1475  void Brick::assembleGradient(escript::Da
1475                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1476                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1477                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1478                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,NE1));
1479                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1480                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];
1481                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];
# Line 1441  void Brick::assembleGradient(escript::Da Line 1499  void Brick::assembleGradient(escript::Da
1499              } // end of face 0              } // end of face 0
1500              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1501  #pragma omp for nowait  #pragma omp for nowait
1502                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1503                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1504                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1505                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1506                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1451  void Brick::assembleGradient(escript::Da Line 1509  void Brick::assembleGradient(escript::Da
1509                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1510                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1511                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1512                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,NE1));
1513                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1514                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1515                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
# Line 1475  void Brick::assembleGradient(escript::Da Line 1533  void Brick::assembleGradient(escript::Da
1533              } // end of face 1              } // end of face 1
1534              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1535  #pragma omp for nowait  #pragma omp for nowait
1536                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1537                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1538                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1539                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1540                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1485  void Brick::assembleGradient(escript::Da Line 1543  void Brick::assembleGradient(escript::Da
1543                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1544                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1545                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1546                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,NE0));
1547                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1548                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];
1549                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];
# Line 1508  void Brick::assembleGradient(escript::Da Line 1566  void Brick::assembleGradient(escript::Da
1566              } // end of face 2              } // end of face 2
1567              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1568  #pragma omp for nowait  #pragma omp for nowait
1569                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1570                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1571                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1572                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1573                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1518  void Brick::assembleGradient(escript::Da Line 1576  void Brick::assembleGradient(escript::Da
1576                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1577                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1578                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1579                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,NE0));
1580                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1581                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1582                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
# Line 1542  void Brick::assembleGradient(escript::Da Line 1600  void Brick::assembleGradient(escript::Da
1600              } // end of face 3              } // end of face 3
1601              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1602  #pragma omp for nowait  #pragma omp for nowait
1603                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1604                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1605                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1606                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1607                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1552  void Brick::assembleGradient(escript::Da Line 1610  void Brick::assembleGradient(escript::Da
1610                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1611                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1612                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1613                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,NE0));
1614                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1615                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];
1616                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];
# Line 1576  void Brick::assembleGradient(escript::Da Line 1634  void Brick::assembleGradient(escript::Da
1634              } // end of face 4              } // end of face 4
1635              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1636  #pragma omp for nowait  #pragma omp for nowait
1637                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1638                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1639                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1640                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1641                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1586  void Brick::assembleGradient(escript::Da Line 1644  void Brick::assembleGradient(escript::Da
1644                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1645                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1646                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1647                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,NE0));
1648                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1649                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1650                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
# Line 1623  void Brick::assembleGradient(escript::Da Line 1681  void Brick::assembleGradient(escript::Da
1681              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1682              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1683  #pragma omp for nowait  #pragma omp for nowait
1684                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1685                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1686                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1687                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1688                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1633  void Brick::assembleGradient(escript::Da Line 1691  void Brick::assembleGradient(escript::Da
1691                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1692                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1693                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1694                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,NE1));
1695                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1696                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1697                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];
# Line 1644  void Brick::assembleGradient(escript::Da Line 1702  void Brick::assembleGradient(escript::Da
1702              } // end of face 0              } // end of face 0
1703              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1704  #pragma omp for nowait  #pragma omp for nowait
1705                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1706                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1707                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1708                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1709                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1654  void Brick::assembleGradient(escript::Da Line 1712  void Brick::assembleGradient(escript::Da
1712                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1713                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1714                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1715                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,NE1));
1716                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1717                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1718                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];
# Line 1665  void Brick::assembleGradient(escript::Da Line 1723  void Brick::assembleGradient(escript::Da
1723              } // end of face 1              } // end of face 1
1724              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1725  #pragma omp for nowait  #pragma omp for nowait
1726                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1727                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1728                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1729                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1730                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1675  void Brick::assembleGradient(escript::Da Line 1733  void Brick::assembleGradient(escript::Da
1733                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1734                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1735                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1736                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,NE0));
1737                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1738                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];
1739                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1686  void Brick::assembleGradient(escript::Da Line 1744  void Brick::assembleGradient(escript::Da
1744              } // end of face 2              } // end of face 2
1745              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1746  #pragma omp for nowait  #pragma omp for nowait
1747                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1748                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1749                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1750                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1751                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1696  void Brick::assembleGradient(escript::Da Line 1754  void Brick::assembleGradient(escript::Da
1754                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1755                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1756                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1757                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,NE0));
1758                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1759                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];
1760                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1707  void Brick::assembleGradient(escript::Da Line 1765  void Brick::assembleGradient(escript::Da
1765              } // end of face 3              } // end of face 3
1766              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1767  #pragma omp for nowait  #pragma omp for nowait
1768                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1769                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1770                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1771                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1772                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1717  void Brick::assembleGradient(escript::Da Line 1775  void Brick::assembleGradient(escript::Da
1775                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1776                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1777                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1778                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,NE0));
1779                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1780                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];
1781                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];
# Line 1728  void Brick::assembleGradient(escript::Da Line 1786  void Brick::assembleGradient(escript::Da
1786              } // end of face 4              } // end of face 4
1787              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1788  #pragma omp for nowait  #pragma omp for nowait
1789                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1790                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1791                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1792                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1793                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1738  void Brick::assembleGradient(escript::Da Line 1796  void Brick::assembleGradient(escript::Da
1796                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1797                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1798                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1799                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,NE0));
1800                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1801                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];
1802                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];
# Line 2056  void Brick::dofToNodes(escript::Data& ou Line 2114  void Brick::dofToNodes(escript::Data& ou
2114      coupler->startCollect(in.getDataRO());      coupler->startCollect(in.getDataRO());
2115    
2116      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2117        const dim_t numNodes = getNumNodes();
2118      out.requireWrite();      out.requireWrite();
2119      const double* buffer = coupler->finishCollect();      const double* buffer = coupler->finishCollect();
2120    
2121  #pragma omp parallel for  #pragma omp parallel for
2122      for (index_t i=0; i<getNumNodes(); i++) {      for (index_t i=0; i<numNodes; i++) {
2123          const double* src=(m_dofMap[i]<numDOF ?          const double* src=(m_dofMap[i]<numDOF ?
2124                  in.getSampleDataRO(m_dofMap[i])                  in.getSampleDataRO(m_dofMap[i])
2125                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
# Line 2068  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 2086  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 2127  void Brick::populateSampleIds() Line 2389  void Brick::populateSampleIds()
2389      else      else
2390          m_faceCount[5]=0;          m_faceCount[5]=0;
2391    
2392      m_faceId.resize(getNumFaceElements());      const dim_t NFE = getNumFaceElements();
2393        m_faceId.resize(NFE);
2394    
2395      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
2396      const index_t bottom = (m_offset[1]==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
# Line 2135  void Brick::populateSampleIds() Line 2398  void Brick::populateSampleIds()
2398      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2399      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2400      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2401        const dim_t NN0 = m_NN[0];
2402        const dim_t NN1 = m_NN[1];
2403        const dim_t NN2 = m_NN[2];
2404        const dim_t NE0 = m_NE[0];
2405        const dim_t NE1 = m_NE[1];
2406        const dim_t NE2 = m_NE[2];
2407    
2408      // the following is a compromise between efficiency and code length to      // the following is a compromise between efficiency and code length to
2409      // set the node id's according to the order mentioned above.      // set the node id's according to the order mentioned above.
# Line 2153  void Brick::populateSampleIds() Line 2422  void Brick::populateSampleIds()
2422          // set edge id's          // set edge id's
2423          // edges in x-direction, including corners          // edges in x-direction, including corners
2424  #pragma omp for nowait  #pragma omp for nowait
2425          for (dim_t i=0; i<m_NN[0]; i++) {          for (dim_t i=0; i<NN0; i++) {
2426              m_nodeId[i] = globalNodeId(i, 0, 0); // LF              m_nodeId[i] = globalNodeId(i, 0, 0); // LF
2427              m_nodeId[m_NN[0]*(m_NN[1]-1)+i] = globalNodeId(i, m_NN[1]-1, 0); // UF              m_nodeId[NN0*(NN1-1)+i] = globalNodeId(i, NN1-1, 0); // UF
2428              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+i] = globalNodeId(i, 0, m_NN[2]-1); // LB              m_nodeId[NN0*NN1*(NN2-1)+i] = globalNodeId(i, 0, NN2-1); // LB
2429              m_nodeId[m_NN[0]*m_NN[1]*m_NN[2]-m_NN[0]+i] = globalNodeId(i, m_NN[1]-1, m_NN[2]-1); // UB              m_nodeId[NN0*NN1*NN2-NN0+i] = globalNodeId(i, NN1-1, NN2-1); // UB
2430          }          }
2431          // edges in y-direction, without corners          // edges in y-direction, without corners
2432  #pragma omp for nowait  #pragma omp for nowait
2433          for (dim_t i=1; i<m_NN[1]-1; i++) {          for (dim_t i=1; i<NN1-1; i++) {
2434              m_nodeId[m_NN[0]*i] = globalNodeId(0, i, 0); // FL              m_nodeId[NN0*i] = globalNodeId(0, i, 0); // FL
2435              m_nodeId[m_NN[0]*(i+1)-1] = globalNodeId(m_NN[0]-1, i, 0); // FR              m_nodeId[NN0*(i+1)-1] = globalNodeId(NN0-1, i, 0); // FR
2436              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+m_NN[0]*i] = globalNodeId(0, i, m_NN[2]-1); // BL              m_nodeId[NN0*NN1*(NN2-1)+NN0*i] = globalNodeId(0, i, NN2-1); // BL
2437              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+m_NN[0]*(i+1)-1] = globalNodeId(m_NN[0]-1, i, m_NN[2]-1); // BR              m_nodeId[NN0*NN1*(NN2-1)+NN0*(i+1)-1] = globalNodeId(NN0-1, i, NN2-1); // BR
2438          }          }
2439          // edges in z-direction, without corners          // edges in z-direction, without corners
2440  #pragma omp for  #pragma omp for
2441          for (dim_t i=1; i<m_NN[2]-1; i++) {          for (dim_t i=1; i<NN2-1; i++) {
2442              m_nodeId[m_NN[0]*m_NN[1]*i] = globalNodeId(0, 0, i); // LL              m_nodeId[NN0*NN1*i] = globalNodeId(0, 0, i); // LL
2443              m_nodeId[m_NN[0]*m_NN[1]*i+m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0, i); // LR              m_nodeId[NN0*NN1*i+NN0-1] = globalNodeId(NN0-1, 0, i); // LR
2444              m_nodeId[m_NN[0]*m_NN[1]*(i+1)-m_NN[0]] = globalNodeId(0, m_NN[1]-1, i); // UL              m_nodeId[NN0*NN1*(i+1)-NN0] = globalNodeId(0, NN1-1, i); // UL
2445              m_nodeId[m_NN[0]*m_NN[1]*(i+1)-1] = globalNodeId(m_NN[0]-1, m_NN[1]-1, i); // UR              m_nodeId[NN0*NN1*(i+1)-1] = globalNodeId(NN0-1, NN1-1, i); // UR
2446          }          }
2447          // implicit barrier here because some node IDs will be overwritten          // implicit barrier here because some node IDs will be overwritten
2448          // below          // below
# Line 2183  void Brick::populateSampleIds() Line 2452  void Brick::populateSampleIds()
2452          for (dim_t i=0; i<nDOF2; i++) {          for (dim_t i=0; i<nDOF2; i++) {
2453              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2454                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2455                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=k+left+(j+bottom)*NN0+(i+front)*NN0*NN1;
2456                      const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;                      const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;
2457                      m_dofId[dofIdx] = m_nodeId[nodeIdx]                      m_dofId[dofIdx] = m_nodeId[nodeIdx]
2458                          = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;                          = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
# Line 2196  void Brick::populateSampleIds() Line 2465  void Brick::populateSampleIds()
2465  #pragma omp for nowait  #pragma omp for nowait
2466              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2467                  for (dim_t j=0; j<nDOF1; j++) {                  for (dim_t j=0; j<nDOF1; j++) {
2468                      const index_t nodeIdx=(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=(j+bottom)*NN0+(i+front)*NN0*NN1;
2469                      const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;                      const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;
2470                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2471                          = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
# Line 2207  void Brick::populateSampleIds() Line 2476  void Brick::populateSampleIds()
2476  #pragma omp for nowait  #pragma omp for nowait
2477              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2478                  for (dim_t j=0; j<nDOF1; j++) {                  for (dim_t j=0; j<nDOF1; j++) {
2479                      const index_t nodeIdx=(j+bottom+1)*m_NN[0]-1+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=(j+bottom+1)*NN0-1+(i+front)*NN0*NN1;
2480                      const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;                      const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;
2481                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2482                          = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
# Line 2218  void Brick::populateSampleIds() Line 2487  void Brick::populateSampleIds()
2487  #pragma omp for nowait  #pragma omp for nowait
2488              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2489                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2490                      const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=k+left+(i+front)*NN0*NN1;
2491                      const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;                      const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;
2492                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2493                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
# Line 2229  void Brick::populateSampleIds() Line 2498  void Brick::populateSampleIds()
2498  #pragma omp for nowait  #pragma omp for nowait
2499              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2500                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2501                      const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1]+m_NN[0]*(m_NN[1]-1);                      const index_t nodeIdx=k+left+(i+front)*NN0*NN1+NN0*(NN1-1);
2502                      const index_t dofId=k+i*nDOF0*nDOF1;                      const index_t dofId=k+i*nDOF0*nDOF1;
2503                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2504                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
# Line 2240  void Brick::populateSampleIds() Line 2509  void Brick::populateSampleIds()
2509  #pragma omp for nowait  #pragma omp for nowait
2510              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2511                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2512                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0];                      const index_t nodeIdx=k+left+(j+bottom)*NN0;
2513                      const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);                      const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);
2514                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2515                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;
# Line 2251  void Brick::populateSampleIds() Line 2520  void Brick::populateSampleIds()
2520  #pragma omp for nowait  #pragma omp for nowait
2521              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2522                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2523                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1);                      const index_t nodeIdx=k+left+(j+bottom)*NN0+NN0*NN1*(NN2-1);
2524                      const index_t dofId=k+j*nDOF0;                      const index_t dofId=k+j*nDOF0;
2525                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2526                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;
# Line 2261  void Brick::populateSampleIds() Line 2530  void Brick::populateSampleIds()
2530    
2531          // populate element id's          // populate element id's
2532  #pragma omp for nowait  #pragma omp for nowait
2533          for (dim_t i2=0; i2<m_NE[2]; i2++) {          for (dim_t i2=0; i2<NE2; i2++) {
2534              for (dim_t i1=0; i1<m_NE[1]; i1++) {              for (dim_t i1=0; i1<NE1; i1++) {
2535                  for (dim_t i0=0; i0<m_NE[0]; i0++) {                  for (dim_t i0=0; i0<NE0; i0++) {
2536                      m_elementId[i0+i1*m_NE[0]+i2*m_NE[0]*m_NE[1]] =                      m_elementId[i0+i1*NE0+i2*NE0*NE1] =
2537                          (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]                          (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]
2538                          +(m_offset[1]+i1)*m_gNE[0]                          +(m_offset[1]+i1)*m_gNE[0]
2539                          +m_offset[0]+i0;                          +m_offset[0]+i0;
# Line 2274  void Brick::populateSampleIds() Line 2543  void Brick::populateSampleIds()
2543    
2544          // face elements          // face elements
2545  #pragma omp for  #pragma omp for
2546          for (dim_t k=0; k<getNumFaceElements(); k++)          for (dim_t k=0; k<NFE; k++)
2547              m_faceId[k]=k;              m_faceId[k]=k;
2548      } // end parallel section      } // end parallel section
2549    
# Line 2306  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  //private
2583  void Brick::createPattern()  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
2615    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 2330  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 2391  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 2424  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 2458  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 2561  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 2594  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 2659  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 2989  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 3061  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 3084  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 3135  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 3315  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 3348  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 3363  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.4988  
changed lines
  Added in v.5136

  ViewVC Help
Powered by ViewVC 1.1.26