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

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

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

branches/ripleygmg_from_3668/ripley/src/Brick.cpp revision 3691 by caltinay, Wed Nov 23 23:07:37 2011 UTC 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() };
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]=(m_l0*(i0+m_offset0))/m_gNE0;              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]=(m_l1*(i1+m_offset1))/m_gNE1;              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]=(m_l2*(i2+m_offset2))/m_gNE2;              coords[2][i2]=getLocalCoordinate(i2, 2);
707          }          }
708      }      }
709      int dims[] = { m_N0, m_N1, m_N2 };      int* dims = const_cast<int*>(getNumNodesPerDim());
710    
711        // write mesh
712      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,      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, 3, NULL, 0,      // write node ids
716              DB_INT, DB_NODECENT, NULL);      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 3,
717                NULL, 0, DB_INT, DB_NODECENT, NULL);
718    
719        // write element ids
720        dims = const_cast<int*>(getNumElementsPerDim());
721        DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
722                dims, 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
723    
724        // rank 0 writes multimesh and multivar
725      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
726          vector<string> tempstrings;          vector<string> tempstrings;
727          vector<char*> names;          vector<char*> names;
# Line 193  void Brick::dump(const string& fileName) Line 746  void Brick::dump(const string& fileName)
746          types.assign(m_mpiInfo->size, DB_QUADVAR);          types.assign(m_mpiInfo->size, DB_QUADVAR);
747          DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],          DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
748                 &types[0], NULL);                 &types[0], NULL);
749            tempstrings.clear();
750            names.clear();
751            for (dim_t i=0; i<m_mpiInfo->size; i++) {
752                stringstream path;
753                path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
754                tempstrings.push_back(path.str());
755                names.push_back((char*)tempstrings.back().c_str());
756            }
757            DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
758                   &types[0], NULL);
759      }      }
760    
761      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 205  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    
775  const int* Brick::borrowSampleReferenceIDs(int fsType) const  const int* Brick::borrowSampleReferenceIDs(int fsType) const
776  {  {
777      if (fsType == Nodes)      switch (fsType) {
778          return &m_nodeId[0];          case Nodes:
779            case ReducedNodes: //FIXME: reduced
780                return &m_nodeId[0];
781            case DegreesOfFreedom:
782            case ReducedDegreesOfFreedom: //FIXME: reduced
783                return &m_dofId[0];
784            case Elements:
785            case ReducedElements:
786                return &m_elementId[0];
787            case FaceElements:
788            case ReducedFaceElements:
789                return &m_faceId[0];
790            case Points:
791                return &m_diracPointNodeIDs[0];
792            default:
793                break;
794        }
795    
796      throw RipleyException("borrowSampleReferenceIDs() only implemented for Nodes");      stringstream msg;
797        msg << "borrowSampleReferenceIDs: invalid function space type "<<fsType;
798        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=getNumNodes()*m_mpiInfo->rank;  
806          const index_t myLast=getNumNodes()*(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::interpolateOnDomain(escript::Data& target,  void Brick::setToNormal(escript::Data& out) const
                                 const escript::Data& in) const  
855  {  {
856      const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
857      const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));          out.requireWrite();
858      if (inDomain != *this)  #pragma omp parallel
859          throw RipleyException("Illegal domain of interpolant");          {
860      if (targetDomain != *this)              if (m_faceOffset[0] > -1) {
861          throw RipleyException("Illegal domain of interpolation target");  #pragma omp for nowait
862                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
863                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
864                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
865                            // set vector at four quadrature points
866                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
867                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
868                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
869                            *o++ = -1.; *o++ = 0.; *o = 0.;
870                        }
871                    }
872                }
873    
874                if (m_faceOffset[1] > -1) {
875    #pragma omp for nowait
876                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
877                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
878                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
879                            // set vector at four quadrature points
880                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
881                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
882                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
883                            *o++ = 1.; *o++ = 0.; *o = 0.;
884                        }
885                    }
886                }
887    
888                if (m_faceOffset[2] > -1) {
889    #pragma omp for nowait
890                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
891                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
892                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
893                            // set vector at four quadrature points
894                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
895                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
896                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
897                            *o++ = 0.; *o++ = -1.; *o = 0.;
898                        }
899                    }
900                }
901    
902                if (m_faceOffset[3] > -1) {
903    #pragma omp for nowait
904                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
905                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
906                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
907                            // set vector at four quadrature points
908                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
909                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
910                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
911                            *o++ = 0.; *o++ = 1.; *o = 0.;
912                        }
913                    }
914                }
915    
916                if (m_faceOffset[4] > -1) {
917    #pragma omp for nowait
918                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
919                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
920                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
921                            // set vector at four quadrature points
922                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
923                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
924                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
925                            *o++ = 0.; *o++ = 0.; *o = -1.;
926                        }
927                    }
928                }
929    
930                if (m_faceOffset[5] > -1) {
931    #pragma omp for nowait
932                    for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
933                        for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {
934                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
935                            // set vector at four quadrature points
936                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
937                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
938                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
939                            *o++ = 0.; *o++ = 0.; *o = 1.;
940                        }
941                    }
942                }
943            } // end of parallel section
944        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
945            out.requireWrite();
946    #pragma omp parallel
947            {
948                if (m_faceOffset[0] > -1) {
949    #pragma omp for nowait
950                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
951                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
952                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
953                            *o++ = -1.;
954                            *o++ = 0.;
955                            *o = 0.;
956                        }
957                    }
958                }
959    
960                if (m_faceOffset[1] > -1) {
961    #pragma omp for nowait
962                    for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {
963                        for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {
964                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
965                            *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      throw RipleyException("interpolateOnDomain() not implemented");              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 {
1022            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                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      throw RipleyException("getPattern() not implemented");              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 262  void Brick::Print_Mesh_Info(const bool f Line 1128  void Brick::Print_Mesh_Info(const bool f
1128          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
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                  << "  " << (m_l0*(i%m_N0+m_offset0))/m_gNE0                  << "  " << getLocalCoordinate(i%m_NN[0], 0)
1132                  << "  " << (m_l1*(i%(m_N0*m_N1)/m_N0+m_offset1))/m_gNE1                  << "  " << getLocalCoordinate(i%(m_NN[0]*m_NN[1])/m_NN[0], 1)
1133                  << "  " << (m_l2*(i/(m_N0*m_N1)+m_offset2))/m_gNE2 << endl;                  << "  " << getLocalCoordinate(i/(m_NN[0]*m_NN[1]), 2) << endl;
1134          }          }
1135      }      }
1136  }  }
1137    
 //protected  
 dim_t Brick::getNumFaceElements() const  
 {  
     dim_t n=0;  
     //left  
     if (m_offset0==0)  
         n+=m_NE1*m_NE2;  
     //right  
     if (m_mpiInfo->rank%m_NX==m_NX-1)  
         n+=m_NE1*m_NE2;  
     //bottom  
     if (m_offset1==0)  
         n+=m_NE0*m_NE2;  
     //top  
     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)  
         n+=m_NE0*m_NE2;  
     //front  
     if (m_offset2==0)  
         n+=m_NE0*m_NE1;  
     //back  
     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)  
         n+=m_NE0*m_NE1;  
   
     return n;  
 }  
