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