/[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

branches/ripleygmg_from_3668/ripley/src/Brick.cpp revision 3702 by caltinay, Fri Dec 2 06:12:32 2011 UTC trunk/ripley/src/Brick.cpp revision 4861 by sshaw, Thu Apr 10 05:17:47 2014 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14    *
15    *****************************************************************************/
16    
17    #include <limits>
18    
19  #include <ripley/Brick.h>  #include <ripley/Brick.h>
20  extern "C" {  #include <paso/SystemMatrix.h>
21  #include "paso/SystemMatrixPattern.h"  #include <esysUtils/esysFileWriter.h>
22  }  #include <ripley/DefaultAssembler3D.h>
23    #include <ripley/WaveAssembler3D.h>
24    #include <ripley/LameAssembler3D.h>
25    #include <ripley/domainhelpers.h>
26    #include <boost/scoped_array.hpp>
27    
28    #ifdef USE_NETCDF
29    #include <netcdfcpp.h>
30    #endif
31    
32  #if USE_SILO  #if USE_SILO
33  #include <silo.h>  #include <silo.h>
# Line 25  extern "C" { Line 38  extern "C" {
38    
39  #include <iomanip>  #include <iomanip>
40    
41    #include "esysUtils/EsysRandom.h"
42    #include "blocktools.h"
43    
44    
45  using namespace std;  using namespace std;
46    using esysUtils::FileWriter;
47    
48  namespace ripley {  namespace ripley {
49    
50  Brick::Brick(int n0, int n1, int n2, double l0, double l1, double l2, int d0,  int indexOfMax(int a, int b, int c) {
51               int d1, int d2) :      if (a > b) {
52      RipleyDomain(3),          if (c > a) {
53      m_gNE0(n0),              return 2;
54      m_gNE1(n1),          }
55      m_gNE2(n2),          return 0;
56      m_l0(l0),      } else if (b > c) {
57      m_l1(l1),          return 1;
58      m_l2(l2),      }
59      m_NX(d0),      return 2;
60      m_NY(d1),  }
61      m_NZ(d2)  
62    Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
63                 double x1, double y1, double z1, int d0, int d1, int d2,
64                 const std::vector<double>& points, const std::vector<int>& tags,
65                 const simap_t& tagnamestonums) :
66        RipleyDomain(3)
67  {  {
68        if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)
69                * static_cast<long>(n2 + 1) > std::numeric_limits<int>::max())
70            throw RipleyException("The number of elements has overflowed, this "
71                    "limit may be raised in future releases.");
72    
73        if (n0 <= 0 || n1 <= 0 || n2 <= 0)
74            throw RipleyException("Number of elements in each spatial dimension "
75                    "must be positive");
76    
77        // ignore subdivision parameters for serial run
78        if (m_mpiInfo->size == 1) {
79            d0=1;
80            d1=1;
81            d2=1;
82        }
83        bool warn=false;
84    
85        std::vector<int> factors;
86        int ranks = m_mpiInfo->size;
87        int epr[3] = {n0,n1,n2};
88        int d[3] = {d0,d1,d2};
89        if (d0<=0 || d1<=0 || d2<=0) {
90            for (int i = 0; i < 3; i++) {
91                if (d[i] < 1) {
92                    d[i] = 1;
93                    continue;
94                }
95                epr[i] = -1; // can no longer be max
96                if (ranks % d[i] != 0) {
97                    throw RipleyException("Invalid number of spatial subdivisions");
98                }
99                //remove
100                ranks /= d[i];
101            }
102            factorise(factors, ranks);
103            if (factors.size() != 0) {
104                warn = true;
105            }
106        }
107        while (factors.size() > 0) {
108            int i = indexOfMax(epr[0],epr[1],epr[2]);
109            int f = factors.back();
110            factors.pop_back();
111            d[i] *= f;
112            epr[i] /= f;
113        }
114        d0 = d[0]; d1 = d[1]; d2 = d[2];
115    
116      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
117      // among number of ranks      // among number of ranks
118      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (d0*d1*d2 != m_mpiInfo->size){
119          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
120        }
121        if (warn) {
122            cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
123                << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
124        }
125    
126        double l0 = x1-x0;
127        double l1 = y1-y0;
128        double l2 = z1-z0;
129        m_dx[0] = l0/n0;
130        m_dx[1] = l1/n1;
131        m_dx[2] = l2/n2;
132    
133        if ((n0+1)%d0 > 0) {
134            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
135            l0=m_dx[0]*n0;
136            cout << "Warning: Adjusted number of elements and length. N0="
137                << n0 << ", l0=" << l0 << endl;
138        }
139        if ((n1+1)%d1 > 0) {
140            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
141            l1=m_dx[1]*n1;
142            cout << "Warning: Adjusted number of elements and length. N1="
143                << n1 << ", l1=" << l1 << endl;
144        }
145        if ((n2+1)%d2 > 0) {
146            n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;
147            l2=m_dx[2]*n2;
148            cout << "Warning: Adjusted number of elements and length. N2="
149                << n2 << ", l2=" << l2 << endl;
150        }
151    
152        if ((d0 > 1 && (n0+1)/d0<2) || (d1 > 1 && (n1+1)/d1<2) || (d2 > 1 && (n2+1)/d2<2))
153            throw RipleyException("Too few elements for the number of ranks");
154    
155      if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)      m_gNE[0] = n0;
156          throw RipleyException("Number of elements must be separable into number of ranks in each dimension");      m_gNE[1] = n1;
157        m_gNE[2] = n2;
158        m_origin[0] = x0;
159        m_origin[1] = y0;
160        m_origin[2] = z0;
161        m_length[0] = l0;
162        m_length[1] = l1;
163        m_length[2] = l2;
164        m_NX[0] = d0;
165        m_NX[1] = d1;
166        m_NX[2] = d2;
167    
168        // local number of elements (including overlap)
169        m_NE[0] = m_ownNE[0] = (d0>1 ? (n0+1)/d0 : n0);
170        if (m_mpiInfo->rank%d0>0 && m_mpiInfo->rank%d0<d0-1)
171            m_NE[0]++;
172        else if (d0>1 && m_mpiInfo->rank%d0==d0-1)
173            m_ownNE[0]--;
174    
175        m_NE[1] = m_ownNE[1] = (d1>1 ? (n1+1)/d1 : n1);
176        if (m_mpiInfo->rank%(d0*d1)/d0>0 && m_mpiInfo->rank%(d0*d1)/d0<d1-1)
177            m_NE[1]++;
178        else if (d1>1 && m_mpiInfo->rank%(d0*d1)/d0==d1-1)
179            m_ownNE[1]--;
180    
181        m_NE[2] = m_ownNE[2] = (d2>1 ? (n2+1)/d2 : n2);
182        if (m_mpiInfo->rank/(d0*d1)>0 && m_mpiInfo->rank/(d0*d1)<d2-1)
183            m_NE[2]++;
184        else if (d2>1 && m_mpiInfo->rank/(d0*d1)==d2-1)
185            m_ownNE[2]--;
186    
187        // local number of nodes
188        m_NN[0] = m_NE[0]+1;
189        m_NN[1] = m_NE[1]+1;
190        m_NN[2] = m_NE[2]+1;
191    
     // local number of elements  
     m_NE0 = n0/m_NX;  
     m_NE1 = n1/m_NY;  
     m_NE2 = n2/m_NZ;  
     // local number of nodes (not necessarily owned)  
     m_N0 = m_NE0+1;  
     m_N1 = m_NE1+1;  
     m_N2 = m_NE2+1;  
192      // bottom-left-front node is at (offset0,offset1,offset2) in global mesh      // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
193      m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);      m_offset[0] = (n0+1)/d0*(m_mpiInfo->rank%d0);
194      m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);      if (m_offset[0] > 0)
195      m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));          m_offset[0]--;
196        m_offset[1] = (n1+1)/d1*(m_mpiInfo->rank%(d0*d1)/d0);
197        if (m_offset[1] > 0)
198            m_offset[1]--;
199        m_offset[2] = (n2+1)/d2*(m_mpiInfo->rank/(d0*d1));
200        if (m_offset[2] > 0)
201            m_offset[2]--;
202    
203      populateSampleIds();      populateSampleIds();
204        createPattern();
205        
206        assembler = new DefaultAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
207        for (map<string, int>::const_iterator i = tagnamestonums.begin();
208                i != tagnamestonums.end(); i++) {
209            setTagMap(i->first, i->second);
210        }
211        addPoints(tags.size(), &points[0], &tags[0]);
212  }  }
213    
214    
215  Brick::~Brick()  Brick::~Brick()
216  {  {
217        delete assembler;
218  }  }
219    
220  string Brick::getDescription() const  string Brick::getDescription() const
# Line 77  string Brick::getDescription() const Line 224  string Brick::getDescription() const
224    
225  bool Brick::operator==(const AbstractDomain& other) const  bool Brick::operator==(const AbstractDomain& other) const
226  {  {
227      if (dynamic_cast<const Brick*>(&other))      const Brick* o=dynamic_cast<const Brick*>(&other);
228          return this==&other;      if (o) {
229            return (RipleyDomain::operator==(other) &&
230                    m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1] && m_gNE[2]==o->m_gNE[2]
231                    && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1] && m_origin[2]==o->m_origin[2]
232                    && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1] && m_length[2]==o->m_length[2]
233                    && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1] && m_NX[2]==o->m_NX[2]);
234        }
235    
236      return false;      return false;
237  }  }
238    
239    void Brick::readNcGrid(escript::Data& out, string filename, string varname,
240                const ReaderParameters& params) const
241    {
242    #ifdef USE_NETCDF
243        // check destination function space
244        int myN0, myN1, myN2;
245        if (out.getFunctionSpace().getTypeCode() == Nodes) {
246            myN0 = m_NN[0];
247            myN1 = m_NN[1];
248            myN2 = m_NN[2];
249        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
250                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
251            myN0 = m_NE[0];
252            myN1 = m_NE[1];
253            myN2 = m_NE[2];
254        } else
255            throw RipleyException("readNcGrid(): invalid function space for output data object");
256    
257        if (params.first.size() != 3)
258            throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
259    
260        if (params.numValues.size() != 3)
261            throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
262    
263        if (params.multiplier.size() != 3)
264            throw RipleyException("readNcGrid(): argument 'multiplier' must have 3 entries");
265        for (size_t i=0; i<params.multiplier.size(); i++)
266            if (params.multiplier[i]<1)
267                throw RipleyException("readNcGrid(): all multipliers must be positive");
268    
269        // check file existence and size
270        NcFile f(filename.c_str(), NcFile::ReadOnly);
271        if (!f.is_valid())
272            throw RipleyException("readNcGrid(): cannot open file");
273    
274        NcVar* var = f.get_var(varname.c_str());
275        if (!var)
276            throw RipleyException("readNcGrid(): invalid variable name");
277    
278        // TODO: rank>0 data support
279        const int numComp = out.getDataPointSize();
280        if (numComp > 1)
281            throw RipleyException("readNcGrid(): only scalar data supported");
282    
283        const int dims = var->num_dims();
284        boost::scoped_array<long> edges(var->edges());
285    
286        // is this a slice of the data object (dims!=3)?
287        // note the expected ordering of edges (as in numpy: z,y,x)
288        if ( (dims==3 && (params.numValues[2] > edges[0] ||
289                          params.numValues[1] > edges[1] ||
290                          params.numValues[0] > edges[2]))
291                || (dims==2 && params.numValues[2]>1)
292                || (dims==1 && (params.numValues[2]>1 || params.numValues[1]>1)) ) {
293            throw RipleyException("readNcGrid(): not enough data in file");
294        }
295    
296        // check if this rank contributes anything
297        if (params.first[0] >= m_offset[0]+myN0 ||
298                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
299                params.first[1] >= m_offset[1]+myN1 ||
300                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
301                params.first[2] >= m_offset[2]+myN2 ||
302                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
303            return;
304        }
305    
306        // now determine how much this rank has to write
307    
308        // first coordinates in data object to write to
309        const int first0 = max(0, params.first[0]-m_offset[0]);
310        const int first1 = max(0, params.first[1]-m_offset[1]);
311        const int first2 = max(0, params.first[2]-m_offset[2]);
312        // indices to first value in file (not accounting for reverse yet)
313        int idx0 = max(0, m_offset[0]-params.first[0]);
314        int idx1 = max(0, m_offset[1]-params.first[1]);
315        int idx2 = max(0, m_offset[2]-params.first[2]);
316        // number of values to read
317        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
318        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
319        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
320    
321        // make sure we read the right block if going backwards through file
322        if (params.reverse[0])
323            idx0 = edges[dims-1]-num0-idx0;
324        if (dims>1 && params.reverse[1])
325            idx1 = edges[dims-2]-num1-idx1;
326        if (dims>2 && params.reverse[2])
327            idx2 = edges[dims-3]-num2-idx2;
328    
329    
330        vector<double> values(num0*num1*num2);
331        if (dims==3) {
332            var->set_cur(idx2, idx1, idx0);
333            var->get(&values[0], num2, num1, num0);
334        } else if (dims==2) {
335            var->set_cur(idx1, idx0);
336            var->get(&values[0], num1, num0);
337        } else {
338            var->set_cur(idx0);
339            var->get(&values[0], num0);
340        }
341    
342        const int dpp = out.getNumDataPointsPerSample();
343        out.requireWrite();
344    
345        // helpers for reversing
346        const int x0 = (params.reverse[0] ? num0-1 : 0);
347        const int x_mult = (params.reverse[0] ? -1 : 1);
348        const int y0 = (params.reverse[1] ? num1-1 : 0);
349        const int y_mult = (params.reverse[1] ? -1 : 1);
350        const int z0 = (params.reverse[2] ? num2-1 : 0);
351        const int z_mult = (params.reverse[2] ? -1 : 1);
352    
353        for (index_t z=0; z<num2; z++) {
354            for (index_t y=0; y<num1; y++) {
355    #pragma omp parallel for
356                for (index_t x=0; x<num0; x++) {
357                    const int baseIndex = first0+x*params.multiplier[0]
358                                         +(first1+y*params.multiplier[1])*myN0
359                                         +(first2+z*params.multiplier[2])*myN0*myN1;
360                    const int srcIndex=(z0+z_mult*z)*num1*num0
361                                      +(y0+y_mult*y)*num0
362                                      +(x0+x_mult*x);
363                    if (!isnan(values[srcIndex])) {
364                        for (index_t m2=0; m2<params.multiplier[2]; m2++) {
365                            for (index_t m1=0; m1<params.multiplier[1]; m1++) {
366                                for (index_t m0=0; m0<params.multiplier[0]; m0++) {
367                                    const int dataIndex = baseIndex+m0
368                                                   +m1*myN0
369                                                   +m2*myN0*myN1;
370                                    double* dest = out.getSampleDataRW(dataIndex);
371                                    for (index_t q=0; q<dpp; q++) {
372                                        *dest++ = values[srcIndex];
373                                    }
374                                }
375                            }
376                        }
377                    }
378                }
379            }
380        }
381    #else
382        throw RipleyException("readNcGrid(): not compiled with netCDF support");
383    #endif
384    }
385    
386    #ifdef USE_BOOSTIO
387    void Brick::readBinaryGridFromZipped(escript::Data& out, string filename,
388                               const ReaderParameters& params) const
389    {
390        // the mapping is not universally correct but should work on our
391        // supported platforms
392        switch (params.dataType) {
393            case DATATYPE_INT32:
394                readBinaryGridZippedImpl<int>(out, filename, params);
395                break;
396            case DATATYPE_FLOAT32:
397                readBinaryGridZippedImpl<float>(out, filename, params);
398                break;
399            case DATATYPE_FLOAT64:
400                readBinaryGridZippedImpl<double>(out, filename, params);
401                break;
402            default:
403                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
404        }
405    }
406    #endif
407    
408    void Brick::readBinaryGrid(escript::Data& out, string filename,
409                               const ReaderParameters& params) const
410    {
411        // the mapping is not universally correct but should work on our
412        // supported platforms
413        switch (params.dataType) {
414            case DATATYPE_INT32:
415                readBinaryGridImpl<int>(out, filename, params);
416                break;
417            case DATATYPE_FLOAT32:
418                readBinaryGridImpl<float>(out, filename, params);
419                break;
420            case DATATYPE_FLOAT64:
421                readBinaryGridImpl<double>(out, filename, params);
422                break;
423            default:
424                throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
425        }
426    }
427    
428    template<typename ValueType>
429    void Brick::readBinaryGridImpl(escript::Data& out, const string& filename,
430                                   const ReaderParameters& params) const
431    {
432        // check destination function space
433        int myN0, myN1, myN2;
434        if (out.getFunctionSpace().getTypeCode() == Nodes) {
435            myN0 = m_NN[0];
436            myN1 = m_NN[1];
437            myN2 = m_NN[2];
438        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
439                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
440            myN0 = m_NE[0];
441            myN1 = m_NE[1];
442            myN2 = m_NE[2];
443        } else
444            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
445    
446        if (params.first.size() != 3)
447            throw RipleyException("readBinaryGrid(): argument 'first' must have 3 entries");
448    
449        if (params.numValues.size() != 3)
450            throw RipleyException("readBinaryGrid(): argument 'numValues' must have 3 entries");
451    
452        if (params.multiplier.size() != 3)
453            throw RipleyException("readBinaryGrid(): argument 'multiplier' must have 3 entries");
454        for (size_t i=0; i<params.multiplier.size(); i++)
455            if (params.multiplier[i]<1)
456                throw RipleyException("readBinaryGrid(): all multipliers must be positive");
457    
458        // check file existence and size
459        ifstream f(filename.c_str(), ifstream::binary);
460        if (f.fail()) {
461            throw RipleyException("readBinaryGrid(): cannot open file");
462        }
463        f.seekg(0, ios::end);
464        const int numComp = out.getDataPointSize();
465        const int filesize = f.tellg();
466        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
467        if (filesize < reqsize) {
468            f.close();
469            throw RipleyException("readBinaryGrid(): not enough data in file");
470        }
471    
472        // check if this rank contributes anything
473        if (params.first[0] >= m_offset[0]+myN0 ||
474                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
475                params.first[1] >= m_offset[1]+myN1 ||
476                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
477                params.first[2] >= m_offset[2]+myN2 ||
478                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
479            f.close();
480            return;
481        }
482    
483        // now determine how much this rank has to write
484    
485        // first coordinates in data object to write to
486        const int first0 = max(0, params.first[0]-m_offset[0]);
487        const int first1 = max(0, params.first[1]-m_offset[1]);
488        const int first2 = max(0, params.first[2]-m_offset[2]);
489        // indices to first value in file
490        const int idx0 = max(0, m_offset[0]-params.first[0]);
491        const int idx1 = max(0, m_offset[1]-params.first[1]);
492        const int idx2 = max(0, m_offset[2]-params.first[2]);
493        // number of values to read
494        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
495        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
496        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
497    
498        out.requireWrite();
499        vector<ValueType> values(num0*numComp);
500        const int dpp = out.getNumDataPointsPerSample();
501    
502        for (int z=0; z<num2; z++) {
503            for (int y=0; y<num1; y++) {
504                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
505                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
506                f.seekg(fileofs*sizeof(ValueType));
507                f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
508    
509                for (int x=0; x<num0; x++) {
510                    const int baseIndex = first0+x*params.multiplier[0]
511                                         +(first1+y*params.multiplier[1])*myN0
512                                         +(first2+z*params.multiplier[2])*myN0*myN1;
513                    for (int m2=0; m2<params.multiplier[2]; m2++) {
514                        for (int m1=0; m1<params.multiplier[1]; m1++) {
515                            for (int m0=0; m0<params.multiplier[0]; m0++) {
516                                const int dataIndex = baseIndex+m0
517                                               +m1*myN0
518                                               +m2*myN0*myN1;
519                                double* dest = out.getSampleDataRW(dataIndex);
520                                for (int c=0; c<numComp; c++) {
521                                    ValueType val = values[x*numComp+c];
522    
523                                    if (params.byteOrder != BYTEORDER_NATIVE) {
524                                        char* cval = reinterpret_cast<char*>(&val);
525                                        // this will alter val!!
526                                        byte_swap32(cval);
527                                    }
528                                    if (!std::isnan(val)) {
529                                        for (int q=0; q<dpp; q++) {
530                                            *dest++ = static_cast<double>(val);
531                                        }
532                                    }
533                                }
534                            }
535                        }
536                    }
537                }
538            }
539        }
540    
541        f.close();
542    }
543    
544    #ifdef USE_BOOSTIO
545    template<typename ValueType>
546    void Brick::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
547                                   const ReaderParameters& params) const
548    {
549        // check destination function space
550        int myN0, myN1, myN2;
551        if (out.getFunctionSpace().getTypeCode() == Nodes) {
552            myN0 = m_NN[0];
553            myN1 = m_NN[1];
554            myN2 = m_NN[2];
555        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
556                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
557            myN0 = m_NE[0];
558            myN1 = m_NE[1];
559            myN2 = m_NE[2];
560        } else
561            throw RipleyException("readBinaryGridFromZipped(): invalid function space for output data object");
562    
563        if (params.first.size() != 3)
564            throw RipleyException("readBinaryGridFromZipped(): argument 'first' must have 3 entries");
565    
566        if (params.numValues.size() != 3)
567            throw RipleyException("readBinaryGridFromZipped(): argument 'numValues' must have 3 entries");
568    
569        if (params.multiplier.size() != 3)
570            throw RipleyException("readBinaryGridFromZipped(): argument 'multiplier' must have 3 entries");
571        for (size_t i=0; i<params.multiplier.size(); i++)
572            if (params.multiplier[i]<1)
573                throw RipleyException("readBinaryGridFromZipped(): all multipliers must be positive");
574    
575        // check file existence and size
576        ifstream f(filename.c_str(), ifstream::binary);
577        if (f.fail()) {
578            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
579        }
580        f.seekg(0, ios::end);
581        const int numComp = out.getDataPointSize();
582        int filesize = f.tellg();
583        f.seekg(0, ios::beg);
584        std::vector<char> compressed(filesize);
585        f.read((char*)&compressed[0], filesize);
586        f.close();
587        std::vector<char> decompressed = unzip(compressed);
588        filesize = decompressed.size();
589        const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
590        if (filesize < reqsize) {
591            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
592        }
593    
594        // check if this rank contributes anything
595        if (params.first[0] >= m_offset[0]+myN0 ||
596                params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
597                params.first[1] >= m_offset[1]+myN1 ||
598                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
599                params.first[2] >= m_offset[2]+myN2 ||
600                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
601            return;
602        }
603    
604        // now determine how much this rank has to write
605    
606        // first coordinates in data object to write to
607        const int first0 = max(0, params.first[0]-m_offset[0]);
608        const int first1 = max(0, params.first[1]-m_offset[1]);
609        const int first2 = max(0, params.first[2]-m_offset[2]);
610        // indices to first value in file
611        const int idx0 = max(0, m_offset[0]-params.first[0]);
612        const int idx1 = max(0, m_offset[1]-params.first[1]);
613        const int idx2 = max(0, m_offset[2]-params.first[2]);
614        // number of values to read
615        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
616        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
617        const int num2 = min(params.numValues[2]-idx2, myN2-first2);
618    
619        out.requireWrite();
620        vector<ValueType> values(num0*numComp);
621        const int dpp = out.getNumDataPointsPerSample();
622    
623        for (int z=0; z<num2; z++) {
624            for (int y=0; y<num1; y++) {
625                const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
626                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
627                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
628                
629                for (int x=0; x<num0; x++) {
630                    const int baseIndex = first0+x*params.multiplier[0]
631                                         +(first1+y*params.multiplier[1])*myN0
632                                         +(first2+z*params.multiplier[2])*myN0*myN1;
633                    for (int m2=0; m2<params.multiplier[2]; m2++) {
634                        for (int m1=0; m1<params.multiplier[1]; m1++) {
635                            for (int m0=0; m0<params.multiplier[0]; m0++) {
636                                const int dataIndex = baseIndex+m0
637                                               +m1*myN0
638                                               +m2*myN0*myN1;
639                                double* dest = out.getSampleDataRW(dataIndex);
640                                for (int c=0; c<numComp; c++) {
641                                    ValueType val = values[x*numComp+c];
642    
643                                    if (params.byteOrder != BYTEORDER_NATIVE) {
644                                        char* cval = reinterpret_cast<char*>(&val);
645                                        // this will alter val!!
646                                        byte_swap32(cval);
647                                    }
648                                    if (!std::isnan(val)) {
649                                        for (int q=0; q<dpp; q++) {
650                                            *dest++ = static_cast<double>(val);
651                                        }
652                                    }
653                                }
654                            }
655                        }
656                    }
657                }
658            }
659        }
660    }
661    #endif
662    
663    void Brick::writeBinaryGrid(const escript::Data& in, string filename,
664                                int byteOrder, int dataType) const
665    {
666        // the mapping is not universally correct but should work on our
667        // supported platforms
668        switch (dataType) {
669            case DATATYPE_INT32:
670                writeBinaryGridImpl<int>(in, filename, byteOrder);
671                break;
672            case DATATYPE_FLOAT32:
673                writeBinaryGridImpl<float>(in, filename, byteOrder);
674                break;
675            case DATATYPE_FLOAT64:
676                writeBinaryGridImpl<double>(in, filename, byteOrder);
677                break;
678            default:
679                throw RipleyException("writeBinaryGrid(): invalid or unsupported datatype");
680        }
681    }
682    
683    template<typename ValueType>
684    void Brick::writeBinaryGridImpl(const escript::Data& in,
685                                    const string& filename, int byteOrder) const
686    {
687        // check function space and determine number of points
688        int myN0, myN1, myN2;
689        int totalN0, totalN1, totalN2;
690        if (in.getFunctionSpace().getTypeCode() == Nodes) {
691            myN0 = m_NN[0];
692            myN1 = m_NN[1];
693            myN2 = m_NN[2];
694            totalN0 = m_gNE[0]+1;
695            totalN1 = m_gNE[1]+1;
696            totalN2 = m_gNE[2]+1;
697        } else if (in.getFunctionSpace().getTypeCode() == Elements ||
698                    in.getFunctionSpace().getTypeCode() == ReducedElements) {
699            myN0 = m_NE[0];
700            myN1 = m_NE[1];
701            myN2 = m_NE[2];
702            totalN0 = m_gNE[0];
703            totalN1 = m_gNE[1];
704            totalN2 = m_gNE[2];
705        } else
706            throw RipleyException("writeBinaryGrid(): invalid function space of data object");
707    
708        const int numComp = in.getDataPointSize();
709        const int dpp = in.getNumDataPointsPerSample();
710        const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1*totalN2;
711    
712        if (numComp > 1 || dpp > 1)
713            throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
714    
715        // from here on we know that each sample consists of one value
716        FileWriter fw;
717        fw.openFile(filename, fileSize);
718        MPIBarrier();
719    
720        for (index_t z=0; z<myN2; z++) {
721            for (index_t y=0; y<myN1; y++) {
722                const int fileofs = (m_offset[0]+(m_offset[1]+y)*totalN0
723                                    +(m_offset[2]+z)*totalN0*totalN1)*sizeof(ValueType);
724                ostringstream oss;
725    
726                for (index_t x=0; x<myN0; x++) {
727                    const double* sample = in.getSampleDataRO(z*myN0*myN1+y*myN0+x);
728                    ValueType fvalue = static_cast<ValueType>(*sample);
729                    if (byteOrder == BYTEORDER_NATIVE) {
730                        oss.write((char*)&fvalue, sizeof(fvalue));
731                    } else {
732                        char* value = reinterpret_cast<char*>(&fvalue);
733                        oss.write(byte_swap32(value), sizeof(fvalue));
734                    }
735                }
736                fw.writeAt(oss, fileofs);
737            }
738        }
739        fw.close();
740    }
741    
742  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
743  {  {
744  #if USE_SILO  #if USE_SILO
# Line 91  void Brick::dump(const string& fileName) Line 747  void Brick::dump(const string& fileName)
747          fn+=".silo";          fn+=".silo";
748      }      }
749    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
750      int driver=DB_HDF5;          int driver=DB_HDF5;    
751      string siloPath;      string siloPath;
752      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
753    
754  #ifdef ESYS_MPI  #ifdef ESYS_MPI
755      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
756        const int NUM_SILO_FILES = 1;
757        const char* blockDirFmt = "/block%04d";
758  #endif  #endif
759    
760      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 143  void Brick::dump(const string& fileName) Line 799  void Brick::dump(const string& fileName)
799      }      }
800      */      */
801    
802      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_NN[0]]);
803      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
804      boost::scoped_ptr<double> z(new double[m_N2]);      boost::scoped_ptr<double> z(new double[m_NN[2]]);
805      double* coords[3] = { x.get(), y.get(), z.get() };      double* coords[3] = { x.get(), y.get(), z.get() };
     pair<double,double> xdx = getFirstCoordAndSpacing(0);  
     pair<double,double> ydy = getFirstCoordAndSpacing(1);  
     pair<double,double> zdz = getFirstCoordAndSpacing(2);  
