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