/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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