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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26