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

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

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

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

Legend:
Removed from v.3691  
changed lines
  Added in v.4650

  ViewVC Help
Powered by ViewVC 1.1.26