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