806  #pragma omp parallel  #pragma omp parallel
807      {      {
808  #pragma omp for  #pragma omp for
809          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
810              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=getLocalCoordinate(i0, 0);
811          }          }
812  #pragma omp for  #pragma omp for
813          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
814              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=getLocalCoordinate(i1, 1);
815          }          }
816  #pragma omp for  #pragma omp for
817          for (dim_t i2 = 0; i2 < m_N2; i2++) {          for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {
818              coords[2][i2]=zdz.first+i2*zdz.second;              coords[2][i2]=getLocalCoordinate(i2, 2);
819          }          }
820      }      }
821      IndexVector dims = getNumNodesPerDim();      int* dims = const_cast<int*>(getNumNodesPerDim());
822      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,  
823        // write mesh
824        DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,
825              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
826    
827      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,      // write node ids
828        DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 3,
829              NULL, 0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
830    
831      // write element ids      // write element ids
832      dims = getNumElementsPerDim();      dims = const_cast<int*>(getNumElementsPerDim());
833      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],      DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
834              &dims[0], 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);              dims, 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
835    
836      // rank 0 writes multimesh and multivar      // rank 0 writes multimesh and multivar
837      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
# Line 224  void Brick::dump(const string& fileName) Line 880  void Brick::dump(const string& fileName)
880      }      }
881    
882  #else // USE_SILO  #else // USE_SILO
883      throw RipleyException("dump(): no Silo support");      throw RipleyException("dump: no Silo support");
884  #endif  #endif
885  }  }
886    
# Line 232  const int* Brick::borrowSampleReferenceI Line 888  const int* Brick::borrowSampleReferenceI
888  {  {
889      switch (fsType) {      switch (fsType) {
890          case Nodes:          case Nodes:
891            case ReducedNodes: //FIXME: reduced
892              return &m_nodeId[0];              return &m_nodeId[0];
893            case DegreesOfFreedom:
894            case ReducedDegreesOfFreedom: //FIXME: reduced
895                return &m_dofId[0];
896          case Elements:          case Elements:
897            case ReducedElements:
898              return &m_elementId[0];              return &m_elementId[0];
899          case FaceElements:          case FaceElements:
900            case ReducedFaceElements:
901              return &m_faceId[0];              return &m_faceId[0];
902            case Points:
903                return &m_diracPointNodeIDs[0];
904          default:          default:
905              break;              break;
906      }      }
907    
908      stringstream msg;      stringstream msg;
909      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs: invalid function space type "<<fsType;
         << fsType;  
910      throw RipleyException(msg.str());      throw RipleyException(msg.str());
911  }  }
912    
913  bool Brick::ownSample(int fsCode, index_t id) const  bool Brick::ownSample(int fsType, index_t id) const
914  {  {
915  #ifdef ESYS_MPI      if (getMPISize()==1)
916      if (fsCode == Nodes) {          return true;
917          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];  
918          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;      switch (fsType) {
919          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          case Nodes:
920      } else          case ReducedNodes: //FIXME: reduced
921          throw RipleyException("ownSample() only implemented for Nodes");              return (m_dofMap[id] < getNumDOF());
922  #else          case DegreesOfFreedom:
923      return true;          case ReducedDegreesOfFreedom:
924  #endif              return true;
925            case Elements:
926            case ReducedElements:
927                {
928                    // check ownership of element's _last_ node
929                    const index_t x=id%m_NE[0] + 1;
930                    const index_t y=id%(m_NE[0]*m_NE[1])/m_NE[0] + 1;
931                    const index_t z=id/(m_NE[0]*m_NE[1]) + 1;
932                    return (m_dofMap[x + m_NN[0]*y + m_NN[0]*m_NN[1]*z] < getNumDOF());
933                }
934            case FaceElements:
935            case ReducedFaceElements:
936                {
937                    // check ownership of face element's last node
938                    dim_t n=0;
939                    for (size_t i=0; i<6; i++) {
940                        n+=m_faceCount[i];
941                        if (id<n) {
942                            const index_t j=id-n+m_faceCount[i];
943                            if (i>=4) { // front or back
944                                const index_t first=(i==4 ? 0 : m_NN[0]*m_NN[1]*(m_NN[2]-1));
945                                return (m_dofMap[first+j%m_NE[0]+1+(j/m_NE[0]+1)*m_NN[0]] < getNumDOF());
946                            } else if (i>=2) { // bottom or top
947                                const index_t first=(i==2 ? 0 : m_NN[0]*(m_NN[1]-1));
948                                return (m_dofMap[first+j%m_NE[0]+1+(j/m_NE[0]+1)*m_NN[0]*m_NN[1]] < getNumDOF());
949                            } else { // left or right
950                                const index_t first=(i==0 ? 0 : m_NN[0]-1);
951                                return (m_dofMap[first+(j%m_NE[1]+1)*m_NN[0]+(j/m_NE[1]+1)*m_NN[0]*m_NN[1]] < getNumDOF());
952                            }
953                        }
954                    }
955                    return false;
956                }
957            default:
958                break;
959        }
960    
961        stringstream msg;
962        msg << "ownSample: invalid function space type " << fsType;
963        throw RipleyException(msg.str());
964  }  }
965    
966  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  void Brick::setToNormal(escript::Data& out) const
                                             bool reducedColOrder) const  
