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

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

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

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

Legend:
Removed from v.3702  
changed lines
  Added in v.4657

  ViewVC Help
Powered by ViewVC 1.1.26