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