967  {  {
968      if (reducedRowOrder || reducedColOrder)      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
969          throw RipleyException("getPattern() not implemented for reduced order");          out.requireWrite();
970    #pragma omp parallel
971            {
972                if (m_faceOffset[0] > -1) {
973    #pragma omp for nowait
974                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
975                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
976                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
977                            // set vector at four quadrature points
978                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
979                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
980                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
981                            *o++ = -1.; *o++ = 0.; *o = 0.;
982                        }
983                    }
984                }
985    
986                if (m_faceOffset[1] > -1) {
987    #pragma omp for nowait
988                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
989                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
990                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
991                            // set vector at four quadrature points
992                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
993                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
994                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
995                            *o++ = 1.; *o++ = 0.; *o = 0.;
996                        }
997                    }
998                }
999    
1000                if (m_faceOffset[2] > -1) {
1001    #pragma omp for nowait
1002                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1003                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1004                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1005                            // set vector at four quadrature points
1006                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1007                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1008                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1009                            *o++ = 0.; *o++ = -1.; *o = 0.;
1010                        }
1011                    }
1012                }
1013    
1014                if (m_faceOffset[3] > -1) {
1015    #pragma omp for nowait
1016                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1017                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1018                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1019                            // set vector at four quadrature points
1020                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1021                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1022                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1023                            *o++ = 0.; *o++ = 1.; *o = 0.;
1024                        }
1025                    }
1026                }
1027    
1028                if (m_faceOffset[4] > -1) {
1029    #pragma omp for nowait
1030                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1031                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1032                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1033                            // set vector at four quadrature points
1034                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1035                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1036                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1037                            *o++ = 0.; *o++ = 0.; *o = -1.;
1038                        }
1039                    }
1040                }
1041    
1042                if (m_faceOffset[5] > -1) {
1043    #pragma omp for nowait
1044                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1045                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1046                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1047                            // set vector at four quadrature points
1048                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1049                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1050                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1051                            *o++ = 0.; *o++ = 0.; *o = 1.;
1052                        }
1053                    }
1054                }
1055            } // end of parallel section
1056        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1057            out.requireWrite();
1058    #pragma omp parallel
1059            {
1060                if (m_faceOffset[0] > -1) {
1061    #pragma omp for nowait
1062                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1063                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1064                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1065                            *o++ = -1.;
1066                            *o++ = 0.;
1067                            *o = 0.;
1068                        }
1069                    }
1070                }
1071    
1072                if (m_faceOffset[1] > -1) {
1073    #pragma omp for nowait
1074                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1075                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1076                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1077                            *o++ = 1.;
1078                            *o++ = 0.;
1079                            *o = 0.;
1080                        }
1081                    }
1082                }
1083    
1084                if (m_faceOffset[2] > -1) {
1085    #pragma omp for nowait
1086                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1087                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1088                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1089                            *o++ = 0.;
1090                            *o++ = -1.;
1091                            *o = 0.;
1092                        }
1093                    }
1094                }
1095    
1096                if (m_faceOffset[3] > -1) {
1097    #pragma omp for nowait
1098                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1099                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1100                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1101                            *o++ = 0.;
1102                            *o++ = 1.;
1103                            *o = 0.;
1104                        }
1105                    }
1106                }
1107    
1108                if (m_faceOffset[4] > -1) {
1109    #pragma omp for nowait
1110                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1111                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1112                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1113                            *o++ = 0.;
1114                            *o++ = 0.;
1115                            *o = -1.;
1116                        }
1117                    }
1118                }
1119    
1120      throw RipleyException("getPattern() not implemented");              if (m_faceOffset[5] > -1) {
1121    #pragma omp for nowait
1122                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1123                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1124                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1125                            *o++ = 0.;
1126                            *o++ = 0.;
1127                            *o = 1.;
1128                        }
1129                    }
1130                }
1131            } // end of parallel section
1132    
1133        } else {
1134            stringstream msg;
1135            msg << "setToNormal: invalid function space type "
1136                << out.getFunctionSpace().getTypeCode();
1137            throw RipleyException(msg.str());
1138        }
1139    }
1140    
1141    void Brick::setToSize(escript::Data& out) const
1142    {
1143        if (out.getFunctionSpace().getTypeCode() == Elements
1144                || out.getFunctionSpace().getTypeCode() == ReducedElements) {
1145            out.requireWrite();
1146            const dim_t numQuad=out.getNumDataPointsPerSample();
1147            const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);
1148    #pragma omp parallel for
1149            for (index_t k = 0; k < getNumElements(); ++k) {
1150                double* o = out.getSampleDataRW(k);
1151                fill(o, o+numQuad, size);
1152            }
1153        } else if (out.getFunctionSpace().getTypeCode() == FaceElements
1154                || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1155            out.requireWrite();
1156            const dim_t numQuad=out.getNumDataPointsPerSample();
1157    #pragma omp parallel
1158            {
1159                if (m_faceOffset[0] > -1) {
1160                    const double size=min(m_dx[1],m_dx[2]);
1161    #pragma omp for nowait
1162                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1163                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1164                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1165                            fill(o, o+numQuad, size);
1166                        }
1167                    }
1168                }
1169    
1170                if (m_faceOffset[1] > -1) {
1171                    const double size=min(m_dx[1],m_dx[2]);
1172    #pragma omp for nowait
1173                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1174                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1175                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1176                            fill(o, o+numQuad, size);
1177                        }
1178                    }
1179                }
1180    
1181                if (m_faceOffset[2] > -1) {
1182                    const double size=min(m_dx[0],m_dx[2]);
1183    #pragma omp for nowait
1184                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1185                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1186                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1187                            fill(o, o+numQuad, size);
1188                        }
1189                    }
1190                }
1191    
1192                if (m_faceOffset[3] > -1) {
1193                    const double size=min(m_dx[0],m_dx[2]);
1194    #pragma omp for nowait
1195                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
1196                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1197                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1198                            fill(o, o+numQuad, size);
1199                        }
1200                    }
1201                }
1202    
1203                if (m_faceOffset[4] > -1) {
1204                    const double size=min(m_dx[0],m_dx[1]);
1205    #pragma omp for nowait
1206                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1207                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1208                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1209                            fill(o, o+numQuad, size);
1210                        }
1211                    }
1212                }
1213    
1214                if (m_faceOffset[5] > -1) {
1215                    const double size=min(m_dx[0],m_dx[1]);
1216    #pragma omp for nowait
1217                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
1218                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
1219                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1220                            fill(o, o+numQuad, size);
1221                        }
1222                    }
1223                }
1224            } // end of parallel section
1225    
1226        } else {
1227            stringstream msg;
1228            msg << "setToSize: invalid function space type "
1229                << out.getFunctionSpace().getTypeCode();
1230            throw RipleyException(msg.str());
1231        }
1232  }  }
1233    
1234  void Brick::Print_Mesh_Info(const bool full) const  void Brick::Print_Mesh_Info(const bool full) const
# Line 277  void Brick::Print_Mesh_Info(const bool f Line 1238  void Brick::Print_Mesh_Info(const bool f
1238          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
1239          cout.precision(15);          cout.precision(15);
1240          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
         pair<double,double> xdx = getFirstCoordAndSpacing(0);  
         pair<double,double> ydy = getFirstCoordAndSpacing(1);  
         pair<double,double> zdz = getFirstCoordAndSpacing(2);  
1241          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
1242              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
1243                  << "  " << xdx.first+(i%m_N0)*xdx.second                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
1244                  << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second                  << "  " << getLocalCoordinate(i%(m_NN[0]*m_NN[1])/m_NN[0], 1)
1245                  << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;                  << "  " << getLocalCoordinate(i/(m_NN[0]*m_NN[1]), 2) << endl;
1246          }          }
1247      }      }
1248  }  }
1249    
 IndexVector Brick::getNumNodesPerDim() const  
 {  
     IndexVector ret;  
     ret.push_back(m_N0);  
     ret.push_back(m_N1);  
     ret.push_back(m_N2);  
     return ret;  
 }  
