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