/[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

revision 3796 by gross, Thu Feb 2 06:26:15 2012 UTC revision 4174 by caltinay, Wed Jan 30 03:21:27 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 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" {  extern "C" {
18  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
19  }  }
20    
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24    
25  #if USE_SILO  #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 32  namespace ripley { Line 38  namespace ripley {
38  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
39               double x1, double y1, double z1, int d0, int d1, int d2) :               double x1, double y1, double z1, int d0, int d1, int d2) :
40      RipleyDomain(3),      RipleyDomain(3),
     m_gNE0(n0),  
     m_gNE1(n1),  
     m_gNE2(n2),  
41      m_x0(x0),      m_x0(x0),
42      m_y0(y0),      m_y0(y0),
43      m_z0(z0),      m_z0(z0),
44      m_l0(x1-x0),      m_l0(x1-x0),
45      m_l1(y1-y0),      m_l1(y1-y0),
46      m_l2(z1-z0),      m_l2(z1-z0)
     m_NX(d0),  
     m_NY(d1),  
     m_NZ(d2)  
47  {  {
48        // ignore subdivision parameters for serial run
49        if (m_mpiInfo->size == 1) {
50            d0=1;
51            d1=1;
52            d2=1;
53        }
54    
55        bool warn=false;
56        // if number of subdivisions is non-positive, try to subdivide by the same
57        // ratio as the number of elements
58        if (d0<=0 && d1<=0 && d2<=0) {
59            warn=true;
60            d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));
61            d0=max(1, d0);
62            d1=max(1, (int)(d0*n1/(float)n0));
63            d2=m_mpiInfo->size/(d0*d1);
64            if (d0*d1*d2 != m_mpiInfo->size) {
65                // ratios not the same so leave "smallest" side undivided and try
66                // dividing 2 sides only
67                if (n0>=n1) {
68                    if (n1>=n2) {
69                        d0=d1=0;
70                        d2=1;
71                    } else {
72                        d0=d2=0;
73                        d1=1;
74                    }
75                } else {
76                    if (n0>=n2) {
77                        d0=d1=0;
78                        d2=1;
79                    } else {
80                        d0=1;
81                        d1=d2=0;
82                    }
83                }
84            }
85        }
86        if (d0<=0 && d1<=0) {
87            warn=true;
88            d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
89            d1=m_mpiInfo->size/d0;
90            if (d0*d1*d2 != m_mpiInfo->size) {
91                // ratios not the same so subdivide side with more elements only
92                if (n0>n1) {
93                    d0=0;
94                    d1=1;
95                } else {
96                    d0=1;
97                    d1=0;
98                }
99            }
100        } else if (d0<=0 && d2<=0) {
101            warn=true;
102            d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));
103            d2=m_mpiInfo->size/d0;
104            if (d0*d1*d2 != m_mpiInfo->size) {
105                // ratios not the same so subdivide side with more elements only
106                if (n0>n2) {
107                    d0=0;
108                    d2=1;
109                } else {
110                    d0=1;
111                    d2=0;
112                }
113            }
114        } else if (d1<=0 && d2<=0) {
115            warn=true;
116            d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));
117            d2=m_mpiInfo->size/d1;
118            if (d0*d1*d2 != m_mpiInfo->size) {
119                // ratios not the same so subdivide side with more elements only
120                if (n1>n2) {
121                    d1=0;
122                    d2=1;
123                } else {
124                    d1=1;
125                    d2=0;
126                }
127            }
128        }
129        if (d0<=0) {
130            // d1,d2 are preset, determine d0
131            d0=m_mpiInfo->size/(d1*d2);
132        } else if (d1<=0) {
133            // d0,d2 are preset, determine d1
134            d1=m_mpiInfo->size/(d0*d2);
135        } else if (d2<=0) {
136            // d0,d1 are preset, determine d2
137            d2=m_mpiInfo->size/(d0*d1);
138        }
139    
140        m_NX=d0;
141        m_NY=d1;
142        m_NZ=d2;
143    
144      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
145      // among number of ranks      // among number of ranks
146      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
147          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
148    
149      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0 || (n2+1)%m_NZ > 0)      if (warn) {
150          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
151                << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
152        }
153    
154        if ((n0+1)%m_NX > 0) {
155            double Dx=m_l0/n0;
156            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
157            m_l0=Dx*n0;
158            cout << "Warning: Adjusted number of elements and length. N0="
159                << n0 << ", l0=" << m_l0 << endl;
160        }
161        if ((n1+1)%m_NY > 0) {
162            double Dy=m_l1/n1;
163            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
164            m_l1=Dy*n1;
165            cout << "Warning: Adjusted number of elements and length. N1="
166                << n1 << ", l1=" << m_l1 << endl;
167        }
168        if ((n2+1)%m_NZ > 0) {
169            double Dz=m_l2/n2;
170            n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;
171            m_l2=Dz*n2;
172            cout << "Warning: Adjusted number of elements and length. N2="
173                << n2 << ", l2=" << m_l2 << endl;
174        }
175    
176        m_gNE0=n0;
177        m_gNE1=n1;
178        m_gNE2=n2;
179    
180      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2) || (m_NZ > 1 && (n2+1)/m_NZ<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2) || (m_NZ > 1 && (n2+1)/m_NZ<2))
181          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
# Line 121  bool Brick::operator==(const AbstractDom Line 245  bool Brick::operator==(const AbstractDom
245      return false;      return false;
246  }  }
247    
248    void Brick::readBinaryGrid(escript::Data& out, string filename,
249                               const vector<int>& first,
250                               const vector<int>& numValues) const
251    {
252        // check destination function space
253        int myN0, myN1, myN2;
254        if (out.getFunctionSpace().getTypeCode() == Nodes) {
255            myN0 = m_N0;
256            myN1 = m_N1;
257            myN2 = m_N2;
258        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
259                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
260            myN0 = m_NE0;
261            myN1 = m_NE1;
262            myN2 = m_NE2;
263        } else
264            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
265    
266        // check file existence and size
267        ifstream f(filename.c_str(), ifstream::binary);
268        if (f.fail()) {
269            throw RipleyException("readBinaryGrid(): cannot open file");
270        }
271        f.seekg(0, ios::end);
272        const int numComp = out.getDataPointSize();
273        const int filesize = f.tellg();
274        const int reqsize = numValues[0]*numValues[1]*numValues[2]*numComp*sizeof(float);
275        if (filesize < reqsize) {
276            f.close();
277            throw RipleyException("readBinaryGrid(): not enough data in file");
278        }
279    
280        // check if this rank contributes anything
281        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
282                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
283                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
284            f.close();
285            return;
286        }
287    
288        // now determine how much this rank has to write
289    
290        // first coordinates in data object to write to
291        const int first0 = max(0, first[0]-m_offset0);
292        const int first1 = max(0, first[1]-m_offset1);
293        const int first2 = max(0, first[2]-m_offset2);
294        // indices to first value in file
295        const int idx0 = max(0, m_offset0-first[0]);
296        const int idx1 = max(0, m_offset1-first[1]);
297        const int idx2 = max(0, m_offset2-first[2]);
298        // number of values to write
299        const int num0 = min(numValues[0]-idx0, myN0-first0);
300        const int num1 = min(numValues[1]-idx1, myN1-first1);
301        const int num2 = min(numValues[2]-idx2, myN2-first2);
302    
303        out.requireWrite();
304        vector<float> values(num0*numComp);
305        const int dpp = out.getNumDataPointsPerSample();
306    
307        for (index_t z=0; z<num2; z++) {
308            for (index_t y=0; y<num1; y++) {
309                const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
310                f.seekg(fileofs*sizeof(float));
311                f.read((char*)&values[0], num0*numComp*sizeof(float));
312                const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1;
313    
314                for (index_t x=0; x<num0; x++) {
315                    double* dest = out.getSampleDataRW(dataIndex+x);
316                    for (index_t c=0; c<numComp; c++) {
317                        if (!isnan(values[x*numComp+c])) {
318                            for (index_t q=0; q<dpp; q++) {
319                                *dest++ = static_cast<double>(values[x*numComp+c]);
320                            }
321                        }
322                    }
323                }
324            }
325        }
326    
327        f.close();
328    }
329    
330    void Brick::readNcGrid(escript::Data& out, string filename, string varname,
331                const vector<int>& first, const vector<int>& numValues) const
332    {
333    #ifdef USE_NETCDF
334        // check destination function space
335        int myN0, myN1, myN2;
336        if (out.getFunctionSpace().getTypeCode() == Nodes) {
337            myN0 = m_N0;
338            myN1 = m_N1;
339            myN2 = m_N2;
340        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
341                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
342            myN0 = m_NE0;
343            myN1 = m_NE1;
344            myN2 = m_NE2;
345        } else
346            throw RipleyException("readNcGrid(): invalid function space for output data object");
347    
348        if (first.size() != 3)
349            throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
350    
351        if (numValues.size() != 3)
352            throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
353    
354        // check file existence and size
355        NcFile f(filename.c_str(), NcFile::ReadOnly);
356        if (!f.is_valid())
357            throw RipleyException("readNcGrid(): cannot open file");
358    
359        NcVar* var = f.get_var(varname.c_str());
360        if (!var)
361            throw RipleyException("readNcGrid(): invalid variable name");
362    
363        // TODO: rank>0 data support
364        const int numComp = out.getDataPointSize();
365        if (numComp > 1)
366            throw RipleyException("readNcGrid(): only scalar data supported");
367    
368        const int dims = var->num_dims();
369        const long *edges = var->edges();
370    
371        // is this a slice of the data object (dims!=3)?
372        // note the expected ordering of edges (as in numpy: z,y,x)
373        if ( (dims==3 && (numValues[2] > edges[0] || numValues[1] > edges[1]
374                          || numValues[0] > edges[2]))
375                || (dims==2 && numValues[2]>1)
376                || (dims==1 && (numValues[2]>1 || numValues[1]>1)) ) {
377            throw RipleyException("readNcGrid(): not enough data in file");
378        }
379    
380        // check if this rank contributes anything
381        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
382                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
383                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
384            return;
385        }
386    
387        // now determine how much this rank has to write
388    
389        // first coordinates in data object to write to
390        const int first0 = max(0, first[0]-m_offset0);
391        const int first1 = max(0, first[1]-m_offset1);
392        const int first2 = max(0, first[2]-m_offset2);
393        // indices to first value in file
394        const int idx0 = max(0, m_offset0-first[0]);
395        const int idx1 = max(0, m_offset1-first[1]);
396        const int idx2 = max(0, m_offset2-first[2]);
397        // number of values to write
398        const int num0 = min(numValues[0]-idx0, myN0-first0);
399        const int num1 = min(numValues[1]-idx1, myN1-first1);
400        const int num2 = min(numValues[2]-idx2, myN2-first2);
401    
402        vector<double> values(num0*num1*num2);
403        if (dims==3) {
404            var->set_cur(idx2, idx1, idx0);
405            var->get(&values[0], num2, num1, num0);
406        } else if (dims==2) {
407            var->set_cur(idx1, idx0);
408            var->get(&values[0], num1, num0);
409        } else {
410            var->set_cur(idx0);
411            var->get(&values[0], num0);
412        }
413    
414        const int dpp = out.getNumDataPointsPerSample();
415        out.requireWrite();
416    
417        for (index_t z=0; z<num2; z++) {
418            for (index_t y=0; y<num1; y++) {
419    #pragma omp parallel for
420                for (index_t x=0; x<num0; x++) {
421                    const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1+x;
422                    const int srcIndex=z*num1*num0+y*num0+x;
423                    if (!isnan(values[srcIndex])) {
424                        double* dest = out.getSampleDataRW(dataIndex);
425                        for (index_t q=0; q<dpp; q++) {
426                            *dest++ = values[srcIndex];
427                        }
428                    }
429                }
430            }
431        }
432    #else
433        throw RipleyException("readNcGrid(): not compiled with netCDF support");
434    #endif
435    }
436    
437  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
438  {  {
439  #if USE_SILO  #if USE_SILO
# Line 528  void Brick::setToSize(escript::Data& out Line 841  void Brick::setToSize(escript::Data& out
841          const double xSize=getFirstCoordAndSpacing(0).second;          const double xSize=getFirstCoordAndSpacing(0).second;
842          const double ySize=getFirstCoordAndSpacing(1).second;          const double ySize=getFirstCoordAndSpacing(1).second;
843          const double zSize=getFirstCoordAndSpacing(2).second;          const double zSize=getFirstCoordAndSpacing(2).second;
844          const double size=min(min(xSize,ySize),zSize);          const double size=sqrt(xSize*xSize+ySize*ySize+zSize*zSize);
845  #pragma omp parallel for  #pragma omp parallel for
846          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
847              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 770  void Brick::assembleGradient(escript::Da Line 1083  void Brick::assembleGradient(escript::Da
1083    
1084      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
1085          out.requireWrite();          out.requireWrite();
1086  #pragma omp parallel for  #pragma omp parallel
1087          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
1088              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
1089                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
1090                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
1091                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
1092                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
1093                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
1094                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
1095                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_111(numComp);
1096                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
1097                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
1098                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1099                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1100                          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) / h0;                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1101                          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) / h0;                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1102                          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) / h0;                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1103                          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) / h0;                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1104                          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) / h1;                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1105                          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) / h1;                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1106                          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) / h1;                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1107                          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) / h1;                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1108                          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) / h2;                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1109                          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) / h2;                          for (index_t i=0; i < numComp; ++i) {
1110                          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) / h2;                              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) / h0;
1111                          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) / h2;                              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) / h0;
1112                          o[INDEX3(i,0,0,numComp,3)] = V0;                              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) / h0;
1113                          o[INDEX3(i,1,0,numComp,3)] = V4;                              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) / h0;
1114                          o[INDEX3(i,2,0,numComp,3)] = V8;                              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) / h1;
1115                          o[INDEX3(i,0,1,numComp,3)] = V0;                              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) / h1;
1116                          o[INDEX3(i,1,1,numComp,3)] = V5;                              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) / h1;
1117                          o[INDEX3(i,2,1,numComp,3)] = V9;                              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) / h1;
1118                          o[INDEX3(i,0,2,numComp,3)] = V1;                              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) / h2;
1119                          o[INDEX3(i,1,2,numComp,3)] = V4;                              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) / h2;
1120                          o[INDEX3(i,2,2,numComp,3)] = V10;                              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) / h2;
1121                          o[INDEX3(i,0,3,numComp,3)] = V1;                              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) / h2;
1122                          o[INDEX3(i,1,3,numComp,3)] = V5;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1123                          o[INDEX3(i,2,3,numComp,3)] = V11;                              o[INDEX3(i,1,0,numComp,3)] = V4;
1124                          o[INDEX3(i,0,4,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V8;
1125                          o[INDEX3(i,1,4,numComp,3)] = V6;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1126                          o[INDEX3(i,2,4,numComp,3)] = V8;                              o[INDEX3(i,1,1,numComp,3)] = V5;
1127                          o[INDEX3(i,0,5,numComp,3)] = V2;                              o[INDEX3(i,2,1,numComp,3)] = V9;
1128                          o[INDEX3(i,1,5,numComp,3)] = V7;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1129                          o[INDEX3(i,2,5,numComp,3)] = V9;                              o[INDEX3(i,1,2,numComp,3)] = V4;
1130                          o[INDEX3(i,0,6,numComp,3)] = V3;                              o[INDEX3(i,2,2,numComp,3)] = V10;
1131                          o[INDEX3(i,1,6,numComp,3)] = V6;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1132                          o[INDEX3(i,2,6,numComp,3)] = V10;                              o[INDEX3(i,1,3,numComp,3)] = V5;
1133                          o[INDEX3(i,0,7,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V11;
1134                          o[INDEX3(i,1,7,numComp,3)] = V7;                              o[INDEX3(i,0,4,numComp,3)] = V2;
1135                          o[INDEX3(i,2,7,numComp,3)] = V11;                              o[INDEX3(i,1,4,numComp,3)] = V6;
1136                      } // end of component loop i                              o[INDEX3(i,2,4,numComp,3)] = V8;
1137                  } // end of k0 loop                              o[INDEX3(i,0,5,numComp,3)] = V2;
1138              } // end of k1 loop                              o[INDEX3(i,1,5,numComp,3)] = V7;
1139          } // end of k2 loop                              o[INDEX3(i,2,5,numComp,3)] = V9;
1140                                o[INDEX3(i,0,6,numComp,3)] = V3;
1141                                o[INDEX3(i,1,6,numComp,3)] = V6;
1142                                o[INDEX3(i,2,6,numComp,3)] = V10;
1143                                o[INDEX3(i,0,7,numComp,3)] = V3;
1144                                o[INDEX3(i,1,7,numComp,3)] = V7;
1145                                o[INDEX3(i,2,7,numComp,3)] = V11;
1146                            } // end of component loop i
1147                        } // end of k0 loop
1148                    } // end of k1 loop
1149                } // end of k2 loop
1150            } // end of parallel section
1151      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
1152          out.requireWrite();          out.requireWrite();
1153  #pragma omp parallel for  #pragma omp parallel
1154          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
1155              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
1156                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
1157                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
1158                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
1159                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
1160                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
1161                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
1162                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
1163                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
1164                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
1165                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1166                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1167                          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 / h0;                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1168                          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 / h1;                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1169                          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 / h2;                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1170                      } // end of component loop i                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1171                  } // end of k0 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1172              } // end of k1 loop                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1173          } // end of k2 loop                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1174                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1175                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1176                            for (index_t i=0; i < numComp; ++i) {
1177                                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 / h0;
1178                                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 / h1;
1179                                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 / h2;
1180                            } // end of component loop i
1181                        } // end of k0 loop
1182                    } // end of k1 loop
1183                } // end of k2 loop
1184            } // end of parallel section
1185      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1186          out.requireWrite();          out.requireWrite();
1187  #pragma omp parallel  #pragma omp parallel
1188          {          {
1189                vector<double> f_000(numComp);
1190                vector<double> f_001(numComp);
1191                vector<double> f_010(numComp);
1192                vector<double> f_011(numComp);
1193                vector<double> f_100(numComp);
1194                vector<double> f_101(numComp);
1195                vector<double> f_110(numComp);
1196                vector<double> f_111(numComp);
1197              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1198  #pragma omp for nowait  #pragma omp for nowait
1199                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1200                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1201                          const double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1202                          const double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1203                          const double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1204                          const double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1205                          const double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1206                          const double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1207                          const double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1208                          const double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1209                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1210                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1211                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;
# Line 889  void Brick::assembleGradient(escript::Da Line 1232  void Brick::assembleGradient(escript::Da
1232  #pragma omp for nowait  #pragma omp for nowait
1233                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1234                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1235                          const double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1236                          const double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1237                          const double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1238                          const double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1239                          const double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1240                          const double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1241                          const double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1242                          const double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1243                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1244                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1245                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
# Line 923  void Brick::assembleGradient(escript::Da Line 1266  void Brick::assembleGradient(escript::Da
1266  #pragma omp for nowait  #pragma omp for nowait
1267                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1268                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1269                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1)), numComp*sizeof(double));
1270                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1271                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1)), numComp*sizeof(double));
1272                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1273                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1)), numComp*sizeof(double));
1274                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1275                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1)), numComp*sizeof(double));
1276                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1277                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1278                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1279                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;
# Line 956  void Brick::assembleGradient(escript::Da Line 1299  void Brick::assembleGradient(escript::Da
1299  #pragma omp for nowait  #pragma omp for nowait
1300                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1301                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1302                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1)), numComp*sizeof(double));
1303                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1304                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
1305                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1306                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1)), numComp*sizeof(double));
1307                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1308                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
1309                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1310                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1311                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1312                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
# Line 990  void Brick::assembleGradient(escript::Da Line 1333  void Brick::assembleGradient(escript::Da
1333  #pragma omp for nowait  #pragma omp for nowait
1334                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1335                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1336                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1)), numComp*sizeof(double));
1337                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1)), numComp*sizeof(double));
1338                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
1339                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1)), numComp*sizeof(double));
1340                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1)), numComp*sizeof(double));
1341                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1)), numComp*sizeof(double));
1342                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
1343                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1)), numComp*sizeof(double));
1344                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1345                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1346                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;
# Line 1024  void Brick::assembleGradient(escript::Da Line 1367  void Brick::assembleGradient(escript::Da
1367  #pragma omp for nowait  #pragma omp for nowait
1368                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1369                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1370                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1371                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1372                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1373                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1374                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1375                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1376                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1377                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1378                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1379                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1380                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
# Line 1059  void Brick::assembleGradient(escript::Da Line 1402  void Brick::assembleGradient(escript::Da
1402          out.requireWrite();          out.requireWrite();
1403  #pragma omp parallel  #pragma omp parallel
1404          {          {
1405                vector<double> f_000(numComp);
1406                vector<double> f_001(numComp);
1407                vector<double> f_010(numComp);
1408                vector<double> f_011(numComp);
1409                vector<double> f_100(numComp);
1410                vector<double> f_101(numComp);
1411                vector<double> f_110(numComp);
1412                vector<double> f_111(numComp);
1413              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1414  #pragma omp for nowait  #pragma omp for nowait
1415                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1416                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1417                          const double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1418                          const double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1419                          const double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1420                          const double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1421                          const double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1422                          const double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1423                          const double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1424                          const double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1425                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1426                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1427                              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 / h0;                              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 / h0;
# Line 1084  void Brick::assembleGradient(escript::Da Line 1435  void Brick::assembleGradient(escript::Da
1435  #pragma omp for nowait  #pragma omp for nowait
1436                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1437                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1438                          const double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1439                          const double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1440                          const double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1441                          const double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1442                          const double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1443                          const double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1444                          const double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1445                          const double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1446                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1447                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1448                              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 / h0;                              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 / h0;
# Line 1105  void Brick::assembleGradient(escript::Da Line 1456  void Brick::assembleGradient(escript::Da
1456  #pragma omp for nowait  #pragma omp for nowait
1457                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1458                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1459                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1)), numComp*sizeof(double));
1460                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1461                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1)), numComp*sizeof(double));
1462                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1463                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1)), numComp*sizeof(double));
1464                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1465                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1)), numComp*sizeof(double));
1466                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1467                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1468                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1469                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
# Line 1126  void Brick::assembleGradient(escript::Da Line 1477  void Brick::assembleGradient(escript::Da
1477  #pragma omp for nowait  #pragma omp for nowait
1478                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1479                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1480                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1)), numComp*sizeof(double));
1481                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1482                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
1483                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1484                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1)), numComp*sizeof(double));
1485                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1486                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
1487                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1488                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1489                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1490                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
# Line 1147  void Brick::assembleGradient(escript::Da Line 1498  void Brick::assembleGradient(escript::Da
1498  #pragma omp for nowait  #pragma omp for nowait
1499                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1500                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1501                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1)), numComp*sizeof(double));
1502                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1)), numComp*sizeof(double));
1503                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
1504                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1)), numComp*sizeof(double));
1505                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1)), numComp*sizeof(double));
1506                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1)), numComp*sizeof(double));
1507                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
1508                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1)), numComp*sizeof(double));
1509                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1510                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1511                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
# Line 1168  void Brick::assembleGradient(escript::Da Line 1519  void Brick::assembleGradient(escript::Da
1519  #pragma omp for nowait  #pragma omp for nowait
1520                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1521                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1522                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1523                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1524                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1525                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1526                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1527                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1528                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1)), numComp*sizeof(double));
1529                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
1530                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1531                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1532                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
# Line 1199  void Brick::assembleIntegrate(vector<dou Line 1550  void Brick::assembleIntegrate(vector<dou
1550      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1551      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1552      const index_t front = (m_offset2==0 ? 0 : 1);      const index_t front = (m_offset2==0 ? 0 : 1);
1553      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs = arg.getFunctionSpace().getTypeCode();
1554        if (fs == Elements && arg.actsExpanded()) {
1555          const double w_0 = h0*h1*h2/8.;          const double w_0 = h0*h1*h2/8.;
1556  #pragma omp parallel  #pragma omp parallel
1557          {          {
# Line 1228  void Brick::assembleIntegrate(vector<dou Line 1580  void Brick::assembleIntegrate(vector<dou
1580              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1581                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1582          } // end of parallel section          } // end of parallel section
1583      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1584        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1585          const double w_0 = h0*h1*h2;          const double w_0 = h0*h1*h2;
1586  #pragma omp parallel  #pragma omp parallel
1587          {          {
# Line 1249  void Brick::assembleIntegrate(vector<dou Line 1602  void Brick::assembleIntegrate(vector<dou
1602              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1603                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1604          } // end of parallel section          } // end of parallel section
1605      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1606        } else if (fs == FaceElements && arg.actsExpanded()) {
1607          const double w_0 = h1*h2/4.;          const double w_0 = h1*h2/4.;
1608          const double w_1 = h0*h2/4.;          const double w_1 = h0*h2/4.;
1609          const double w_2 = h0*h1/4.;          const double w_2 = h0*h1/4.;
# Line 1357  void Brick::assembleIntegrate(vector<dou Line 1711  void Brick::assembleIntegrate(vector<dou
1711                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1712          } // end of parallel section          } // end of parallel section
1713    
1714      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1715          const double w_0 = h1*h2;          const double w_0 = h1*h2;
1716          const double w_1 = h0*h2;          const double w_1 = h0*h2;
1717          const double w_2 = h0*h1;          const double w_2 = h0*h1;
# Line 1440  void Brick::assembleIntegrate(vector<dou Line 1794  void Brick::assembleIntegrate(vector<dou
1794              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1795                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1796          } // end of parallel section          } // end of parallel section
1797        } // function space selector
     }  