1250    
1251  IndexVector Brick::getNumElementsPerDim() const  //protected
1252    void Brick::assembleCoordinates(escript::Data& arg) const
1253  {  {
1254      IndexVector ret;      escriptDataC x = arg.getDataC();
1255      ret.push_back(m_NE0);      int numDim = m_numDim;
1256      ret.push_back(m_NE1);      if (!isDataPointShapeEqual(&x, 1, &numDim))
1257      ret.push_back(m_NE2);          throw RipleyException("setToX: Invalid Data object shape");
1258      return ret;      if (!numSamplesEqual(&x, 1, getNumNodes()))
1259            throw RipleyException("setToX: Illegal number of samples in Data object");
1260    
1261        arg.requireWrite();
1262    #pragma omp parallel for
1263        for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {
1264            for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
1265                for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
1266                    double* point = arg.getSampleDataRW(i0+m_NN[0]*i1+m_NN[0]*m_NN[1]*i2);
1267                    point[0] = getLocalCoordinate(i0, 0);
1268                    point[1] = getLocalCoordinate(i1, 1);
1269                    point[2] = getLocalCoordinate(i2, 2);
1270                }
1271            }
1272        }
1273  }  }
1274    
1275  IndexVector Brick::getNumFacesPerBoundary() const  //protected
1276    void Brick::assembleGradient(escript::Data& out, const escript::Data& in) const
1277  {  {
1278      IndexVector ret(6, 0);      const dim_t numComp = in.getDataPointSize();
1279      //left      const double C0 = .044658198738520451079;
1280      if (m_offset0==0)      const double C1 = .16666666666666666667;
1281          ret[0]=m_NE1*m_NE2;      const double C2 = .21132486540518711775;
1282      //right      const double C3 = .25;
1283      if (m_mpiInfo->rank%m_NX==m_NX-1)      const double C4 = .5;
1284          ret[1]=m_NE1*m_NE2;      const double C5 = .62200846792814621559;
1285      //bottom      const double C6 = .78867513459481288225;
1286      if (m_offset1==0)  
1287          ret[2]=m_NE0*m_NE2;      if (out.getFunctionSpace().getTypeCode() == Elements) {
1288      //top          out.requireWrite();
1289      if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)  #pragma omp parallel
1290          ret[3]=m_NE0*m_NE2;          {
1291      //front              vector<double> f_000(numComp);
1292      if (m_offset2==0)              vector<double> f_001(numComp);
1293          ret[4]=m_NE0*m_NE1;              vector<double> f_010(numComp);
1294      //back              vector<double> f_011(numComp);
1295      if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)              vector<double> f_100(numComp);
1296          ret[5]=m_NE0*m_NE1;              vector<double> f_101(numComp);
1297      return ret;              vector<double> f_110(numComp);
1298                vector<double> f_111(numComp);
1299    #pragma omp for
1300                for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1301                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1302                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1303                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1304                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1305                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1306                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1307                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1308                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1309                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1310                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1311                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
1312                            for (index_t i=0; i < numComp; ++i) {
1313                                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];
1314                                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];
1315                                const double V2=((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];
1316                                const double V3=((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];
1317                                const double V4=((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1318                                const double V5=((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1319                                const double V6=((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1320                                const double V7=((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1321                                const double V8=((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1322                                const double V9=((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1323                                const double V10=((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1324                                const double V11=((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1325                                o[INDEX3(i,0,0,numComp,3)] = V0;
1326                                o[INDEX3(i,1,0,numComp,3)] = V4;
1327                                o[INDEX3(i,2,0,numComp,3)] = V8;
1328                                o[INDEX3(i,0,1,numComp,3)] = V0;
1329                                o[INDEX3(i,1,1,numComp,3)] = V5;
1330                                o[INDEX3(i,2,1,numComp,3)] = V9;
1331                                o[INDEX3(i,0,2,numComp,3)] = V1;
1332                                o[INDEX3(i,1,2,numComp,3)] = V4;
1333                                o[INDEX3(i,2,2,numComp,3)] = V10;
1334                                o[INDEX3(i,0,3,numComp,3)] = V1;
1335                                o[INDEX3(i,1,3,numComp,3)] = V5;
1336                                o[INDEX3(i,2,3,numComp,3)] = V11;
1337                                o[INDEX3(i,0,4,numComp,3)] = V2;
1338                                o[INDEX3(i,1,4,numComp,3)] = V6;
1339                                o[INDEX3(i,2,4,numComp,3)] = V8;
1340                                o[INDEX3(i,0,5,numComp,3)] = V2;
1341                                o[INDEX3(i,1,5,numComp,3)] = V7;
1342                                o[INDEX3(i,2,5,numComp,3)] = V9;
1343                                o[INDEX3(i,0,6,numComp,3)] = V3;
1344                                o[INDEX3(i,1,6,numComp,3)] = V6;
1345                                o[INDEX3(i,2,6,numComp,3)] = V10;
1346                                o[INDEX3(i,0,7,numComp,3)] = V3;
1347                                o[INDEX3(i,1,7,numComp,3)] = V7;
1348                                o[INDEX3(i,2,7,numComp,3)] = V11;
1349                            } // end of component loop i
1350                        } // end of k0 loop
1351                    } // end of k1 loop
1352                } // end of k2 loop
1353            } // end of parallel section
1354        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
1355            out.requireWrite();
1356    #pragma omp parallel
1357            {
1358                vector<double> f_000(numComp);
1359                vector<double> f_001(numComp);
1360                vector<double> f_010(numComp);
1361                vector<double> f_011(numComp);
1362                vector<double> f_100(numComp);
1363                vector<double> f_101(numComp);
1364                vector<double> f_110(numComp);
1365                vector<double> f_111(numComp);
1366    #pragma omp for
1367                for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1368                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1369                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1370                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1371                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1372                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1373                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1374                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1375                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1376                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1377                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1378                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
1379                            for (index_t i=0; i < numComp; ++i) {
1380                                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];
1381                                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];
1382                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / m_dx[2];
1383                            } // end of component loop i
1384                        } // end of k0 loop
1385                    } // end of k1 loop
1386                } // end of k2 loop
1387            } // end of parallel section
1388        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1389            out.requireWrite();
1390    #pragma omp parallel
1391            {
1392                vector<double> f_000(numComp);
1393                vector<double> f_001(numComp);
1394                vector<double> f_010(numComp);
1395                vector<double> f_011(numComp);
1396                vector<double> f_100(numComp);
1397                vector<double> f_101(numComp);
1398                vector<double> f_110(numComp);
1399                vector<double> f_111(numComp);
1400                if (m_faceOffset[0] > -1) {
1401    #pragma omp for nowait
1402                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1403                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1404                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1405                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1406                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1407                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1408                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1409                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1410                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1411                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1412                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1413                            for (index_t i=0; i < numComp; ++i) {
1414                                const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];
1415                                const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];
1416                                const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / m_dx[2];
1417                                const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / m_dx[2];
1418                                o[INDEX3(i,0,0,numComp,3)] = ((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];
1419                                o[INDEX3(i,1,0,numComp,3)] = V0;
1420                                o[INDEX3(i,2,0,numComp,3)] = V2;
1421                                o[INDEX3(i,0,1,numComp,3)] = ((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];
1422                                o[INDEX3(i,1,1,numComp,3)] = V0;
1423                                o[INDEX3(i,2,1,numComp,3)] = V3;
1424                                o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];
1425                                o[INDEX3(i,1,2,numComp,3)] = V1;
1426                                o[INDEX3(i,2,2,numComp,3)] = V2;
1427                                o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];
1428                                o[INDEX3(i,1,3,numComp,3)] = V1;
1429                                o[INDEX3(i,2,3,numComp,3)] = V3;
1430                            } // end of component loop i
1431                        } // end of k1 loop
1432                    } // end of k2 loop
1433                } // end of face 0
1434                if (m_faceOffset[1] > -1) {
1435    #pragma omp for nowait
1436                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1437                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1438                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1439                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1440                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1441                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1442                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1443                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1444                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1445                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1446                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1447                            for (index_t i=0; i < numComp; ++i) {
1448                                const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1449                                const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
1450                                const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / m_dx[2];
1451                                const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / m_dx[2];
1452                                o[INDEX3(i,0,0,numComp,3)] = ((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];
1453                                o[INDEX3(i,1,0,numComp,3)] = V0;
1454                                o[INDEX3(i,2,0,numComp,3)] = V2;
1455                                o[INDEX3(i,0,1,numComp,3)] = ((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];
1456                                o[INDEX3(i,1,1,numComp,3)] = V0;
1457                                o[INDEX3(i,2,1,numComp,3)] = V3;
1458                                o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];
1459                                o[INDEX3(i,1,2,numComp,3)] = V1;
1460                                o[INDEX3(i,2,2,numComp,3)] = V2;
1461                                o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];
1462                                o[INDEX3(i,1,3,numComp,3)] = V1;
1463                                o[INDEX3(i,2,3,numComp,3)] = V3;
1464                            } // end of component loop i
1465                        } // end of k1 loop
1466                    } // end of k2 loop
1467                } // end of face 1
1468                if (m_faceOffset[2] > -1) {
1469    #pragma omp for nowait
1470                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1471                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1472                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1473                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1474                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1475                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1476                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1477                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1478                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1479                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1480                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1481                            for (index_t i=0; i < numComp; ++i) {
1482                                const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];
1483                                const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];
1484                                const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / m_dx[2];
1485                                o[INDEX3(i,0,0,numComp,3)] = V0;
1486                                o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1487                                o[INDEX3(i,2,0,numComp,3)] = V1;
1488                                o[INDEX3(i,0,1,numComp,3)] = V0;
1489                                o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1490                                o[INDEX3(i,2,1,numComp,3)] = V2;
1491                                o[INDEX3(i,0,2,numComp,3)] = V0;
1492                                o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1493                                o[INDEX3(i,2,2,numComp,3)] = V1;
1494                                o[INDEX3(i,0,3,numComp,3)] = V0;
1495                                o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1496                                o[INDEX3(i,2,3,numComp,3)] = V2;
1497                            } // end of component loop i
1498                        } // end of k0 loop
1499                    } // end of k2 loop
1500                } // end of face 2
1501                if (m_faceOffset[3] > -1) {
1502    #pragma omp for nowait
1503                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1504                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1505                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1506                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1507                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1508                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1509                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1510                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1511                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1512                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1513                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1514                            for (index_t i=0; i < numComp; ++i) {
1515                                const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1516                                const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
1517                                const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / m_dx[2];
1518                                const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / m_dx[2];
1519                                o[INDEX3(i,0,0,numComp,3)] = V0;
1520                                o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1521                                o[INDEX3(i,2,0,numComp,3)] = V2;
1522                                o[INDEX3(i,0,1,numComp,3)] = V0;
1523                                o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1524                                o[INDEX3(i,2,1,numComp,3)] = V3;
1525                                o[INDEX3(i,0,2,numComp,3)] = V1;
1526                                o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / m_dx[1];
1527                                o[INDEX3(i,2,2,numComp,3)] = V2;
1528                                o[INDEX3(i,0,3,numComp,3)] = V1;
1529                                o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / m_dx[1];
1530                                o[INDEX3(i,2,3,numComp,3)] = V3;
1531                            } // end of component loop i
1532                        } // end of k0 loop
1533                    } // end of k2 loop
1534                } // end of face 3
1535                if (m_faceOffset[4] > -1) {
1536    #pragma omp for nowait
1537                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1538                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1539                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1540                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1541                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1542                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1543                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1544                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1545                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1546                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1547                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1548                            for (index_t i=0; i < numComp; ++i) {
1549                                const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];
1550                                const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];
1551                                const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / m_dx[1];
1552                                const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / m_dx[1];
1553                                o[INDEX3(i,0,0,numComp,3)] = V0;
1554                                o[INDEX3(i,1,0,numComp,3)] = V2;
1555                                o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1556                                o[INDEX3(i,0,1,numComp,3)] = V0;
1557                                o[INDEX3(i,1,1,numComp,3)] = V3;
1558                                o[INDEX3(i,2,1,numComp,3)] = ((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1559                                o[INDEX3(i,0,2,numComp,3)] = V1;
1560                                o[INDEX3(i,1,2,numComp,3)] = V2;
1561                                o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1562                                o[INDEX3(i,0,3,numComp,3)] = V1;
1563                                o[INDEX3(i,1,3,numComp,3)] = V3;
1564                                o[INDEX3(i,2,3,numComp,3)] = ((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1565                            } // end of component loop i
1566                        } // end of k0 loop
1567                    } // end of k1 loop
1568                } // end of face 4
1569                if (m_faceOffset[5] > -1) {
1570    #pragma omp for nowait
1571                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1572                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1573                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1574                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1575                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1576                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1577                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1578                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1579                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1580                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1581                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1582                            for (index_t i=0; i < numComp; ++i) {
1583                                const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1584                                const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
1585                                const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1586                                const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
1587                                o[INDEX3(i,0,0,numComp,3)] = V0;
1588                                o[INDEX3(i,1,0,numComp,3)] = V2;
1589                                o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1590                                o[INDEX3(i,0,1,numComp,3)] = V0;
1591                                o[INDEX3(i,1,1,numComp,3)] = V3;
1592                                o[INDEX3(i,2,1,numComp,3)] = ((f_011[i]-f_010[i])*C0 + (f_101[i]-f_100[i])*C5 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1593                                o[INDEX3(i,0,2,numComp,3)] = V1;
1594                                o[INDEX3(i,1,2,numComp,3)] = V2;
1595                                o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / m_dx[2];
1596                                o[INDEX3(i,0,3,numComp,3)] = V1;
1597                                o[INDEX3(i,1,3,numComp,3)] = V3;
1598                                o[INDEX3(i,2,3,numComp,3)] = ((f_001[i]-f_000[i])*C0 + (f_111[i]-f_110[i])*C5 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / m_dx[2];
1599                            } // end of component loop i
1600                        } // end of k0 loop
1601                    } // end of k1 loop
1602                } // end of face 5
1603            } // end of parallel section
1604        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1605            out.requireWrite();
1606    #pragma omp parallel
1607            {
1608                vector<double> f_000(numComp);
1609                vector<double> f_001(numComp);
1610                vector<double> f_010(numComp);
1611                vector<double> f_011(numComp);
1612                vector<double> f_100(numComp);
1613                vector<double> f_101(numComp);
1614                vector<double> f_110(numComp);
1615                vector<double> f_111(numComp);
1616                if (m_faceOffset[0] > -1) {
1617    #pragma omp for nowait
1618                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1619                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1620                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1621                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1622                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1623                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1624                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1625                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1626                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1627                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1628                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1629                            for (index_t i=0; i < numComp; ++i) {
1630                                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];
1631                                o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];
1632                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / m_dx[2];
1633                            } // end of component loop i
1634                        } // end of k1 loop
1635                    } // end of k2 loop
1636                } // end of face 0
1637                if (m_faceOffset[1] > -1) {
1638    #pragma omp for nowait
1639                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1640                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1641                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1642                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1643                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1644                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1645                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1646                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1647                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1648                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1649                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1650                            for (index_t i=0; i < numComp; ++i) {
1651                                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];
1652                                o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];
1653                                o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / m_dx[2];
1654                            } // end of component loop i
1655                        } // end of k1 loop
1656                    } // end of k2 loop
1657                } // end of face 1
1658                if (m_faceOffset[2] > -1) {
1659    #pragma omp for nowait
1660                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1661                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1662                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1663                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1664                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1665                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1666                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1667                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1668                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1669                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1670                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1671                            for (index_t i=0; i < numComp; ++i) {
1672                                o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];
1673                                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];
1674                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / m_dx[2];
1675                            } // end of component loop i
1676                        } // end of k0 loop
1677                    } // end of k2 loop
1678                } // end of face 2
1679                if (m_faceOffset[3] > -1) {
1680    #pragma omp for nowait
1681                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
1682                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1683                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1684                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1685                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1686                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1687                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1688                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1689                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1690                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1691                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1692                            for (index_t i=0; i < numComp; ++i) {
1693                                o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];
1694                                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];
1695                                o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / m_dx[2];
1696                            } // end of component loop i
1697                        } // end of k0 loop
1698                    } // end of k2 loop
1699                } // end of face 3
1700                if (m_faceOffset[4] > -1) {
1701    #pragma omp for nowait
1702                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1703                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1704                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1705                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1706                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1707                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1708                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1709                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1710                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1711                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1712                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1713                            for (index_t i=0; i < numComp; ++i) {
1714                                o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];
1715                                o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];
1716                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C4 / m_dx[2];
1717                            } // end of component loop i
1718                        } // end of k0 loop
1719                    } // end of k1 loop
1720                } // end of face 4
1721                if (m_faceOffset[5] > -1) {
1722    #pragma omp for nowait
1723                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
1724                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
1725                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1726                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1727                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1728                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1729                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1730                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1731                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1732                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1733                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1734                            for (index_t i=0; i < numComp; ++i) {
1735                                o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];
1736                                o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];
1737                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / m_dx[2];
1738                            } // end of component loop i
1739                        } // end of k0 loop
1740                    } // end of k1 loop
1741                } // end of face 5
1742            } // end of parallel section
1743        }
1744  }  }
1745    
1746  pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const  //protected
1747    void Brick::assembleIntegrate(vector<double>& integrals, const escript::Data& arg) const
1748  {  {
1749      if (dim==0)      const dim_t numComp = arg.getDataPointSize();
1750          return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);      const index_t left = (m_offset[0]==0 ? 0 : 1);
1751      else if (dim==1)      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1752          return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);      const index_t front = (m_offset[2]==0 ? 0 : 1);
1753      else if (dim==2)      const int fs = arg.getFunctionSpace().getTypeCode();
1754          return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);      if (fs == Elements && arg.actsExpanded()) {
1755            const double w_0 = m_dx[0]*m_dx[1]*m_dx[2]/8.;
1756    #pragma omp parallel
1757            {
1758                vector<double> int_local(numComp, 0);
1759    #pragma omp for nowait
1760                for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1761                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1762                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1763                            const double* f = arg.getSampleDataRO(INDEX3(k0, k1, k2, m_NE[0], m_NE[1]));
1764                            for (index_t i=0; i < numComp; ++i) {
1765                                const double f_0 = f[INDEX2(i,0,numComp)];
1766                                const double f_1 = f[INDEX2(i,1,numComp)];
1767                                const double f_2 = f[INDEX2(i,2,numComp)];
1768                                const double f_3 = f[INDEX2(i,3,numComp)];
1769                                const double f_4 = f[INDEX2(i,4,numComp)];
1770                                const double f_5 = f[INDEX2(i,5,numComp)];
1771                                const double f_6 = f[INDEX2(i,6,numComp)];
1772                                const double f_7 = f[INDEX2(i,7,numComp)];
1773                                int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
1774                            }  // end of component loop i
1775                        } // end of k0 loop
1776                    } // end of k1 loop
1777                } // end of k2 loop
1778    
1779    #pragma omp critical
1780                for (index_t i=0; i<numComp; i++)
1781                    integrals[i]+=int_local[i];
1782            } // end of parallel section
1783    
1784      throw RipleyException("getFirstCoordAndSpacing(): invalid argument");      } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1785  }          const double w_0 = m_dx[0]*m_dx[1]*m_dx[2];
1786    #pragma omp parallel
1787            {
1788                vector<double> int_local(numComp, 0);
1789    #pragma omp for nowait
1790                for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1791                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1792                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1793                            const double* f = arg.getSampleDataRO(INDEX3(k0, k1, k2, m_NE[0], m_NE[1]));
1794                            for (index_t i=0; i < numComp; ++i) {
1795                                int_local[i]+=f[i]*w_0;
1796                            }  // end of component loop i
1797                        } // end of k0 loop
1798                    } // end of k1 loop
1799                } // end of k2 loop
1800    
1801    #pragma omp critical
1802                for (index_t i=0; i<numComp; i++)
1803                    integrals[i]+=int_local[i];
1804            } // end of parallel section
1805    
1806        } else if (fs == FaceElements && arg.actsExpanded()) {
1807            const double w_0 = m_dx[1]*m_dx[2]/4.;
1808            const double w_1 = m_dx[0]*m_dx[2]/4.;
1809            const double w_2 = m_dx[0]*m_dx[1]/4.;
1810    #pragma omp parallel
1811            {
1812                vector<double> int_local(numComp, 0);
1813                if (m_faceOffset[0] > -1) {
1814    #pragma omp for nowait
1815                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1816                        for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1817                            const double* f = arg.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1818                            for (index_t i=0; i < numComp; ++i) {
1819                                const double f_0 = f[INDEX2(i,0,numComp)];
1820                                const double f_1 = f[INDEX2(i,1,numComp)];
1821                                const double f_2 = f[INDEX2(i,2,numComp)];
1822                                const double f_3 = f[INDEX2(i,3,numComp)];
1823                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1824                            }  // end of component loop i
1825                        } // end of k1 loop
1826                    } // end of k2 loop
1827                }
1828    
1829                if (m_faceOffset[1] > -1) {
1830    #pragma omp for nowait
1831                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1832                        for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1833                            const double* f = arg.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1834                            for (index_t i=0; i < numComp; ++i) {
1835                                const double f_0 = f[INDEX2(i,0,numComp)];
1836                                const double f_1 = f[INDEX2(i,1,numComp)];
1837                                const double f_2 = f[INDEX2(i,2,numComp)];
1838                                const double f_3 = f[INDEX2(i,3,numComp)];
1839                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1840                            }  // end of component loop i
1841                        } // end of k1 loop
1842                    } // end of k2 loop
1843                }
1844    
1845                if (m_faceOffset[2] > -1) {
1846    #pragma omp for nowait
1847                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1848                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1849                            const double* f = arg.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1850                            for (index_t i=0; i < numComp; ++i) {
1851                                const double f_0 = f[INDEX2(i,0,numComp)];
1852                                const double f_1 = f[INDEX2(i,1,numComp)];
1853                                const double f_2 = f[INDEX2(i,2,numComp)];
1854                                const double f_3 = f[INDEX2(i,3,numComp)];
1855                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1856                            }  // end of component loop i
1857                        } // end of k1 loop
1858                    } // end of k2 loop
1859                }
1860    
1861                if (m_faceOffset[3] > -1) {
1862    #pragma omp for nowait
1863                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1864                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1865                            const double* f = arg.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1866                            for (index_t i=0; i < numComp; ++i) {
1867                                const double f_0 = f[INDEX2(i,0,numComp)];
1868                                const double f_1 = f[INDEX2(i,1,numComp)];
1869                                const double f_2 = f[INDEX2(i,2,numComp)];
1870                                const double f_3 = f[INDEX2(i,3,numComp)];
1871                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1872                            }  // end of component loop i
1873                        } // end of k1 loop
1874                    } // end of k2 loop
1875                }
1876    
1877                if (m_faceOffset[4] > -1) {
1878    #pragma omp for nowait
1879                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1880                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1881                            const double* f = arg.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1882                            for (index_t i=0; i < numComp; ++i) {
1883                                const double f_0 = f[INDEX2(i,0,numComp)];
1884                                const double f_1 = f[INDEX2(i,1,numComp)];
1885                                const double f_2 = f[INDEX2(i,2,numComp)];
1886                                const double f_3 = f[INDEX2(i,3,numComp)];
1887                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1888                            }  // end of component loop i
1889                        } // end of k1 loop
1890                    } // end of k2 loop
1891                }
1892    
1893                if (m_faceOffset[5] > -1) {
1894    #pragma omp for nowait
1895                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1896                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1897                            const double* f = arg.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1898                            for (index_t i=0; i < numComp; ++i) {
1899                                const double f_0 = f[INDEX2(i,0,numComp)];
1900                                const double f_1 = f[INDEX2(i,1,numComp)];
1901                                const double f_2 = f[INDEX2(i,2,numComp)];
1902                                const double f_3 = f[INDEX2(i,3,numComp)];
1903                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1904                            }  // end of component loop i
1905                        } // end of k1 loop
1906                    } // end of k2 loop
1907                }
1908    
1909    #pragma omp critical
1910                for (index_t i=0; i<numComp; i++)
1911                    integrals[i]+=int_local[i];
1912            } // end of parallel section
1913    
1914        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1915            const double w_0 = m_dx[1]*m_dx[2];
1916            const double w_1 = m_dx[0]*m_dx[2];
1917            const double w_2 = m_dx[0]*m_dx[1];
1918    #pragma omp parallel
1919            {
1920                vector<double> int_local(numComp, 0);
1921                if (m_faceOffset[0] > -1) {
1922    #pragma omp for nowait
1923                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1924                        for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1925                            const double* f = arg.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1926                            for (index_t i=0; i < numComp; ++i) {
1927                                int_local[i]+=f[i]*w_0;
1928                            }  // end of component loop i
1929                        } // end of k1 loop
1930                    } // end of k2 loop
1931                }
1932    
1933                if (m_faceOffset[1] > -1) {
1934    #pragma omp for nowait
1935                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1936                        for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1937                            const double* f = arg.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1938                            for (index_t i=0; i < numComp; ++i) {
1939                                int_local[i]+=f[i]*w_0;
1940                            }  // end of component loop i
1941                        } // end of k1 loop
1942                    } // end of k2 loop
1943                }
1944    
1945                if (m_faceOffset[2] > -1) {
1946    #pragma omp for nowait
1947                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1948                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1949                            const double* f = arg.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1950                            for (index_t i=0; i < numComp; ++i) {
1951                                int_local[i]+=f[i]*w_1;
1952                            }  // end of component loop i
1953                        } // end of k1 loop
1954                    } // end of k2 loop
1955                }
1956    
1957                if (m_faceOffset[3] > -1) {
1958    #pragma omp for nowait
1959                    for (index_t k2 = front; k2 < front+m_ownNE[2]; ++k2) {
1960                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1961                            const double* f = arg.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1962                            for (index_t i=0; i < numComp; ++i) {
1963                                int_local[i]+=f[i]*w_1;
1964                            }  // end of component loop i
1965                        } // end of k1 loop
1966                    } // end of k2 loop
1967                }
1968    
1969                if (m_faceOffset[4] > -1) {
1970    #pragma omp for nowait
1971                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1972                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1973                            const double* f = arg.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1974                            for (index_t i=0; i < numComp; ++i) {
1975                                int_local[i]+=f[i]*w_2;
1976                            }  // end of component loop i
1977                        } // end of k1 loop
1978                    } // end of k2 loop
1979                }
1980    
1981                if (m_faceOffset[5] > -1) {
1982    #pragma omp for nowait
1983                    for (index_t k1 = bottom; k1 < bottom+m_ownNE[1]; ++k1) {
1984                        for (index_t k0 = left; k0 < left+m_ownNE[0]; ++k0) {
1985                            const double* f = arg.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1986                            for (index_t i=0; i < numComp; ++i) {
1987                                int_local[i]+=f[i]*w_2;
1988                            }  // end of component loop i
1989                        } // end of k1 loop
1990                    } // end of k2 loop
1991                }
1992    
1993    #pragma omp critical
1994                for (index_t i=0; i<numComp; i++)
1995                    integrals[i]+=int_local[i];
1996            } // end of parallel section
1997        } // function space selector
1998    }
1999    
2000  //protected  //protected
2001  dim_t Brick::getNumFaceElements() const  dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const
2002  {  {
2003      dim_t n=0;      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2004      //left      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2005      if (m_offset0==0)      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2006          n+=m_NE1*m_NE2;      const int x=node%nDOF0;
2007      //right      const int y=node%(nDOF0*nDOF1)/nDOF0;
2008      if (m_mpiInfo->rank%m_NX==m_NX-1)      const int z=node/(nDOF0*nDOF1);
2009          n+=m_NE1*m_NE2;      int num=0;
2010      //bottom      // loop through potential neighbours and add to index if positions are
2011      if (m_offset1==0)      // within bounds
2012          n+=m_NE0*m_NE2;      for (int i2=-1; i2<2; i2++) {
2013      //top          for (int i1=-1; i1<2; i1++) {
2014      if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)              for (int i0=-1; i0<2; i0++) {
2015          n+=m_NE0*m_NE2;                  // skip node itself
2016      //front                  if (i0==0 && i1==0 && i2==0)
2017      if (m_offset2==0)                      continue;
2018          n+=m_NE0*m_NE1;                  // location of neighbour node
2019      //back                  const int nx=x+i0;
2020      if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)                  const int ny=y+i1;
2021          n+=m_NE0*m_NE1;                  const int nz=z+i2;
2022                    if (nx>=0 && ny>=0 && nz>=0
2023                            && nx<nDOF0 && ny<nDOF1 && nz<nDOF2) {
2024                        index.push_back(nz*nDOF0*nDOF1+ny*nDOF0+nx);
2025                        num++;
2026                    }
2027                }
2028            }
2029        }
2030    
2031      return n;      return num;
2032  }  }
2033    
2034  //protected  //protected
2035  void Brick::assembleCoordinates(escript::Data& arg) const  void Brick::nodesToDOF(escript::Data& out, const escript::Data& in) const
2036  {  {
2037      escriptDataC x = arg.getDataC();      const dim_t numComp = in.getDataPointSize();
2038      int numDim = m_numDim;      out.requireWrite();
     if (!isDataPointShapeEqual(&x, 1, &numDim))  
         throw RipleyException("setToX: Invalid Data object shape");  
     if (!numSamplesEqual(&x, 1, getNumNodes()))  
         throw RipleyException("setToX: Illegal number of samples in Data object");  
2039    
2040      pair<double,double> xdx = getFirstCoordAndSpacing(0);      const index_t left = (m_offset[0]==0 ? 0 : 1);
2041      pair<double,double> ydy = getFirstCoordAndSpacing(1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
2042      pair<double,double> zdz = getFirstCoordAndSpacing(2);      const index_t front = (m_offset[2]==0 ? 0 : 1);
2043      arg.requireWrite();      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2044        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2045        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2046  #pragma omp parallel for  #pragma omp parallel for
2047      for (dim_t i2 = 0; i2 < m_N2; i2++) {      for (index_t i=0; i<nDOF2; i++) {
2048          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (index_t j=0; j<nDOF1; j++) {
2049              for (dim_t i0 = 0; i0 < m_N0; i0++) {              for (index_t k=0; k<nDOF0; k++) {
2050                  double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);                  const index_t n=k+left+(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];
2051                  point[0] = xdx.first+i0*xdx.second;                  const double* src=in.getSampleDataRO(n);
2052                  point[1] = ydy.first+i1*ydy.second;                  copy(src, src+numComp, out.getSampleDataRW(k+j*nDOF0+i*nDOF0*nDOF1));
                 point[2] = zdz.first+i2*zdz.second;  
2053              }              }
2054          }          }
2055      }      }
2056  }  }
2057    
2058    //protected
2059    void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
2060    {
2061        const dim_t numComp = in.getDataPointSize();
2062        paso::Coupler_ptr coupler(new paso::Coupler(m_connector, numComp));
2063        // expand data object if necessary to be able to grab the whole data
2064        const_cast<escript::Data*>(&in)->expand();
2065        coupler->startCollect(in.getSampleDataRO(0));
2066    
2067        const dim_t numDOF = getNumDOF();
2068        out.requireWrite();
2069        const double* buffer = coupler->finishCollect();
2070    
2071    #pragma omp parallel for
2072        for (index_t i=0; i<getNumNodes(); i++) {
2073            const double* src=(m_dofMap[i]<numDOF ?
2074                    in.getSampleDataRO(m_dofMap[i])
2075                    : &buffer[(m_dofMap[i]-numDOF)*numComp]);
2076            copy(src, src+numComp, out.getSampleDataRW(i));
2077        }
2078    }
2079    
2080  //private  //private
2081  void Brick::populateSampleIds()  void Brick::populateSampleIds()
2082  {  {
2083      // identifiers are ordered from left to right, bottom to top, front to back      // degrees of freedom are numbered from left to right, bottom to top, front
2084      // on each rank, except for the shared nodes which are owned by the rank      // to back in each rank, continuing on the next rank (ranks also go
2085      // below / to the left / to the front of the current rank      // left-right, bottom-top, front-back).
2086        // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
2087        // helps when writing out data rank after rank.
2088    
2089      // build node distribution vector first.      // build node distribution vector first.
2090      // m_nodeDistribution[i] is the first node id on rank i, that is      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
2091      // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes      // constant for all ranks in this implementation
2092      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);      m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
2093      m_nodeDistribution[1]=getNumNodes();      const dim_t numDOF=getNumDOF();
2094      for (dim_t k=1; k<m_mpiInfo->size-1; k++) {      for (dim_t k=1; k<m_mpiInfo->size; k++) {
2095          const index_t x = k%m_NX;          m_nodeDistribution[k]=k*numDOF;
         const index_t y = k%(m_NX*m_NY)/m_NX;  
         const index_t z = k/(m_NX*m_NY);  
         index_t numNodes=getNumNodes();  
         if (x>0)  
             numNodes-=m_N1*m_N2;  
         if (y>0)  
             numNodes-=m_N0*m_N2;  
         if (z>0)  
             numNodes-=m_N0*m_N1;  
         // if an edge was subtracted twice add it back  
         if (x>0 && y>0)  
             numNodes+=m_N2;  
         if (x>0 && z>0)  
             numNodes+=m_N1;  
         if (y>0 && z>0)  
             numNodes+=m_N0;  
         // the corner node was removed 3x and added back 3x, so subtract it  
         if (x>0 && y>0 && z>0)  
             numNodes--;  
         m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;  
2096      }      }
2097      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();      m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
2098        
2099        try {
2100            m_nodeId.resize(getNumNodes());
2101            m_dofId.resize(numDOF);
2102            m_elementId.resize(getNumElements());
2103        } catch (const std::length_error& le) {
2104            throw RipleyException("The system does not have sufficient memory for a domain of this size.");
2105        }
2106        
2107        // populate face element counts
2108        //left
2109        if (m_offset[0]==0)
2110            m_faceCount[0]=m_NE[1]*m_NE[2];
2111        else
2112            m_faceCount[0]=0;
2113        //right
2114        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
2115            m_faceCount[1]=m_NE[1]*m_NE[2];
2116        else
2117            m_faceCount[1]=0;
2118        //bottom
2119        if (m_offset[1]==0)
2120            m_faceCount[2]=m_NE[0]*m_NE[2];
2121        else
2122            m_faceCount[2]=0;
2123        //top
2124        if (m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0]==m_NX[1]-1)
2125            m_faceCount[3]=m_NE[0]*m_NE[2];
2126        else
2127            m_faceCount[3]=0;
2128        //front
2129        if (m_offset[2]==0)
2130            m_faceCount[4]=m_NE[0]*m_NE[1];
2131        else
2132            m_faceCount[4]=0;
2133        //back
2134        if (m_mpiInfo->rank/(m_NX[0]*m_NX[1])==m_NX[2]-1)
2135            m_faceCount[5]=m_NE[0]*m_NE[1];
2136        else
2137            m_faceCount[5]=0;
2138    
2139      m_nodeId.resize(getNumNodes());      m_faceId.resize(getNumFaceElements());
2140    
2141      // the bottom, left and front planes are not owned by this rank so the      const index_t left = (m_offset[0]==0 ? 0 : 1);
2142      // identifiers need to be computed accordingly      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
2143      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t front = (m_offset[2]==0 ? 0 : 1);
2144      const index_t bottom = (m_offset1==0 ? 0 : 1);      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2145      const index_t front = (m_offset2==0 ? 0 : 1);      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2146        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2147      // case 1: all nodes on left plane are owned by rank on the left  
2148      if (left>0) {      // the following is a compromise between efficiency and code length to
2149          const int neighbour=m_mpiInfo->rank-1;      // set the node id's according to the order mentioned above.
2150          const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // First we set all the edge and corner id's in a rather slow way since
2151          const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);      // they might or might not be owned by this rank. Next come the own
2152  #pragma omp parallel for      // node id's which are identical to the DOF id's (simple loop), and finally
2153          for (dim_t i2=front; i2<m_N2; i2++) {      // the 6 faces are set but only if required...
2154              for (dim_t i1=bottom; i1<m_N1; i1++) {  
2155                  m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]  #define globalNodeId(x,y,z) \
2156                      + (i1-bottom+1)*leftN0      ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1*nDOF2+(m_offset[0]+x)%nDOF0\
2157                      + (i2-front)*leftN0*leftN1 - 1;      + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*nDOF2*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0\
2158        + ((m_offset[2]+z)/nDOF2)*nDOF0*nDOF1*nDOF2*m_NX[0]*m_NX[1]+((m_offset[2]+z)%nDOF2)*nDOF0*nDOF1
2159    
2160    #pragma omp parallel
2161        {
2162            // set edge id's
2163            // edges in x-direction, including corners
2164    #pragma omp for nowait
2165            for (dim_t i=0; i<m_NN[0]; i++) {
2166                m_nodeId[i] = globalNodeId(i, 0, 0); // LF
2167                m_nodeId[m_NN[0]*(m_NN[1]-1)+i] = globalNodeId(i, m_NN[1]-1, 0); // UF
2168                m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+i] = globalNodeId(i, 0, m_NN[2]-1); // LB
2169                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
2170            }
2171            // edges in y-direction, without corners
2172    #pragma omp for nowait
2173            for (dim_t i=1; i<m_NN[1]-1; i++) {
2174                m_nodeId[m_NN[0]*i] = globalNodeId(0, i, 0); // FL
2175                m_nodeId[m_NN[0]*(i+1)-1] = globalNodeId(m_NN[0]-1, i, 0); // FR
2176                m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+m_NN[0]*i] = globalNodeId(0, i, m_NN[2]-1); // BL
2177                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
2178            }
2179            // edges in z-direction, without corners
2180    #pragma omp for
2181            for (dim_t i=1; i<m_NN[2]-1; i++) {
2182                m_nodeId[m_NN[0]*m_NN[1]*i] = globalNodeId(0, 0, i); // LL
2183                m_nodeId[m_NN[0]*m_NN[1]*i+m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0, i); // LR
2184                m_nodeId[m_NN[0]*m_NN[1]*(i+1)-m_NN[0]] = globalNodeId(0, m_NN[1]-1, i); // UL
2185                m_nodeId[m_NN[0]*m_NN[1]*(i+1)-1] = globalNodeId(m_NN[0]-1, m_NN[1]-1, i); // UR
2186            }
2187            // implicit barrier here because some node IDs will be overwritten
2188            // below
2189    
2190            // populate degrees of freedom and own nodes (identical id)
2191    #pragma omp for nowait
2192            for (dim_t i=0; i<nDOF2; i++) {
2193                for (dim_t j=0; j<nDOF1; j++) {
2194                    for (dim_t k=0; k<nDOF0; k++) {
2195                        const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];
2196                        const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;
2197                        m_dofId[dofIdx] = m_nodeId[nodeIdx]
2198                            = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
2199                    }
2200                }
2201            }
2202    
2203            // populate the rest of the nodes (shared with other ranks)
2204            if (m_faceCount[0]==0) { // left plane
2205    #pragma omp for nowait
2206                for (dim_t i=0; i<nDOF2; i++) {
2207                    for (dim_t j=0; j<nDOF1; j++) {
2208                        const index_t nodeIdx=(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];
2209                        const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;
2210                        m_nodeId[nodeIdx]
2211                            = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
2212                    }
2213                }
2214            }
2215            if (m_faceCount[1]==0) { // right plane
2216    #pragma omp for nowait
2217                for (dim_t i=0; i<nDOF2; i++) {
2218                    for (dim_t j=0; j<nDOF1; j++) {
2219                        const index_t nodeIdx=(j+bottom+1)*m_NN[0]-1+(i+front)*m_NN[0]*m_NN[1];
2220                        const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;
2221                        m_nodeId[nodeIdx]
2222                            = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
2223                    }
2224                }
2225            }
2226            if (m_faceCount[2]==0) { // bottom plane
2227    #pragma omp for nowait
2228                for (dim_t i=0; i<nDOF2; i++) {
2229                    for (dim_t k=0; k<nDOF0; k++) {
2230                        const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1];
2231                        const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;
2232                        m_nodeId[nodeIdx]
2233                            = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
2234                    }
2235                }
2236            }
2237            if (m_faceCount[3]==0) { // top plane
2238    #pragma omp for nowait
2239                for (dim_t i=0; i<nDOF2; i++) {
2240                    for (dim_t k=0; k<nDOF0; k++) {
2241                        const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1]+m_NN[0]*(m_NN[1]-1);
2242                        const index_t dofId=k+i*nDOF0*nDOF1;
2243                        m_nodeId[nodeIdx]
2244                            = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
2245                    }
2246                }
2247            }
2248            if (m_faceCount[4]==0) { // front plane
2249    #pragma omp for nowait
2250                for (dim_t j=0; j<nDOF1; j++) {
2251                    for (dim_t k=0; k<nDOF0; k++) {
2252                        const index_t nodeIdx=k+left+(j+bottom)*m_NN[0];
2253                        const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);
2254                        m_nodeId[nodeIdx]
2255                            = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;
2256                    }
2257                }
2258            }
2259            if (m_faceCount[5]==0) { // back plane
2260    #pragma omp for nowait
2261                for (dim_t j=0; j<nDOF1; j++) {
2262                    for (dim_t k=0; k<nDOF0; k++) {
2263                        const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1);
2264                        const index_t dofId=k+j*nDOF0;
2265                        m_nodeId[nodeIdx]
2266                            = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;
2267                    }
2268                }
2269            }
2270    
2271            // populate element id's
2272    #pragma omp for nowait
2273            for (dim_t i2=0; i2<m_NE[2]; i2++) {
2274                for (dim_t i1=0; i1<m_NE[1]; i1++) {
2275                    for (dim_t i0=0; i0<m_NE[0]; i0++) {
2276                        m_elementId[i0+i1*m_NE[0]+i2*m_NE[0]*m_NE[1]] =
2277                            (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]
2278                            +(m_offset[1]+i1)*m_gNE[0]
2279                            +m_offset[0]+i0;
2280                    }
2281              }              }
2282          }          }
2283    
2284            // face elements
2285    #pragma omp for
2286            for (dim_t k=0; k<getNumFaceElements(); k++)
2287                m_faceId[k]=k;
2288        } // end parallel section
2289    
2290    #undef globalNodeId
2291    
2292        m_nodeTags.assign(getNumNodes(), 0);
2293        updateTagsInUse(Nodes);
2294    
2295        m_elementTags.assign(getNumElements(), 0);
2296        updateTagsInUse(Elements);
2297    
2298        // generate face offset vector and set face tags
2299        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
2300        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
2301        m_faceOffset.assign(6, -1);
2302        m_faceTags.clear();
2303        index_t offset=0;
2304        for (size_t i=0; i<6; i++) {
2305            if (m_faceCount[i]>0) {
2306                m_faceOffset[i]=offset;
2307                offset+=m_faceCount[i];
2308                m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
2309            }
2310      }      }
2311      // case 2: all nodes on bottom plane are owned by rank below      setTagMap("left", LEFT);
2312      if (bottom>0) {      setTagMap("right", RIGHT);
2313          const int neighbour=m_mpiInfo->rank-m_NX;      setTagMap("bottom", BOTTOM);
2314          const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      setTagMap("top", TOP);
2315          const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);      setTagMap("front", FRONT);
2316        setTagMap("back", BACK);
2317        updateTagsInUse(FaceElements);
2318    }
2319    
2320    //private
2321    void Brick::createPattern()
2322    {
2323        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2324        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2325        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2326        const index_t left = (m_offset[0]==0 ? 0 : 1);
2327        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
2328        const index_t front = (m_offset[2]==0 ? 0 : 1);
2329    
2330        // populate node->DOF mapping with own degrees of freedom.
2331        // The rest is assigned in the loop further down
2332        m_dofMap.assign(getNumNodes(), 0);
2333  #pragma omp parallel for  #pragma omp parallel for
2334          for (dim_t i2=front; i2<m_N2; i2++) {      for (index_t i=front; i<front+nDOF2; i++) {
2335              for (dim_t i0=left; i0<m_N0; i0++) {          for (index_t j=bottom; j<bottom+nDOF1; j++) {
2336                  m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]              for (index_t k=left; k<left+nDOF0; k++) {
2337                      + bottomN0*(bottomN1-1)                  m_dofMap[i*m_NN[0]*m_NN[1]+j*m_NN[0]+k]=(i-front)*nDOF0*nDOF1+(j-bottom)*nDOF0+k-left;
                     + (i2-front)*bottomN0*bottomN1 + i0-left;  
2338              }              }
2339          }          }
2340      }      }
2341      // case 3: all nodes on front plane are owned by rank in front  
2342      if (front>0) {      // build list of shared components and neighbours by looping through
2343          const int neighbour=m_mpiInfo->rank-m_NX*m_NY;      // all potential neighbouring ranks and checking if positions are
2344          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      // within bounds
2345          const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);      const dim_t numDOF=nDOF0*nDOF1*nDOF2;
2346          const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);      vector<IndexVector> colIndices(numDOF); // for the couple blocks
2347  #pragma omp parallel for      RankVector neighbour;
2348          for (dim_t i1=bottom; i1<m_N1; i1++) {      IndexVector offsetInShared(1,0);
2349              for (dim_t i0=left; i0<m_N0; i0++) {      IndexVector sendShared, recvShared;
2350                  m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]      int numShared=0, expectedShared=0;;
2351                      + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;      const int x=m_mpiInfo->rank%m_NX[0];
2352        const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2353        const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
2354        for (int i2=-1; i2<2; i2++) {
2355            for (int i1=-1; i1<2; i1++) {
2356                for (int i0=-1; i0<2; i0++) {
2357                    // skip this rank
2358                    if (i0==0 && i1==0 && i2==0)
2359                        continue;
2360                    // location of neighbour rank
2361                    const int nx=x+i0;
2362                    const int ny=y+i1;
2363                    const int nz=z+i2;
2364                    if (!(nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2])) {
2365                        continue;
2366                    }
2367                    if (i0==0 && i1==0)
2368                        expectedShared += nDOF0*nDOF1;
2369                    else if (i0==0 && i2==0)
2370                        expectedShared += nDOF0*nDOF2;
2371                    else if (i1==0 && i2==0)
2372                        expectedShared += nDOF1*nDOF2;
2373                    else if (i0==0)
2374                        expectedShared += nDOF0;
2375                    else if (i1==0)
2376                        expectedShared += nDOF1;
2377                    else if (i2==0)
2378                        expectedShared += nDOF2;
2379                    else
2380                        expectedShared++;
2381              }              }
2382          }          }
2383      }      }
2384      // case 4: nodes on front bottom edge are owned by the corresponding rank      
2385      if (front>0 && bottom>0) {      vector<IndexVector> rowIndices(expectedShared);
2386          const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);      
2387          const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);      for (int i2=-1; i2<2; i2++) {
2388          const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);          for (int i1=-1; i1<2; i1++) {
2389          const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);              for (int i0=-1; i0<2; i0++) {
2390  #pragma omp parallel for                  // skip this rank
2391          for (dim_t i0=left; i0<m_N0; i0++) {                  if (i0==0 && i1==0 && i2==0)
2392              m_nodeId[i0]=m_nodeDistribution[neighbour]                      continue;
2393                  + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;                  // location of neighbour rank
2394                    const int nx=x+i0;
2395                    const int ny=y+i1;
2396                    const int nz=z+i2;
2397                    if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {
2398                        neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);
2399                        if (i0==0 && i1==0) {
2400                            // sharing front or back plane
2401                            offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF1);
2402                            for (dim_t i=0; i<nDOF1; i++) {
2403                                const int firstDOF=(i2==-1 ? i*nDOF0
2404                                        : i*nDOF0 + nDOF0*nDOF1*(nDOF2-1));
2405                                const int firstNode=(i2==-1 ? left+(i+bottom)*m_NN[0]
2406                                        : left+(i+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1));
2407                                for (dim_t j=0; j<nDOF0; j++, numShared++) {
2408                                    sendShared.push_back(firstDOF+j);
2409                                    recvShared.push_back(numDOF+numShared);
2410                                    if (j>0) {
2411                                        if (i>0)
2412                                            doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0, numShared);
2413                                        doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2414                                        if (i<nDOF1-1)
2415                                            doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0, numShared);
2416                                    }
2417                                    if (i>0)
2418                                        doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0, numShared);
2419                                    doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2420                                    if (i<nDOF1-1)
2421                                        doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0, numShared);
2422                                    if (j<nDOF0-1) {
2423                                        if (i>0)
2424                                            doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0, numShared);
2425                                        doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2426                                        if (i<nDOF1-1)
2427                                            doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0, numShared);
2428                                    }
2429                                    m_dofMap[firstNode+j]=numDOF+numShared;
2430                                }
2431                            }
2432                        } else if (i0==0 && i2==0) {
2433                            // sharing top or bottom plane
2434                            offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF2);
2435                            for (dim_t i=0; i<nDOF2; i++) {
2436                                const int firstDOF=(i1==-1 ? i*nDOF0*nDOF1
2437                                        : nDOF0*((i+1)*nDOF1-1));
2438                                const int firstNode=(i1==-1 ?
2439                                        left+(i+front)*m_NN[0]*m_NN[1]
2440                                        : left+m_NN[0]*((i+1+front)*m_NN[1]-1));
2441                                for (dim_t j=0; j<nDOF0; j++, numShared++) {
2442                                    sendShared.push_back(firstDOF+j);
2443                                    recvShared.push_back(numDOF+numShared);
2444                                    if (j>0) {
2445                                        if (i>0)
2446                                            doublyLink(colIndices, rowIndices, firstDOF+j-1-nDOF0*nDOF1, numShared);
2447                                        doublyLink(colIndices, rowIndices, firstDOF+j-1, numShared);
2448                                        if (i<nDOF2-1)
2449                                            doublyLink(colIndices, rowIndices, firstDOF+j-1+nDOF0*nDOF1, numShared);
2450                                    }
2451                                    if (i>0)
2452                                        doublyLink(colIndices, rowIndices, firstDOF+j-nDOF0*nDOF1, numShared);
2453                                    doublyLink(colIndices, rowIndices, firstDOF+j, numShared);
2454                                    if (i<nDOF2-1)
2455                                        doublyLink(colIndices, rowIndices, firstDOF+j+nDOF0*nDOF1, numShared);
2456                                    if (j<nDOF0-1) {
2457                                        if (i>0)
2458                                            doublyLink(colIndices, rowIndices, firstDOF+j+1-nDOF0*nDOF1, numShared);
2459                                        doublyLink(colIndices, rowIndices, firstDOF+j+1, numShared);
2460                                        if (i<nDOF2-1)
2461                                            doublyLink(colIndices, rowIndices, firstDOF+j+1+nDOF0*nDOF1, numShared);
2462                                    }
2463                                    m_dofMap[firstNode+j]=numDOF+numShared;
2464                                }
2465                            }
2466                        } else if (i1==0 && i2==0) {
2467                            // sharing left or right plane
2468                            offsetInShared.push_back(offsetInShared.back()+nDOF1*nDOF2);
2469                            for (dim_t i=0; i<nDOF2; i++) {
2470                                const int firstDOF=(i0==-1 ? i*nDOF0*nDOF1
2471                                        : nDOF0*(1+i*nDOF1)-1);
2472                                const int firstNode=(i0==-1 ?
2473                                        (bottom+(i+front)*m_NN[1])*m_NN[0]
2474                                        : (bottom+1+(i+front)*m_NN[1])*m_NN[0]-1);
2475                                for (dim_t j=0; j<nDOF1; j++, numShared++) {
2476                                    sendShared.push_back(firstDOF+j*nDOF0);
2477                                    recvShared.push_back(numDOF+numShared);
2478                                    if (j>0) {
2479                                        if (i>0)
2480                                            doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0-nDOF0*nDOF1, numShared);
2481                                        doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0, numShared);
2482                                        if (i<nDOF2-1)
2483                                            doublyLink(colIndices, rowIndices, firstDOF+(j-1)*nDOF0+nDOF0*nDOF1, numShared);
2484                                    }
2485                                    if (i>0)
2486                                        doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0-nDOF0*nDOF1, numShared);
2487                                    doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0, numShared);
2488                                    if (i<nDOF2-1)
2489                                        doublyLink(colIndices, rowIndices, firstDOF+j*nDOF0+nDOF0*nDOF1, numShared);
2490                                    if (j<nDOF1-1) {
2491                                        if (i>0)
2492                                            doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0-nDOF0*nDOF1, numShared);
2493                                        doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0, numShared);
2494                                        if (i<nDOF2-1)
2495                                            doublyLink(colIndices, rowIndices, firstDOF+(j+1)*nDOF0+nDOF0*nDOF1, numShared);
2496                                    }
2497                                    m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
2498                                }
2499                            }
2500                        } else if (i0==0) {
2501                            // sharing an edge in x direction
2502                            offsetInShared.push_back(offsetInShared.back()+nDOF0);
2503                            const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
2504                                               +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2505                            const int firstNode=left+(i1+1)/2*m_NN[0]*(m_NN[1]-1)
2506                                                +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
2507                            for (dim_t i=0; i<nDOF0; i++, numShared++) {
2508                                sendShared.push_back(firstDOF+i);
2509                                recvShared.push_back(numDOF+numShared);
2510                                if (i>0)
2511                                    doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
2512                                doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
2513                                if (i<nDOF0-1)
2514                                    doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
2515                                m_dofMap[firstNode+i]=numDOF+numShared;
2516                            }
2517                        } else if (i1==0) {
2518                            // sharing an edge in y direction
2519                            offsetInShared.push_back(offsetInShared.back()+nDOF1);
2520                            const int firstDOF=(i0+1)/2*(nDOF0-1)
2521                                               +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2522                            const int firstNode=bottom*m_NN[0]
2523                                                +(i0+1)/2*(m_NN[0]-1)
2524                                                +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
2525                            for (dim_t i=0; i<nDOF1; i++, numShared++) {
2526                                sendShared.push_back(firstDOF+i*nDOF0);
2527                                recvShared.push_back(numDOF+numShared);
2528                                if (i>0)
2529                                    doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
2530                                doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
2531                                if (i<nDOF1-1)
2532                                    doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
2533                                m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
2534                            }
2535                        } else if (i2==0) {
2536                            // sharing an edge in z direction
2537                            offsetInShared.push_back(offsetInShared.back()+nDOF2);
2538                            const int firstDOF=(i0+1)/2*(nDOF0-1)
2539                                               +(i1+1)/2*nDOF0*(nDOF1-1);
2540                            const int firstNode=front*m_NN[0]*m_NN[1]
2541                                                +(i0+1)/2*(m_NN[0]-1)
2542                                                +(i1+1)/2*m_NN[0]*(m_NN[1]-1);
2543                            for (dim_t i=0; i<nDOF2; i++, numShared++) {
2544                                sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
2545                                recvShared.push_back(numDOF+numShared);
2546                                if (i>0)
2547                                    doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0*nDOF1, numShared);
2548                                doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0*nDOF1, numShared);
2549                                if (i<nDOF2-1)
2550                                    doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0*nDOF1, numShared);
2551                                m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
2552                            }
2553                        } else {
2554                            // sharing a node
2555                            const int dof=(i0+1)/2*(nDOF0-1)
2556                                          +(i1+1)/2*nDOF0*(nDOF1-1)
2557                                          +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2558                            const int node=(i0+1)/2*(m_NN[0]-1)
2559                                           +(i1+1)/2*m_NN[0]*(m_NN[1]-1)
2560                                           +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
2561                            offsetInShared.push_back(offsetInShared.back()+1);
2562                            sendShared.push_back(dof);
2563                            recvShared.push_back(numDOF+numShared);
2564                            doublyLink(colIndices, rowIndices, dof, numShared);
2565                            m_dofMap[node]=numDOF+numShared;
2566                            ++numShared;
2567                        }
2568                    }
2569                }
2570          }          }
2571      }      }
2572      // case 5: nodes on left bottom edge are owned by the corresponding rank  
     if (left>0 && bottom>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
2573  #pragma omp parallel for  #pragma omp parallel for
2574          for (dim_t i2=front; i2<m_N2; i2++) {      for (int i = 0; i < numShared; i++) {
2575              m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]          std::sort(rowIndices[i].begin(), rowIndices[i].end());
2576                  + (1+i2-front)*N0*N1-1;      }
2577    
2578        // create connector
2579        paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
2580                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2581                &offsetInShared[0], 1, 0, m_mpiInfo));
2582        paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
2583                numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
2584                &offsetInShared[0], 1, 0, m_mpiInfo));
2585        m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
2586    
2587        // create main and couple blocks
2588        paso::Pattern_ptr mainPattern = createMainPattern();
2589        paso::Pattern_ptr colPattern, rowPattern;
2590        createCouplePatterns(colIndices, rowIndices, numShared, colPattern, rowPattern);
2591    
2592        // allocate paso distribution
2593        paso::Distribution_ptr distribution(new paso::Distribution(m_mpiInfo,
2594                const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0));
2595    
2596        // finally create the system matrix
2597        m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
2598                distribution, distribution, mainPattern, colPattern, rowPattern,
2599                m_connector, m_connector));
2600    
2601        // useful debug output
2602        /*
2603        cout << "--- rcv_shcomp ---" << endl;
2604        cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
2605        for (size_t i=0; i<neighbour.size(); i++) {
2606            cout << "neighbor[" << i << "]=" << neighbour[i]
2607                << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
2608        }
2609        for (size_t i=0; i<recvShared.size(); i++) {
2610            cout << "shared[" << i << "]=" << recvShared[i] << endl;
2611        }
2612        cout << "--- snd_shcomp ---" << endl;
2613        for (size_t i=0; i<sendShared.size(); i++) {
2614            cout << "shared[" << i << "]=" << sendShared[i] << endl;
2615        }
2616        cout << "--- dofMap ---" << endl;
2617        for (size_t i=0; i<m_dofMap.size(); i++) {
2618            cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << endl;
2619        }
2620        cout << "--- colIndices ---" << endl;
2621        for (size_t i=0; i<colIndices.size(); i++) {
2622            cout << "colIndices[" << i << "].size()=" << colIndices[i].size() << endl;
2623        }
2624        */
2625    
2626        /*
2627        cout << "--- main_pattern ---" << endl;
2628        cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
2629        for (size_t i=0; i<mainPattern->numOutput+1; i++) {
2630            cout << "ptr[" << i << "]=" << mainPattern->ptr[i] << endl;
2631        }
2632        for (size_t i=0; i<mainPattern->ptr[mainPattern->numOutput]; i++) {
2633            cout << "index[" << i << "]=" << mainPattern->index[i] << endl;
2634        }
2635        */
2636    
2637        /*
2638        cout << "--- colCouple_pattern ---" << endl;
2639        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
2640        for (size_t i=0; i<colPattern->numOutput+1; i++) {
2641            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
2642        }
2643        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
2644            cout << "index[" << i << "]=" << colPattern->index[i] << endl;
2645        }
2646        */
2647    
2648        /*
2649        cout << "--- rowCouple_pattern ---" << endl;
2650        cout << "M=" << rowPattern->numOutput << ", N=" << rowPattern->numInput << endl;
2651        for (size_t i=0; i<rowPattern->numOutput+1; i++) {
2652            cout << "ptr[" << i << "]=" << rowPattern->ptr[i] << endl;
2653        }
2654        for (size_t i=0; i<rowPattern->ptr[rowPattern->numOutput]; i++) {
2655            cout << "index[" << i << "]=" << rowPattern->index[i] << endl;
2656        }
2657        */
2658    }
2659    
2660    //private
2661    void Brick::addToMatrixAndRHS(paso::SystemMatrix_ptr S, escript::Data& F,
2662             const vector<double>& EM_S, const vector<double>& EM_F, bool addS,
2663             bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
2664    {
2665        IndexVector rowIndex;
2666        rowIndex.push_back(m_dofMap[firstNode]);
2667        rowIndex.push_back(m_dofMap[firstNode+1]);
2668        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]]);
2669        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]+1]);
2670        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*m_NN[1]]);
2671        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*m_NN[1]+1]);
2672        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)]);
2673        rowIndex.push_back(m_dofMap[firstNode+m_NN[0]*(m_NN[1]+1)+1]);
2674        if (addF) {
2675            double *F_p=F.getSampleDataRW(0);
2676            for (index_t i=0; i<rowIndex.size(); i++) {
2677                if (rowIndex[i]<getNumDOF()) {
2678                    for (index_t eq=0; eq<nEq; eq++) {
2679                        F_p[INDEX2(eq, rowIndex[i], nEq)]+=EM_F[INDEX2(eq,i,nEq)];
2680                    }
2681                }
2682          }          }
2683      }      }
2684      // case 6: nodes on left front edge are owned by the corresponding rank      if (addS) {
2685      if (left>0 && front>0) {          addToSystemMatrix(S, rowIndex, nEq, rowIndex, nComp, EM_S);
         const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);  
 #pragma omp parallel for  
         for (dim_t i1=bottom; i1<m_N1; i1++) {  
             m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]  
                 + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;  
         }  
     }  
     // case 7: bottom-left-front corner node owned by corresponding rank  
     if (left>0 && bottom>0 && front>0) {  
         const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;  
         const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);  
         const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);  
         const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);  
         m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;  
