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