1138    
1139  //protected  //protected
1140  void Brick::assembleCoordinates(escript::Data& arg) const  void Brick::assembleCoordinates(escript::Data& arg) const
# Line 307  void Brick::assembleCoordinates(escript: Line 1148  void Brick::assembleCoordinates(escript:
1148    
1149      arg.requireWrite();      arg.requireWrite();
1150  #pragma omp parallel for  #pragma omp parallel for
1151      for (dim_t i2 = 0; i2 < m_N2; i2++) {      for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {
1152          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {
1153              for (dim_t i0 = 0; i0 < m_N0; i0++) {              for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {
1154                  double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);                  double* point = arg.getSampleDataRW(i0+m_NN[0]*i1+m_NN[0]*m_NN[1]*i2);
1155                  point[0] = (m_l0*(i0+m_offset0))/m_gNE0;                  point[0] = getLocalCoordinate(i0, 0);
1156                  point[1] = (m_l1*(i1+m_offset1))/m_gNE1;                  point[1] = getLocalCoordinate(i1, 1);
1157                  point[2] = (m_l2*(i2+m_offset2))/m_gNE2;                  point[2] = getLocalCoordinate(i2, 2);
1158              }              }
1159          }          }
1160      }      }
1161  }  }
1162    
1163    //protected
1164    void Brick::assembleGradient(escript::Data& out, const escript::Data& in) const
1165    {
1166        const dim_t numComp = in.getDataPointSize();
1167        const double C0 = .044658198738520451079;
1168        const double C1 = .16666666666666666667;
1169        const double C2 = .21132486540518711775;
1170        const double C3 = .25;
1171        const double C4 = .5;
1172        const double C5 = .62200846792814621559;
1173        const double C6 = .78867513459481288225;
1174    
1175        if (out.getFunctionSpace().getTypeCode() == Elements) {
1176            out.requireWrite();
1177    #pragma omp parallel
1178            {
1179                vector<double> f_000(numComp);
1180                vector<double> f_001(numComp);
1181                vector<double> f_010(numComp);
1182                vector<double> f_011(numComp);
1183                vector<double> f_100(numComp);
1184                vector<double> f_101(numComp);
1185                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    //protected
1635    void Brick::assembleIntegrate(vector<double>& integrals, const escript::Data& arg) const
1636    {
1637        const dim_t numComp = arg.getDataPointSize();
1638        const index_t left = (m_offset[0]==0 ? 0 : 1);
1639        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1640        const index_t front = (m_offset[2]==0 ? 0 : 1);
1641        const int fs = arg.getFunctionSpace().getTypeCode();
1642        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        } 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
1889    dim_t Brick::insertNeighbourNodes(IndexVector& index, index_t node) const
1890    {
1891        const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1892        const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1893        const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
1894        const int x=node%nDOF0;
1895        const int y=node%(nDOF0*nDOF1)/nDOF0;
1896        const int z=node/(nDOF0*nDOF1);
1897        int num=0;
1898        // loop through potential neighbours and add to index if positions are
1899        // within bounds
1900        for (int i2=-1; i2<2; i2++) {
1901            for (int i1=-1; i1<2; i1++) {
1902                for (int i0=-1; i0<2; i0++) {
1903                    // skip node itself
1904                    if (i0==0 && i1==0 && i2==0)
1905                        continue;
1906                    // location of neighbour node
1907                    const int nx=x+i0;
1908                    const int ny=y+i1;
1909                    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 num;
1920    }
1921    
1922    //protected
1923    void Brick::nodesToDOF(escript::Data& out, const escript::Data& in) const
1924    {
1925        const dim_t numComp = in.getDataPointSize();
1926        out.requireWrite();
1927    
1928        const index_t left = (m_offset[0]==0 ? 0 : 1);
1929        const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1930        const index_t front = (m_offset[2]==0 ? 0 : 1);
1931        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
1935        for (index_t i=0; i<nDOF2; i++) {
1936            for (index_t j=0; j<nDOF1; j++) {
1937                for (index_t k=0; k<nDOF0; k++) {
1938                    const index_t n=k+left+(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];
1939                    const double* src=in.getSampleDataRO(n);
1940                    copy(src, src+numComp, out.getSampleDataRW(k+j*nDOF0+i*nDOF0*nDOF1));
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      const index_t firstId = getNumNodes()*m_mpiInfo->rank;      // degrees of freedom are numbered from left to right, bottom to top, front
1973      const index_t diff0 = m_N0*(m_N1*m_N2-1)+1;      // to back in each rank, continuing on the next rank (ranks also go
1974      const index_t diff1 = m_N0*m_N1*(m_N2*m_NX-1)+m_N0;      // left-right, bottom-top, front-back).
1975      const index_t diff2 = m_N0*m_N1*m_N2*m_NX*m_NY-m_N0*m_N1*(m_N2-1);      // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1976        // helps when writing out data rank after rank.
1977    
1978        // build node distribution vector first.
1979        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1980        // constant for all ranks in this implementation
1981        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1982        const dim_t numDOF=getNumDOF();
1983        for (dim_t k=1; k<m_mpiInfo->size; k++) {
1984            m_nodeDistribution[k]=k*numDOF;
1985        }
1986        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        // populate face element counts
1992        //left
1993        if (m_offset[0]==0)
1994            m_faceCount[0]=m_NE[1]*m_NE[2];
1995        else
1996            m_faceCount[0]=0;
1997        //right
1998        if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1999            m_faceCount[1]=m_NE[1]*m_NE[2];
2000        else
2001            m_faceCount[1]=0;
2002        //bottom
2003        if (m_offset[1]==0)
2004            m_faceCount[2]=m_NE[0]*m_NE[2];
2005        else
2006            m_faceCount[2]=0;
2007        //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            // populate the rest of the nodes (shared with other ranks)
2088            if (m_faceCount[0]==0) { // left plane
2089    #pragma omp for nowait
2090                for (dim_t i=0; i<nDOF2; i++) {
2091                    for (dim_t j=0; j<nDOF1; j++) {
2092                        const index_t nodeIdx=(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];
2093                        const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;
2094                        m_nodeId[nodeIdx]
2095                            = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
2096                    }
2097                }
2098            }
2099            if (m_faceCount[1]==0) { // right plane
2100    #pragma omp for nowait
2101                for (dim_t i=0; i<nDOF2; i++) {
2102                    for (dim_t j=0; j<nDOF1; j++) {
2103                        const index_t nodeIdx=(j+bottom+1)*m_NN[0]-1+(i+front)*m_NN[0]*m_NN[1];
2104                        const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;
2105                        m_nodeId[nodeIdx]
2106                            = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
2107                    }
2108                }
2109            }
2110            if (m_faceCount[2]==0) { // bottom plane
2111    #pragma omp for nowait
2112                for (dim_t i=0; i<nDOF2; i++) {
2113                    for (dim_t k=0; k<nDOF0; k++) {
2114                        const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1];
2115                        const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;
2116                        m_nodeId[nodeIdx]
2117                            = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
2118                    }
2119                }
2120            }
2121            if (m_faceCount[3]==0) { // top plane
2122    #pragma omp for nowait
2123                for (dim_t i=0; i<nDOF2; i++) {
2124                    for (dim_t k=0; k<nDOF0; k++) {
2125                        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 dofId=k+i*nDOF0*nDOF1;
2127                        m_nodeId[nodeIdx]
2128                            = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
2129                    }
2130                }
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        setTagMap("left", LEFT);
2196        setTagMap("right", RIGHT);
2197        setTagMap("bottom", BOTTOM);
2198        setTagMap("top", TOP);
2199        setTagMap("front", FRONT);
2200        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 k=0; k<getNumNodes(); k++) {      for (index_t i=front; i<front+nDOF2; i++) {
2219          index_t id = firstId+k;          for (index_t j=bottom; j<bottom+nDOF1; j++) {
2220          if (m_offset0 > 0 && k%m_N0==0)              for (index_t k=left; k<left+nDOF0; k++) {
2221              id -= diff0;                  m_dofMap[i*m_NN[0]*m_NN[1]+j*m_NN[0]+k]=(i-front)*nDOF0*nDOF1+(j-bottom)*nDOF0+k-left;
2222          if (m_offset1 > 0 && k%(m_N0*m_N1)<m_N0)              }
2223              id -= diff1;          }
2224          if (m_offset2 > 0 && k/(m_N0*m_N1)==0)      }
2225              id -= diff2;  
2226          m_nodeId[k]=id;      // build list of shared components and neighbours by looping through
2227        // all potential neighbouring ranks and checking if positions are
2228        // within bounds
2229        const dim_t numDOF=nDOF0*nDOF1*nDOF2;
2230        vector<IndexVector> colIndices(numDOF); // for the couple blocks
2231        RankVector neighbour;
2232        IndexVector offsetInShared(1,0);
2233        IndexVector sendShared, recvShared;
2234        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        // create connector
2425        Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
2426                numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
2427                &offsetInShared[0], 1, 0, m_mpiInfo);
2428        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        /*
2477        cout << "--- main_pattern ---" << endl;
2478        cout << "M=" << mainPattern->numOutput << ", N=" << mainPattern->numInput << endl;
2479        for (size_t i=0; i<mainPattern->numOutput+1; i++) {
2480            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        /*
2488        cout << "--- colCouple_pattern ---" << endl;
2489        cout << "M=" << colPattern->numOutput << ", N=" << colPattern->numInput << endl;
2490        for (size_t i=0; i<colPattern->numOutput+1; i++) {
2491            cout << "ptr[" << i << "]=" << colPattern->ptr[i] << endl;
2492        }
2493        for (size_t i=0; i<colPattern->ptr[colPattern->numOutput]; i++) {
2494            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
2544    void Brick::interpolateNodesOnElements(escript::Data& out,
2545                                           const escript::Data& in,
2546                                           bool reduced) const
2547    {
2548        const dim_t numComp = in.getDataPointSize();
2549        if (reduced) {
2550            out.requireWrite();
2551    #pragma omp parallel
2552            {
2553                vector<double> f_000(numComp);
2554                vector<double> f_001(numComp);
2555                vector<double> f_010(numComp);
2556                vector<double> f_011(numComp);
2557                vector<double> f_100(numComp);
2558                vector<double> f_101(numComp);
2559                vector<double> f_110(numComp);
2560                vector<double> f_111(numComp);
2561    #pragma omp for
2562                for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2563                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2564                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2565                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2566                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2567                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2568                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2569                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2570                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2571                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2572                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2573                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));
2574                            for (index_t i=0; i < numComp; ++i) {
2575                                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
2577                        } // end of k0 loop
2578                    } // end of k1 loop
2579                } // end of k2 loop
2580            } // 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
2628    void Brick::interpolateNodesOnFaces(escript::Data& out, const escript::Data& in,
2629                                        bool reduced) const
2630    {
2631        const dim_t numComp = in.getDataPointSize();
2632        if (reduced) {
2633            out.requireWrite();
2634    #pragma omp parallel
2635            {
2636                vector<double> f_000(numComp);
2637                vector<double> f_001(numComp);
2638                vector<double> f_010(numComp);
2639                vector<double> f_011(numComp);
2640                vector<double> f_100(numComp);
2641                vector<double> f_101(numComp);
2642                vector<double> f_110(numComp);
2643                vector<double> f_111(numComp);
2644                if (m_faceOffset[0] > -1) {
2645    #pragma omp for nowait
2646                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2647                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2648                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2649                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2650                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2651                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2652                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
2653                            for (index_t i=0; i < numComp; ++i) {
2654                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_010[i] + f_011[i])/4;
2655                            } // end of component loop i
2656                        } // end of k1 loop
2657                    } // end of k2 loop
2658                } // end of face 0
2659                if (m_faceOffset[1] > -1) {
2660    #pragma omp for nowait
2661                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2662                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2663                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2664                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2665                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2666                            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]+INDEX2(k1,k2,m_NE[1]));
2668                            for (index_t i=0; i < numComp; ++i) {
2669                                o[INDEX2(i,numComp,0)] = (f_100[i] + f_101[i] + f_110[i] + f_111[i])/4;
2670                            } // end of component loop i
2671                        } // end of k1 loop
2672                    } // end of k2 loop
2673                } // end of face 1
2674                if (m_faceOffset[2] > -1) {
2675    #pragma omp for nowait
2676                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2677                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2678                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2679                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2680                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2681                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2682                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
2683                            for (index_t i=0; i < numComp; ++i) {
2684                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_001[i] + f_100[i] + f_101[i])/4;
2685                            } // end of component loop i
2686                        } // end of k0 loop
2687                    } // end of k2 loop
2688                } // end of face 2
2689                if (m_faceOffset[3] > -1) {
2690    #pragma omp for nowait
2691                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2692                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2693                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2694                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2695                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2696                            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                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
2698                            for (index_t i=0; i < numComp; ++i) {
2699                                o[INDEX2(i,numComp,0)] = (f_010[i] + f_011[i] + f_110[i] + f_111[i])/4;
2700                            } // end of component loop i
2701                        } // end of k0 loop
2702                    } // end of k2 loop
2703                } // end of face 3
2704                if (m_faceOffset[4] > -1) {
2705    #pragma omp for nowait
2706                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2707                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2708                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2709                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2710                            memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2711                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
2712                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
2713                            for (index_t i=0; i < numComp; ++i) {
2714                                o[INDEX2(i,numComp,0)] = (f_000[i] + f_010[i] + f_100[i] + f_110[i])/4;
2715                            } // end of component loop i
2716                        } // end of k0 loop
2717                    } // end of k1 loop
2718                } // end of face 4
2719                if (m_faceOffset[5] > -1) {
2720    #pragma omp for nowait
2721                    for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2722                        for (index_t k0=0; k0 < m_NE[0]; ++k0) {
2723                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2724                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2725                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2726                            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                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
2728                            for (index_t i=0; i < numComp; ++i) {
2729                                o[INDEX2(i,numComp,0)] = (f_001[i] + f_011[i] + f_101[i] + f_111[i])/4;
2730                            } // end of component loop i
2731                        } // end of k0 loop
2732                    } // end of k1 loop
2733                } // end of face 5
2734            } // end of parallel section
2735        } else {
2736            out.requireWrite();
2737            const double c0 = 0.044658198738520451079;
2738            const double c1 = 0.16666666666666666667;
2739            const double c2 = 0.62200846792814621559;
2740    #pragma omp parallel
2741            {
2742                vector<double> f_000(numComp);
2743                vector<double> f_001(numComp);
2744                vector<double> f_010(numComp);
2745                vector<double> f_011(numComp);
2746                vector<double> f_100(numComp);
2747                vector<double> f_101(numComp);
2748                vector<double> f_110(numComp);
2749                vector<double> f_111(numComp);
2750                if (m_faceOffset[0] > -1) {
2751    #pragma omp for nowait
2752                    for (index_t k2=0; k2 < m_NE[2]; ++k2) {
2753                        for (index_t k1=0; k1 < m_NE[1]; ++k1) {
2754                            memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2755                            memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2756                            memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
2757                            memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
2758                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
2759                            for (index_t i=0; i < numComp; ++i) {
2760                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);
2761                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);
2762                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);
2763                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);
2764                            } // end of component loop i
2765                        } // 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    
3092        BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3093        size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3094            // there is an overlap of two points both of which need to have "radius" points on either side.
3095        
3096        size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3097        size_t ymidlen=ext[1]-2*inset;  
3098        size_t zmidlen=ext[2]-2*inset;
3099        
3100        Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);    
3101        
3102        MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3103        MPI_Status stats[50];
3104        short rused=0;
3105        
3106        messvec incoms;
3107        messvec outcoms;
3108        
3109        grid.generateInNeighbours(X, Y, Z ,incoms);
3110        grid.generateOutNeighbours(X, Y, Z, outcoms);
3111        
3112        
3113        block.copyAllToBuffer(src);
3114        
3115        
3116        int comserr=0;    
3117        for (size_t i=0;i<incoms.size();++i)
3118        {
3119        message& m=incoms[i];
3120        comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3121        block.setUsed(m.destbuffid);
3122        }
3123    
3124        for (size_t i=0;i<outcoms.size();++i)
3125        {
3126        message& m=outcoms[i];
3127        comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3128        }    
3129        
3130        if (!comserr)
3131        {    
3132            comserr=MPI_Waitall(rused, reqs, stats);
3133        }
3134    
3135        if (comserr)
3136        {
3137        // Yes this is throwing an exception as a result of an MPI error.
3138        // and no we don't inform the other ranks that we are doing this.
3139        // however, we have no reason to believe coms work at this point anyway
3140            throw RipleyException("Error in coms for randomFill");      
3141        }
3142        
3143        block.copyUsedFromBuffer(src);
3144        
3145        
3146        
3147    
3148    #endif    
3149        
3150    /*    
3151    cout << "Pattern (post transfer):\n";    
3152    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3153    {
3154        cout << src[i++] << " ";
3155        if (i%ext[0]==0)
3156          cout << "\n";
3157        if (i%(ext[0]*ext[1])==0)
3158          cout << "\n";
3159    }   */
3160        
3161        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3162        {
3163          
3164        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3165        escript::Data resdat(0, shape, fs , true);
3166        // don't need to check for exwrite because we just made it
3167        escript::DataVector& dv=resdat.getExpandedVectorReference();
3168        
3169        
3170        // now we need to copy values over
3171        for (size_t z=0;z<(internal[2]);++z)
3172        {
3173            for (size_t y=0;y<(internal[1]);++y)    
3174            {
3175            for (size_t x=0;x<(internal[0]);++x)
3176            {
3177                for (unsigned int i=0;i<numvals;++i)
3178                {
3179                    dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3180                }
3181            }
3182            }
3183        }  
3184        delete[] src;
3185        return resdat;      
3186        }
3187        else        // filter enabled  
3188        {
3189        
3190        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3191        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3192        // don't need to check for exwrite because we just made it
3193        escript::DataVector& dv=resdat.getExpandedVectorReference();
3194        double* convolution=get3DGauss(radius, sigma);
3195    
3196    // cout << "Convolution matrix\n";
3197    // size_t di=(radius*2+1);
3198    // double ctot=0;
3199    // for (int i=0;i<di*di*di;++i)
3200    // {
3201    //     cout << convolution[i] << " ";
3202    //     ctot+=convolution[i];
3203    //     if ((i+1)%di==0)
3204    //     {
3205    //  cout << "\n";
3206    //     }
3207    //     if ((i+1)%(di*di)==0)
3208    //     {
3209    //  cout << "\n";
3210    //     }
3211    // }
3212    //
3213    // cout << "Sum of matrix is = " << ctot << endl;
3214    
3215        for (size_t z=0;z<(internal[2]);++z)
3216        {
3217            for (size_t y=0;y<(internal[1]);++y)    
3218            {
3219            for (size_t x=0;x<(internal[0]);++x)
3220            {    
3221                dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3222                
3223            }
3224            }
3225        }
3226        
3227    // cout << "\nResulting matrix:\n";
3228    //     for (size_t z=0;z<(internal[2]);++z)
3229    //     {
3230    //  for (size_t y=0;y<(internal[1]);++y)    
3231    //  {
3232    //      for (size_t x=0;x<(internal[0]);++x)
3233    //      {    
3234    //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";
3235    //      }
3236    //      cout << endl;
3237    //  }
3238    //  cout << endl;
3239    //     }
3240    
3241        
3242        delete[] convolution;
3243        delete[] src;
3244        return resdat;
3245        
3246        }
3247    }
3248    
3249    
3250    
3251    
3252    
3253    
3254    int Brick::findNode(const double *coords) const {
3255        const int NOT_MINE = -1;
3256        //is the found element even owned by this rank
3257        // (inside owned or shared elements but will map to an owned element)
3258        for (int dim = 0; dim < m_numDim; dim++) {
3259            double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3260                    - m_dx[dim]/2.; //allows for point outside mapping onto node
3261            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3262                    + m_dx[dim]/2.;
3263            if (min > coords[dim] || max < coords[dim]) {
3264                return NOT_MINE;
3265            }
3266        }
3267        // get distance from origin
3268        double x = coords[0] - m_origin[0];
3269        double y = coords[1] - m_origin[1];
3270        double z = coords[2] - m_origin[2];
3271        // distance in elements
3272        int ex = (int) floor(x / m_dx[0]);
3273        int ey = (int) floor(y / m_dx[1]);
3274        int ez = (int) floor(z / m_dx[2]);
3275        // set the min distance high enough to be outside the element plus a bit
3276        int closest = NOT_MINE;
3277        double minDist = 1;
3278        for (int dim = 0; dim < m_numDim; dim++) {
3279            minDist += m_dx[dim]*m_dx[dim];
3280        }
3281        //find the closest node
3282        for (int dx = 0; dx < 1; dx++) {
3283            double xdist = x - (ex + dx)*m_dx[0];
3284            for (int dy = 0; dy < 1; dy++) {
3285                double ydist = y - (ey + dy)*m_dx[1];
3286                for (int dz = 0; dz < 1; dz++) {
3287                    double zdist = z - (ez + dz)*m_dx[2];
3288                    double total = xdist*xdist + ydist*ydist + zdist*zdist;
3289                    if (total < minDist) {
3290                        closest = INDEX3(ex+dy-m_offset[0], ey+dy-m_offset[1],
3291                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3292                        minDist = total;