2686      }      }
2687    }
2688    
2689      // the rest of the id's are contiguous  //protected
2690      const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];  void Brick::interpolateNodesOnElements(escript::Data& out,
2691  #pragma omp parallel for                                         const escript::Data& in,
2692      for (dim_t i2=front; i2<m_N2; i2++) {                                         bool reduced) const
2693          for (dim_t i1=bottom; i1<m_N1; i1++) {  {
2694              for (dim_t i0=left; i0<m_N0; i0++) {      const dim_t numComp = in.getDataPointSize();
2695                  m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left      if (reduced) {
2696                      +(i1-bottom)*(m_N0-left)          out.requireWrite();
2697                      +(i2-front)*(m_N0-left)*(m_N1-bottom);  #pragma omp parallel
2698            {
2699                vector<double> f_000(numComp);
2700                vector<double> f_001(numComp);
2701                vector<double> f_010(numComp);
2702                vector<double> f_011(numComp);
2703                vector<double> f_100(numComp);
2704                vector<double> f_101(numComp);
2705                vector<double> f_110(numComp);
2706                vector<double> f_111(numComp);
2707    #pragma omp for
2708                for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2709                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2710                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2711                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2712                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2713                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2714                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2715                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2716                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2717                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2718                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2719                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
2720                            for (index_t i=0; i < numComp; ++i) {
2721                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i])/8;
2722                            } // end of component loop i
2723                        } // end of k0 loop
2724                    } // end of k1 loop
2725                } // end of k2 loop
2726            } // end of parallel section
2727        } else {
2728            out.requireWrite();
2729            const double c0 = .0094373878376559314545;
2730            const double c1 = .035220810900864519624;
2731            const double c2 = .13144585576580214704;
2732            const double c3 = .49056261216234406855;
2733    #pragma omp parallel
2734            {
2735                vector<double> f_000(numComp);
2736                vector<double> f_001(numComp);
2737                vector<double> f_010(numComp);
2738                vector<double> f_011(numComp);
2739                vector<double> f_100(numComp);
2740                vector<double> f_101(numComp);
2741                vector<double> f_110(numComp);
2742                vector<double> f_111(numComp);
2743    #pragma omp for
2744                for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2745                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2746                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2747                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2748                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2749                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2750                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2751                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2752                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2753                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2754                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2755                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
2756                            for (index_t i=0; i < numComp; ++i) {
2757                                o[INDEX2(i,numComp,0)] = f_000[i]*c3 + f_111[i]*c0 + c2*(f_001[i] + f_010[i] + f_100[i]) + c1*(f_011[i] + f_101[i] + f_110[i]);
2758                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_100[i]*c3 + c2*(f_000[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_111[i]);
2759                                o[INDEX2(i,numComp,2)] = f_010[i]*c3 + f_101[i]*c0 + c2*(f_000[i] + f_011[i] + f_110[i]) + c1*(f_001[i] + f_100[i] + f_111[i]);
2760                                o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_110[i]*c3 + c2*(f_010[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_101[i]);
2761                                o[INDEX2(i,numComp,4)] = f_001[i]*c3 + f_110[i]*c0 + c2*(f_000[i] + f_011[i] + f_101[i]) + c1*(f_010[i] + f_100[i] + f_111[i]);
2762                                o[INDEX2(i,numComp,5)] = f_010[i]*c0 + f_101[i]*c3 + c2*(f_001[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_110[i]);
2763                                o[INDEX2(i,numComp,6)] = f_011[i]*c3 + f_100[i]*c0 + c2*(f_001[i] + f_010[i] + f_111[i]) + c1*(f_000[i] + f_101[i] + f_110[i]);
2764                                o[INDEX2(i,numComp,7)] = f_000[i]*c0 + f_111[i]*c3 + c2*(f_011[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_100[i]);
2765                            } // end of component loop i
2766                        } // end of k0 loop
2767                    } // end of k1 loop
2768                } // end of k2 loop
2769            } // end of parallel section
2770        }
2771    }
2772    
2773    //protected
2774    void Brick::interpolateNodesOnFaces(escript::Data& out, const escript::Data& in,
2775                                        bool reduced) const
2776    {
2777        const dim_t numComp = in.getDataPointSize();
2778        if (reduced) {
2779            out.requireWrite();
2780    #pragma omp parallel
2781            {
2782                vector<double> f_000(numComp);
2783                vector<double> f_001(numComp);
2784                vector<double> f_010(numComp);
2785                vector<double> f_011(numComp);
2786                vector<double> f_100(numComp);
2787                vector<double> f_101(numComp);
2788                vector<double> f_110(numComp);
2789                vector<double> f_111(numComp);
2790                if (m_faceOffset[0] > -1) {
2791    #pragma omp for nowait
2792                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2793                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2794                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2795                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2796                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2797                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2798                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
2799                            for (index_t i=0; i < numComp; ++i) {
2800                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_010[i] + f_011[i])/4;
2801                            } // end of component loop i
2802                        } // end of k1 loop
2803                    } // end of k2 loop
2804                } // end of face 0
2805                if (m_faceOffset[1] > -1) {
2806    #pragma omp for nowait
2807                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2808                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2809                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2810                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2811                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2812                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2813                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
2814                            for (index_t i=0; i < numComp; ++i) {
2815                                o[INDEX2(i,numComp,0)] = (f_100[i] + f_101[i] + f_110[i] + f_111[i])/4;
2816                            } // end of component loop i
2817                        } // end of k1 loop
2818                    } // end of k2 loop
2819                } // end of face 1
2820                if (m_faceOffset[2] > -1) {
2821    #pragma omp for nowait
2822                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2823                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2824                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2825                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2826                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2827                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2828                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
2829                            for (index_t i=0; i < numComp; ++i) {
2830                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_100[i] + f_101[i])/4;
2831                            } // end of component loop i
2832                        } // end of k0 loop
2833                    } // end of k2 loop
2834                } // end of face 2
2835                if (m_faceOffset[3] > -1) {
2836    #pragma omp for nowait
2837                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2838                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2839                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2840                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2841                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2842                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2843                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
2844                            for (index_t i=0; i < numComp; ++i) {
2845                                o[INDEX2(i,numComp,0)] = (f_010[i] + f_011[i] + f_110[i] + f_111[i])/4;
2846                            } // end of component loop i
2847                        } // end of k0 loop
2848                    } // end of k2 loop
2849                } // end of face 3
2850                if (m_faceOffset[4] > -1) {
2851    #pragma omp for nowait
2852                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2853                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2854                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2855                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2856                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2857                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2858                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
2859                            for (index_t i=0; i < numComp; ++i) {
2860                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_010[i] + f_100[i] + f_110[i])/4;
2861                            } // end of component loop i
2862                        } // end of k0 loop
2863                    } // end of k1 loop
2864                } // end of face 4
2865                if (m_faceOffset[5] > -1) {
2866    #pragma omp for nowait
2867                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2868                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2869                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2870                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2871                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2872                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2873                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
2874                            for (index_t i=0; i < numComp; ++i) {
2875                                o[INDEX2(i,numComp,0)] = (f_001[i] + f_011[i] + f_101[i] + f_111[i])/4;
2876                            } // end of component loop i
2877                        } // end of k0 loop
2878                    } // end of k1 loop
2879                } // end of face 5
2880            } // end of parallel section
2881        } else {
2882            out.requireWrite();
2883            const double c0 = 0.044658198738520451079;
2884            const double c1 = 0.16666666666666666667;
2885            const double c2 = 0.62200846792814621559;
2886    #pragma omp parallel
2887            {
2888                vector<double> f_000(numComp);
2889                vector<double> f_001(numComp);
2890                vector<double> f_010(numComp);
2891                vector<double> f_011(numComp);
2892                vector<double> f_100(numComp);
2893                vector<double> f_101(numComp);
2894                vector<double> f_110(numComp);
2895                vector<double> f_111(numComp);
2896                if (m_faceOffset[0] > -1) {
2897    #pragma omp for nowait
2898                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2899                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2900                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2901                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2902                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2903                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2904                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
2905                            for (index_t i=0; i < numComp; ++i) {
2906                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);
2907                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);
2908                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);
2909                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);
2910                            } // end of component loop i
2911                        } // end of k1 loop
2912                    } // end of k2 loop
2913                } // end of face 0
2914                if (m_faceOffset[1] > -1) {
2915    #pragma omp for nowait
2916                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2917                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2918                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2919                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2920                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2921                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2922                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
2923                            for (index_t i=0; i < numComp; ++i) {
2924                                o[INDEX2(i,numComp,0)] = f_100[i]*c2 + f_111[i]*c0 + c1*(f_101[i] + f_110[i]);
2925                                o[INDEX2(i,numComp,1)] = f_101[i]*c0 + f_110[i]*c2 + c1*(f_100[i] + f_111[i]);
2926                                o[INDEX2(i,numComp,2)] = f_101[i]*c2 + f_110[i]*c0 + c1*(f_100[i] + f_111[i]);
2927                                o[INDEX2(i,numComp,3)] = f_100[i]*c0 + f_111[i]*c2 + c1*(f_101[i] + f_110[i]);
2928                            } // end of component loop i
2929                        } // end of k1 loop
2930                    } // end of k2 loop
2931                } // end of face 1
2932                if (m_faceOffset[2] > -1) {
2933    #pragma omp for nowait
2934                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2935                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2936                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2937                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2938                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2939                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2940                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
2941                            for (index_t i=0; i < numComp; ++i) {
2942                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_100[i]);
2943                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_101[i]);
2944                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_101[i]);
2945                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_100[i]);
2946                            } // end of component loop i
2947                        } // end of k0 loop
2948                    } // end of k2 loop
2949                } // end of face 2
2950                if (m_faceOffset[3] > -1) {
2951    #pragma omp for nowait
2952                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2953                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2954                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2955                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2956                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2957                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2958                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
2959                            for (index_t i=0; i < numComp; ++i) {
2960                                o[INDEX2(i,numComp,0)] = f_010[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_110[i]);
2961                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_111[i]);
2962                                o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_111[i]);
2963                                o[INDEX2(i,numComp,3)] = f_010[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_110[i]);
2964                            } // end of component loop i
2965                        } // end of k0 loop
2966                    } // end of k2 loop
2967                } // end of face 3
2968                if (m_faceOffset[4] > -1) {
2969    #pragma omp for nowait
2970                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2971                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2972                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2973                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2974                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2975                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2976                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
2977                            for (index_t i=0; i < numComp; ++i) {
2978                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_100[i]);
2979                                o[INDEX2(i,numComp,1)] = f_010[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_110[i]);
2980                                o[INDEX2(i,numComp,2)] = f_010[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_110[i]);
2981                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_100[i]);
2982                            } // end of component loop i
2983                        } // end of k0 loop
2984                    } // end of k1 loop
2985                } // end of face 4
2986                if (m_faceOffset[5] > -1) {
2987    #pragma omp for nowait
2988                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2989                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2990                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2991                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2992                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2993                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2994                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
2995                            for (index_t i=0; i < numComp; ++i) {
2996                                o[INDEX2(i,numComp,0)] = f_001[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_101[i]);
2997                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_111[i]);
2998                                o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_111[i]);
2999                                o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_101[i]);
3000                            } // end of component loop i
3001                        } // end of k0 loop
3002                    } // end of k1 loop
3003                } // end of face 5
3004            } // end of parallel section
3005        }
3006    }
3007    
3008    namespace
3009    {
3010        // Calculates a gaussian blur convolution matrix for 3D
3011        // See wiki article on the subject
3012        double* get3DGauss(unsigned radius, double sigma)
3013        {
3014            double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
3015            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
3016            double total=0;
3017            int r=static_cast<int>(radius);
3018            for (int z=-r;z<=r;++z)
3019            {
3020                for (int y=-r;y<=r;++y)
3021                {
3022                    for (int x=-r;x<=r;++x)
3023                    {        
3024                        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));
3025                        total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];    
3026                    }
3027                }
3028            }
3029            double invtotal=1/total;
3030            for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
3031            {
3032                arr[p]*=invtotal;
3033            }
3034            return arr;
3035        }
3036        
3037        // applies conv to source to get a point.
3038        // (xp, yp) are the coords in the source matrix not the destination matrix
3039        double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
3040        {
3041            size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
3042            size_t sbase=bx+by*width+bz*width*height;
3043            double result=0;
3044            for (int z=0;z<2*radius+1;++z)
3045            {
3046                for (int y=0;y<2*radius+1;++y)
3047                {    
3048                    for (int x=0;x<2*radius+1;++x)
3049                    {
3050                        result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
3051                    }
3052              }              }
3053          }          }
3054            // use this line for "pass-though" (return the centre point value)
3055    //      return source[sbase+(radius)+(radius)*width+(radius)*width*height];
3056            return result;      
3057      }      }
3058    }
3059    
3060      // elements  /* This is a wrapper for filtered (and non-filtered) randoms
3061      m_elementId.resize(getNumElements());   * For detailed doco see randomFillWorker
3062  #pragma omp parallel for  */
3063      for (dim_t k=0; k<getNumElements(); k++) {  escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
3064          m_elementId[k]=k;         const escript::FunctionSpace& what,
3065           long seed, const boost::python::tuple& filter) const
3066    {
3067        int numvals=escript::DataTypes::noValues(shape);
3068        if (len(filter)>0 && (numvals!=1))
3069        {
3070            throw RipleyException("Ripley only supports filters for scalar data.");
3071        }
3072        escript::Data res=randomFillWorker(shape, seed, filter);
3073        if (res.getFunctionSpace()!=what)
3074        {
3075            escript::Data r=escript::Data(res, what);
3076            return r;
3077      }      }
3078        return res;
3079    }
3080    
3081      // face elements  /* This routine produces a Data object filled with smoothed random data.
3082      m_faceId.resize(getNumFaceElements());  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
3083  #pragma omp parallel for  A parameter radius  gives the size of the stencil used for the smoothing.
3084      for (dim_t k=0; k<getNumFaceElements(); k++) {  A point on the left hand edge for example, will still require `radius` extra points to the left
3085          m_faceId[k]=k;  in order to complete the stencil.
3086    
3087    All local calculation is done on an array called `src`, with
3088    dimensions = ext[0] * ext[1] *ext[2].
3089    Where ext[i]= internal[i]+2*radius.
3090    
3091    Now for MPI there is overlap to deal with. We need to share both the overlapping
3092    values themselves but also the external region.
3093    
3094    In a hypothetical 1-D case:
3095    
3096    
3097    1234567
3098    would be split into two ranks thus:
3099    123(4)  (4)567     [4 being a shared element]
3100    
3101    If the radius is 2.   There will be padding elements on the outside:
3102    
3103    pp123(4)  (4)567pp
3104    
3105    To ensure that 4 can be correctly computed on both ranks, values from the other rank
3106    need to be known.
3107    
3108    pp123(4)56   23(4)567pp
3109    
3110    Now in our case, we wout set all the values 23456 on the left rank and send them to the
3111    right hand rank.
3112    
3113    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
3114    inset=2*radius+1    
3115    This is to ensure that values at distance `radius` from the shared/overlapped element
3116    that ripley has.
3117    */
3118    escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
3119    {
3120        unsigned int radius=0;  // these are only used by gaussian
3121        double sigma=0.5;
3122        
3123        unsigned int numvals=escript::DataTypes::noValues(shape);
3124        
3125        if (len(filter)==0)
3126        {
3127        // nothing special required here yet
3128        }
3129        else if (len(filter)==3)
3130        {
3131            boost::python::extract<string> ex(filter[0]);
3132            if (!ex.check() || (ex()!="gaussian"))
3133            {
3134                throw RipleyException("Unsupported random filter for Brick.");
3135            }
3136            boost::python::extract<unsigned int> ex1(filter[1]);
3137            if (!ex1.check())
3138            {
3139                throw RipleyException("Radius of gaussian filter must be a positive integer.");
3140            }
3141            radius=ex1();
3142            sigma=0.5;
3143            boost::python::extract<double> ex2(filter[2]);
3144            if (!ex2.check() || (sigma=ex2())<=0)
3145            {
3146                throw RipleyException("Sigma must be a positive floating point number.");
3147            }            
3148        }
3149        else
3150        {
3151            throw RipleyException("Unsupported random filter");
3152        }
3153    
3154        // number of points in the internal region
3155        // that is, the ones we need smoothed versions of
3156        const dim_t internal[3] = { m_NN[0], m_NN[1], m_NN[2] };
3157        size_t ext[3];
3158        ext[0]=(size_t)internal[0]+2*radius;  // includes points we need as input
3159        ext[1]=(size_t)internal[1]+2*radius;  // for smoothing
3160        ext[2]=(size_t)internal[2]+2*radius;  // for smoothing
3161        
3162        // now we check to see if the radius is acceptable
3163        // That is, would not cross multiple ranks in MPI
3164    
3165        if (2*radius>=internal[0]-4)
3166        {
3167            throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3168        }
3169        if (2*radius>=internal[1]-4)
3170        {
3171            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3172        }
3173        if (2*radius>=internal[2]-4)
3174        {
3175            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3176        }    
3177      
3178        double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3179        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3180        
3181    #ifdef ESYS_MPI
3182      
3183        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3184        {
3185        // since the dimensions are equal for all ranks, this exception
3186        // will be thrown on all ranks
3187        throw RipleyException("Random Data in Ripley requires at least five elements per side per rank.");
3188    
3189        }
3190        dim_t X=m_mpiInfo->rank%m_NX[0];
3191        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3192        dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3193    #endif    
3194    
3195    /*