1798  }  }
1799    
1800  //protected  //protected
# Line 1521  void Brick::dofToNodes(escript::Data& ou Line 1874  void Brick::dofToNodes(escript::Data& ou
1874                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1875          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1876      }      }
1877        Paso_Coupler_free(coupler);
1878  }  }
1879    
1880  //private  //private
# Line 1763  void Brick::createPattern() Line 2117  void Brick::createPattern()
2117                          offsetInShared.push_back(offsetInShared.back()+nDOF0);                          offsetInShared.push_back(offsetInShared.back()+nDOF0);
2118                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
2119                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2120                          const int firstNode=(i1+1)/2*m_N0*(m_N1-1)                          const int firstNode=left+(i1+1)/2*m_N0*(m_N1-1)
2121                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2122                          for (dim_t i=0; i<nDOF0; i++, numShared++) {                          for (dim_t i=0; i<nDOF0; i++, numShared++) {
2123                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
# Line 1780  void Brick::createPattern() Line 2134  void Brick::createPattern()
2134                          offsetInShared.push_back(offsetInShared.back()+nDOF1);                          offsetInShared.push_back(offsetInShared.back()+nDOF1);
2135                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2136                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2137                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=bottom*m_N0
2138                                                +(i0+1)/2*(m_N0-1)
2139                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2140                          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          for (dim_t i=0; i<nDOF1; i++, numShared++) {
2141                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
# Line 1797  void Brick::createPattern() Line 2152  void Brick::createPattern()
2152                          offsetInShared.push_back(offsetInShared.back()+nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF2);
2153                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2154                                             +(i1+1)/2*nDOF0*(nDOF1-1);                                             +(i1+1)/2*nDOF0*(nDOF1-1);
2155                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=front*m_N0*m_N1
2156                                                +(i0+1)/2*(m_N0-1)
2157                                              +(i1+1)/2*m_N0*(m_N1-1);                                              +(i1+1)/2*m_N0*(m_N1-1);
2158                          for (dim_t i=0; i<nDOF2; i++, numShared++) {                          for (dim_t i=0; i<nDOF2; i++, numShared++) {
2159                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
# Line 1956  void Brick::interpolateNodesOnElements(e Line 2312  void Brick::interpolateNodesOnElements(e
2312      if (reduced) {      if (reduced) {
2313          out.requireWrite();          out.requireWrite();
2314          const double c0 = .125;          const double c0 = .125;
2315  #pragma omp parallel for  #pragma omp parallel
2316          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2317              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2318                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2319                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2320                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2321                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2322                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2323                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2324                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2325                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
2326                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
2327                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2328                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2329                          o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2330                      } // end of component loop i                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2331                  } // end of k0 loop                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2332              } // end of k1 loop                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2333          } // end of k2 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2334                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2335                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2336                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2337                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2338                            for (index_t i=0; i < numComp; ++i) {
2339                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);
2340                            } // end of component loop i
2341                        } // end of k0 loop
2342                    } // end of k1 loop
2343                } // end of k2 loop
2344            } // end of parallel section
2345      } else {      } else {
2346          out.requireWrite();          out.requireWrite();
2347          const double c0 = .0094373878376559314545;          const double c0 = .0094373878376559314545;
2348          const double c1 = .035220810900864519624;          const double c1 = .035220810900864519624;
2349          const double c2 = .13144585576580214704;          const double c2 = .13144585576580214704;
2350          const double c3 = .49056261216234406855;          const double c3 = .49056261216234406855;
2351  #pragma omp parallel for  #pragma omp parallel
2352          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2353              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2354                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2355                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2356                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2357                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2358                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2359                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2360                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2361                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));  #pragma omp for
2362                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
2363                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2364                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2365                          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]);                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2366                          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]);                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2367                          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]);                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2368                          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]);                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2369                          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]);                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2370                          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]);                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2371                          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]);                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2372                          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]);                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2373                      } // end of component loop i                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2374                  } // end of k0 loop                          for (index_t i=0; i < numComp; ++i) {
2375              } // end of k1 loop                              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]);
2376          } // end of k2 loop                              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]);
2377                                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]);
2378                                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]);
2379                                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]);
2380                                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]);
2381                                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]);
2382                                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]);
2383                            } // end of component loop i
2384                        } // end of k0 loop
2385                    } // end of k1 loop
2386                } // end of k2 loop
2387            } // end of parallel section
2388      }      }
2389  }  }
2390    
# Line 2020  void Brick::interpolateNodesOnFaces(escr Line 2398  void Brick::interpolateNodesOnFaces(escr
2398          const double c0 = .25;          const double c0 = .25;
2399  #pragma omp parallel  #pragma omp parallel
2400          {          {
2401                vector<double> f_000(numComp);
2402                vector<double> f_001(numComp);
2403                vector<double> f_010(numComp);
2404                vector<double> f_011(numComp);
2405                vector<double> f_100(numComp);
2406                vector<double> f_101(numComp);
2407                vector<double> f_110(numComp);
2408                vector<double> f_111(numComp);
2409              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2410  #pragma omp for nowait  #pragma omp for nowait
2411                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2412                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2413                          const double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2414                          const double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2415                          const double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2416                          const double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2417                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2418                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2419                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);
# Line 2039  void Brick::interpolateNodesOnFaces(escr Line 2425  void Brick::interpolateNodesOnFaces(escr
2425  #pragma omp for nowait  #pragma omp for nowait
2426                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2427                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2428                          const double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2429                          const double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2430                          const double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2431                          const double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2432                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2433                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2434                              o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
# Line 2054  void Brick::interpolateNodesOnFaces(escr Line 2440  void Brick::interpolateNodesOnFaces(escr
2440  #pragma omp for nowait  #pragma omp for nowait
2441                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2442                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2443                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1)), numComp*sizeof(double));
2444                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2445                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1)), numComp*sizeof(double));
2446                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2447                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2448                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2449                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
# Line 2069  void Brick::interpolateNodesOnFaces(escr Line 2455  void Brick::interpolateNodesOnFaces(escr
2455  #pragma omp for nowait  #pragma omp for nowait
2456                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2457                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2458                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
2459                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2460                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
2461                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2462                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2463                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2464                              o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
# Line 2084  void Brick::interpolateNodesOnFaces(escr Line 2470  void Brick::interpolateNodesOnFaces(escr
2470  #pragma omp for nowait  #pragma omp for nowait
2471                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2472                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2473                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1)), numComp*sizeof(double));
2474                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
2475                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1)), numComp*sizeof(double));
2476                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
2477                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2478                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2479                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
# Line 2099  void Brick::interpolateNodesOnFaces(escr Line 2485  void Brick::interpolateNodesOnFaces(escr
2485  #pragma omp for nowait  #pragma omp for nowait
2486                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2487                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2488                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2489                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2490                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2491                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2492                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2493                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2494                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);                              o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
# Line 2118  void Brick::interpolateNodesOnFaces(escr Line 2504  void Brick::interpolateNodesOnFaces(escr
2504          const double c2 = 0.62200846792814621559;          const double c2 = 0.62200846792814621559;
2505  #pragma omp parallel  #pragma omp parallel
2506          {          {
2507                vector<double> f_000(numComp);
2508                vector<double> f_001(numComp);
2509                vector<double> f_010(numComp);
2510                vector<double> f_011(numComp);
2511                vector<double> f_100(numComp);
2512                vector<double> f_101(numComp);
2513                vector<double> f_110(numComp);
2514                vector<double> f_111(numComp);
2515              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2516  #pragma omp for nowait  #pragma omp for nowait
2517                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2518                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2519                          const double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2520                          const double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2521                          const double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2522                          const double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2523                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2524                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2525                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);
# Line 2140  void Brick::interpolateNodesOnFaces(escr Line 2534  void Brick::interpolateNodesOnFaces(escr
2534  #pragma omp for nowait  #pragma omp for nowait
2535                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2536                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2537                          const double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2538                          const double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2539                          const double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2540                          const double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2541                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2542                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2543                              o[INDEX2(i,numComp,0)] = f_100[i]*c2 + f_111[i]*c0 + c1*(f_101[i] + f_110[i]);                              o[INDEX2(i,numComp,0)] = f_100[i]*c2 + f_111[i]*c0 + c1*(f_101[i] + f_110[i]);
# Line 2158  void Brick::interpolateNodesOnFaces(escr Line 2552  void Brick::interpolateNodesOnFaces(escr
2552  #pragma omp for nowait  #pragma omp for nowait
2553                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2554                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2555                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1)), numComp*sizeof(double));
2556                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2557                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1)), numComp*sizeof(double));
2558                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2559                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2560                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2561                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_100[i]);                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_100[i]);
# Line 2176  void Brick::interpolateNodesOnFaces(escr Line 2570  void Brick::interpolateNodesOnFaces(escr
2570  #pragma omp for nowait  #pragma omp for nowait
2571                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2572                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2573                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
2574                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2575                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1)), numComp*sizeof(double));
2576                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2577                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2578                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2579                              o[INDEX2(i,numComp,0)] = f_010[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_110[i]);                              o[INDEX2(i,numComp,0)] = f_010[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_110[i]);
# Line 2194  void Brick::interpolateNodesOnFaces(escr Line 2588  void Brick::interpolateNodesOnFaces(escr
2588  #pragma omp for nowait  #pragma omp for nowait
2589                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2590                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2591                          const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1)), numComp*sizeof(double));
2592                          const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
2593                          const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1)), numComp*sizeof(double));
2594                          const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1)), numComp*sizeof(double));
2595                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2596                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2597                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_100[i]);                              o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_100[i]);
# Line 2212  void Brick::interpolateNodesOnFaces(escr Line 2606  void Brick::interpolateNodesOnFaces(escr
2606  #pragma omp for nowait  #pragma omp for nowait
2607                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2608                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2609                          const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2610                          const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2611                          const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2612                          const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1)), numComp*sizeof(double));
2613                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2614                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2615                              o[INDEX2(i,numComp,0)] = f_001[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_101[i]);                              o[INDEX2(i,numComp,0)] = f_001[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_101[i]);

Legend:
Removed from v.3796  
changed lines
  Added in v.4174

  ViewVC Help
Powered by ViewVC 1.1.26