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