/[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 4010 by caltinay, Tue Oct 2 06:57:11 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2012 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" {
# Line 32  namespace ripley { Line 34  namespace ripley {
34  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,
35               double x1, double y1, double z1, int d0, int d1, int d2) :               double x1, double y1, double z1, int d0, int d1, int d2) :
36      RipleyDomain(3),      RipleyDomain(3),
     m_gNE0(n0),  
     m_gNE1(n1),  
     m_gNE2(n2),  
37      m_x0(x0),      m_x0(x0),
38      m_y0(y0),      m_y0(y0),
39      m_z0(z0),      m_z0(z0),
40      m_l0(x1-x0),      m_l0(x1-x0),
41      m_l1(y1-y0),      m_l1(y1-y0),
42      m_l2(z1-z0),      m_l2(z1-z0)
     m_NX(d0),  
     m_NY(d1),  
     m_NZ(d2)  
43  {  {
44        // ignore subdivision parameters for serial run
45        if (m_mpiInfo->size == 1) {
46            d0=1;
47            d1=1;
48            d2=1;
49        }
50    
51        bool warn=false;
52        // if number of subdivisions is non-positive, try to subdivide by the same
53        // ratio as the number of elements
54        if (d0<=0 && d1<=0 && d2<=0) {
55            warn=true;
56            d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));
57            d0=max(1, d0);
58            d1=max(1, (int)(d0*n1/(float)n0));
59            d2=m_mpiInfo->size/(d0*d1);
60            if (d0*d1*d2 != m_mpiInfo->size) {
61                // ratios not the same so leave "smallest" side undivided and try
62                // dividing 2 sides only
63                if (n0>=n1) {
64                    if (n1>=n2) {
65                        d0=d1=0;
66                        d2=1;
67                    } else {
68                        d0=d2=0;
69                        d1=1;
70                    }
71                } else {
72                    if (n0>=n2) {
73                        d0=d1=0;
74                        d2=1;
75                    } else {
76                        d0=1;
77                        d1=d2=0;
78                    }
79                }
80            }
81        }
82        if (d0<=0 && d1<=0) {
83            warn=true;
84            d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
85            d1=m_mpiInfo->size/d0;
86            if (d0*d1*d2 != m_mpiInfo->size) {
87                // ratios not the same so subdivide side with more elements only
88                if (n0>n1) {
89                    d0=0;
90                    d1=1;
91                } else {
92                    d0=1;
93                    d1=0;
94                }
95            }
96        } else if (d0<=0 && d2<=0) {
97            warn=true;
98            d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));
99            d2=m_mpiInfo->size/d0;
100            if (d0*d1*d2 != m_mpiInfo->size) {
101                // ratios not the same so subdivide side with more elements only
102                if (n0>n2) {
103                    d0=0;
104                    d2=1;
105                } else {
106                    d0=1;
107                    d2=0;
108                }
109            }
110        } else if (d1<=0 && d2<=0) {
111            warn=true;
112            d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));
113            d2=m_mpiInfo->size/d1;
114            if (d0*d1*d2 != m_mpiInfo->size) {
115                // ratios not the same so subdivide side with more elements only
116                if (n1>n2) {
117                    d1=0;
118                    d2=1;
119                } else {
120                    d1=1;
121                    d2=0;
122                }
123            }
124        }
125        if (d0<=0) {
126            // d1,d2 are preset, determine d0
127            d0=m_mpiInfo->size/(d1*d2);
128        } else if (d1<=0) {
129            // d0,d2 are preset, determine d1
130            d1=m_mpiInfo->size/(d0*d2);
131        } else if (d2<=0) {
132            // d0,d1 are preset, determine d2
133            d2=m_mpiInfo->size/(d0*d1);
134        }
135    
136        m_NX=d0;
137        m_NY=d1;
138        m_NZ=d2;
139    
140      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
141      // among number of ranks      // among number of ranks
142      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
143          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
144    
145      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0 || (n2+1)%m_NZ > 0)      if (warn) {
146          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
147                << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
148        }
149    
150        if ((n0+1)%m_NX > 0) {
151            double Dx=m_l0/n0;
152            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
153            m_l0=Dx*n0;
154            cout << "Warning: Adjusted number of elements and length. N0="
155                << n0 << ", l0=" << m_l0 << endl;
156        }
157        if ((n1+1)%m_NY > 0) {
158            double Dy=m_l1/n1;
159            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
160            m_l1=Dy*n1;
161            cout << "Warning: Adjusted number of elements and length. N1="
162                << n1 << ", l1=" << m_l1 << endl;
163        }
164        if ((n2+1)%m_NZ > 0) {
165            double Dz=m_l2/n2;
166            n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;
167            m_l2=Dz*n2;
168            cout << "Warning: Adjusted number of elements and length. N2="
169                << n2 << ", l2=" << m_l2 << endl;
170        }
171    
172        m_gNE0=n0;
173        m_gNE1=n1;
174        m_gNE2=n2;
175    
176      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))
177          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 241  bool Brick::operator==(const AbstractDom
241      return false;      return false;
242  }  }
243    
244    void Brick::readBinaryGrid(escript::Data& out, string filename,
245                               const vector<int>& first,
246                               const vector<int>& numValues) const
247    {
248        // check destination function space
249        int myN0, myN1, myN2;
250        if (out.getFunctionSpace().getTypeCode() == Nodes) {
251            myN0 = m_N0;
252            myN1 = m_N1;
253            myN2 = m_N2;
254        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
255                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
256            myN0 = m_NE0;
257            myN1 = m_NE1;
258            myN2 = m_NE2;
259        } else
260            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
261    
262        // check file existence and size
263        ifstream f(filename.c_str(), ifstream::binary);
264        if (f.fail()) {
265            throw RipleyException("readBinaryGrid(): cannot open file");
266        }
267        f.seekg(0, ios::end);
268        const int numComp = out.getDataPointSize();
269        const int filesize = f.tellg();
270        const int reqsize = numValues[0]*numValues[1]*numValues[2]*numComp*sizeof(float);
271        if (filesize < reqsize) {
272            f.close();
273            throw RipleyException("readBinaryGrid(): not enough data in file");
274        }
275    
276        // check if this rank contributes anything
277        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
278                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
279                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
280            f.close();
281            return;
282        }
283    
284        // now determine how much this rank has to write
285    
286        // first coordinates in data object to write to
287        const int first0 = max(0, first[0]-m_offset0);
288        const int first1 = max(0, first[1]-m_offset1);
289        const int first2 = max(0, first[2]-m_offset2);
290        // indices to first value in file
291        const int idx0 = max(0, m_offset0-first[0]);
292        const int idx1 = max(0, m_offset1-first[1]);
293        const int idx2 = max(0, m_offset2-first[2]);
294        // number of values to write
295        const int num0 = min(numValues[0]-idx0, myN0-first0);
296        const int num1 = min(numValues[1]-idx1, myN1-first1);
297        const int num2 = min(numValues[2]-idx2, myN2-first2);
298    
299        out.requireWrite();
300        vector<float> values(num0*numComp);
301        const int dpp = out.getNumDataPointsPerSample();
302    
303        for (index_t z=0; z<num2; z++) {
304            for (index_t y=0; y<num1; y++) {
305                const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
306                f.seekg(fileofs*sizeof(float));
307                f.read((char*)&values[0], num0*numComp*sizeof(float));
308                const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1;
309    
310                for (index_t x=0; x<num0; x++) {
311                    double* dest = out.getSampleDataRW(dataIndex+x);
312                    for (index_t c=0; c<numComp; c++) {
313                        for (index_t q=0; q<dpp; q++) {
314                            *dest++ = static_cast<double>(values[x*numComp+c]);
315                        }
316                    }
317                }
318            }
319        }
320    
321        f.close();
322    }
323    
324  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
325  {  {
326  #if USE_SILO  #if USE_SILO
# Line 528  void Brick::setToSize(escript::Data& out Line 728  void Brick::setToSize(escript::Data& out
728          const double xSize=getFirstCoordAndSpacing(0).second;          const double xSize=getFirstCoordAndSpacing(0).second;
729          const double ySize=getFirstCoordAndSpacing(1).second;          const double ySize=getFirstCoordAndSpacing(1).second;
730          const double zSize=getFirstCoordAndSpacing(2).second;          const double zSize=getFirstCoordAndSpacing(2).second;
731          const double size=min(min(xSize,ySize),zSize);          const double size=sqrt(xSize*xSize+ySize*ySize+zSize*zSize);
732  #pragma omp parallel for  #pragma omp parallel for
733          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
734              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 770  void Brick::assembleGradient(escript::Da Line 970  void Brick::assembleGradient(escript::Da
970    
971      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
972          out.requireWrite();          out.requireWrite();
973  #pragma omp parallel for  #pragma omp parallel
974          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
975              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
976                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
977                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
978                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
979                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
980                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
981                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
982                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_111(numComp);
983                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
984                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
985                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
986                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
987                          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));
988                          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));
989                          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));
990                          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));
991                          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));
992                          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));
993                          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));
994                          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));
995                          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));
996                          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) {
997                          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;
998                          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;
999                          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;
1000                          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;
1001                          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;
1002                          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;
1003                          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;
1004                          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;
1005                          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;
1006                          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;
1007                          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;
1008                          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;
1009                          o[INDEX3(i,1,3,numComp,3)] = V5;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1010                          o[INDEX3(i,2,3,numComp,3)] = V11;                              o[INDEX3(i,1,0,numComp,3)] = V4;
1011                          o[INDEX3(i,0,4,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V8;
1012                          o[INDEX3(i,1,4,numComp,3)] = V6;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1013                          o[INDEX3(i,2,4,numComp,3)] = V8;                              o[INDEX3(i,1,1,numComp,3)] = V5;
1014                          o[INDEX3(i,0,5,numComp,3)] = V2;                              o[INDEX3(i,2,1,numComp,3)] = V9;
1015                          o[INDEX3(i,1,5,numComp,3)] = V7;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1016                          o[INDEX3(i,2,5,numComp,3)] = V9;                              o[INDEX3(i,1,2,numComp,3)] = V4;
1017                          o[INDEX3(i,0,6,numComp,3)] = V3;                              o[INDEX3(i,2,2,numComp,3)] = V10;
1018                          o[INDEX3(i,1,6,numComp,3)] = V6;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1019                          o[INDEX3(i,2,6,numComp,3)] = V10;                              o[INDEX3(i,1,3,numComp,3)] = V5;
1020                          o[INDEX3(i,0,7,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V11;
1021                          o[INDEX3(i,1,7,numComp,3)] = V7;                              o[INDEX3(i,0,4,numComp,3)] = V2;
1022                          o[INDEX3(i,2,7,numComp,3)] = V11;                              o[INDEX3(i,1,4,numComp,3)] = V6;
1023                      } // end of component loop i                              o[INDEX3(i,2,4,numComp,3)] = V8;
1024                  } // end of k0 loop                              o[INDEX3(i,0,5,numComp,3)] = V2;
1025              } // end of k1 loop                              o[INDEX3(i,1,5,numComp,3)] = V7;
1026          } // end of k2 loop                              o[INDEX3(i,2,5,numComp,3)] = V9;
1027                                o[INDEX3(i,0,6,numComp,3)] = V3;
1028                                o[INDEX3(i,1,6,numComp,3)] = V6;
1029                                o[INDEX3(i,2,6,numComp,3)] = V10;
1030                                o[INDEX3(i,0,7,numComp,3)] = V3;
1031                                o[INDEX3(i,1,7,numComp,3)] = V7;
1032                                o[INDEX3(i,2,7,numComp,3)] = V11;
1033                            } // end of component loop i
1034                        } // end of k0 loop
1035                    } // end of k1 loop
1036                } // end of k2 loop
1037            } // end of parallel section
1038      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
1039          out.requireWrite();          out.requireWrite();
1040  #pragma omp parallel for  #pragma omp parallel
1041          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
1042              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
1043                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
1044                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
1045                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
1046                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
1047                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
1048                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
1049                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
1050                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
1051                      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) {
1052                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1053                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1054                          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));
1055                          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));
1056                          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));
1057                      } // end of component loop i                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1058                  } // end of k0 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1059              } // end of k1 loop                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1060          } // end of k2 loop                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1061                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1062                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1063                            for (index_t i=0; i < numComp; ++i) {
1064                                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;
1065                                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;
1066                                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;
1067                            } // end of component loop i
1068                        } // end of k0 loop
1069                    } // end of k1 loop
1070                } // end of k2 loop
1071            } // end of parallel section
1072      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1073          out.requireWrite();          out.requireWrite();
1074  #pragma omp parallel  #pragma omp parallel
1075          {          {
1076                vector<double> f_000(numComp);
1077                vector<double> f_001(numComp);
1078                vector<double> f_010(numComp);
1079                vector<double> f_011(numComp);
1080                vector<double> f_100(numComp);
1081                vector<double> f_101(numComp);
1082                vector<double> f_110(numComp);
1083                vector<double> f_111(numComp);
1084              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1085  #pragma omp for nowait  #pragma omp for nowait
1086                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1087                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1088                          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));
1089                          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));
1090                          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));
1091                          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));
1092                          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));
1093                          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));
1094                          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));
1095                          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));
1096                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1097                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1098                              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 1119  void Brick::assembleGradient(escript::Da
1119  #pragma omp for nowait  #pragma omp for nowait
1120                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1121                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1122                          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));
1123                          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));
1124                          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));
1125                          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));
1126                          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));
1127                          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));
1128                          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));
1129                          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));
1130                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1131                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1132                              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 1153  void Brick::assembleGradient(escript::Da
1153  #pragma omp for nowait  #pragma omp for nowait
1154                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1155                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1156                          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));
1157                          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));
1158                          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));
1159                          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));
1160                          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));
1161                          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));
1162                          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));
1163                          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));
1164                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1165                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1166                              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 1186  void Brick::assembleGradient(escript::Da
1186  #pragma omp for nowait  #pragma omp for nowait
1187                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1188                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1189                          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));
1190                          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));
1191                          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));
1192                          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));
1193                          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));
1194                          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));
1195                          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));
1196                          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));
1197                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1198                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1199                              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 1220  void Brick::assembleGradient(escript::Da
1220  #pragma omp for nowait  #pragma omp for nowait
1221                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1222                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1223                          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));
1224                          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));
1225                          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));
1226                          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));
1227                          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));
1228                          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));
1229                          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));
1230                          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));
1231                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1232                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1233                              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 1254  void Brick::assembleGradient(escript::Da
1254  #pragma omp for nowait  #pragma omp for nowait
1255                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1256                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1257                          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));
1258                          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));
1259                          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));
1260                          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));
1261                          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));
1262                          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));
1263                          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));
1264                          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));
1265                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1266                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1267                              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 1289  void Brick::assembleGradient(escript::Da
1289          out.requireWrite();          out.requireWrite();
1290  #pragma omp parallel  #pragma omp parallel
1291          {          {
1292                vector<double> f_000(numComp);
1293                vector<double> f_001(numComp);
1294                vector<double> f_010(numComp);
1295                vector<double> f_011(numComp);
1296                vector<double> f_100(numComp);
1297                vector<double> f_101(numComp);
1298                vector<double> f_110(numComp);
1299                vector<double> f_111(numComp);
1300              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1301  #pragma omp for nowait  #pragma omp for nowait
1302                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1303                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1304                          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));
1305                          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));
1306                          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));
1307                          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));
1308                          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));
1309                          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));
1310                          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));
1311                          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));
1312                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1313                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1314                              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 1322  void Brick::assembleGradient(escript::Da
1322  #pragma omp for nowait  #pragma omp for nowait
1323                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1324                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1325                          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));
1326                          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));
1327                          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));
1328                          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));
1329                          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));
1330                          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));
1331                          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));
1332                          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));
1333                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1334                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1335                              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 1343  void Brick::assembleGradient(escript::Da
1343  #pragma omp for nowait  #pragma omp for nowait
1344                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1345                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1346                          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));
1347                          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));
1348                          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));
1349                          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));
1350                          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));
1351                          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));
1352                          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));
1353                          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));
1354                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1355                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1356                              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 1364  void Brick::assembleGradient(escript::Da
1364  #pragma omp for nowait  #pragma omp for nowait
1365                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1366                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1367                          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));
1368                          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));
1369                          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));
1370                          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));
1371                          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));
1372                          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));
1373                          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));
1374                          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));
1375                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1376                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1377                              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 1385  void Brick::assembleGradient(escript::Da
1385  #pragma omp for nowait  #pragma omp for nowait
1386                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1387                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1388                          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));
1389                          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));
1390                          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));
1391                          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));
1392                          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));
1393                          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));
1394                          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));
1395                          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));
1396                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1397                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1398                              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 1406  void Brick::assembleGradient(escript::Da
1406  #pragma omp for nowait  #pragma omp for nowait
1407                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1408                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1409                          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));
1410                          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));
1411                          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));
1412                          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));
1413                          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));
1414                          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));
1415                          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));
1416                          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));
1417                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1418                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1419                              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 1437  void Brick::assembleIntegrate(vector<dou
1437      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1438      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1439      const index_t front = (m_offset2==0 ? 0 : 1);      const index_t front = (m_offset2==0 ? 0 : 1);
1440      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs = arg.getFunctionSpace().getTypeCode();
1441        if (fs == Elements && arg.actsExpanded()) {
1442          const double w_0 = h0*h1*h2/8.;          const double w_0 = h0*h1*h2/8.;
1443  #pragma omp parallel  #pragma omp parallel
1444          {          {
# Line 1228  void Brick::assembleIntegrate(vector<dou Line 1467  void Brick::assembleIntegrate(vector<dou
1467              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1468                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1469          } // end of parallel section          } // end of parallel section
1470      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1471        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1472          const double w_0 = h0*h1*h2;          const double w_0 = h0*h1*h2;
1473  #pragma omp parallel  #pragma omp parallel
1474          {          {
# Line 1249  void Brick::assembleIntegrate(vector<dou Line 1489  void Brick::assembleIntegrate(vector<dou
1489              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1490                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1491          } // end of parallel section          } // end of parallel section
1492      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1493        } else if (fs == FaceElements && arg.actsExpanded()) {
1494          const double w_0 = h1*h2/4.;          const double w_0 = h1*h2/4.;
1495          const double w_1 = h0*h2/4.;          const double w_1 = h0*h2/4.;
1496          const double w_2 = h0*h1/4.;          const double w_2 = h0*h1/4.;
# Line 1357  void Brick::assembleIntegrate(vector<dou Line 1598  void Brick::assembleIntegrate(vector<dou
1598                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1599          } // end of parallel section          } // end of parallel section
1600    
1601      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1602          const double w_0 = h1*h2;          const double w_0 = h1*h2;
1603          const double w_1 = h0*h2;          const double w_1 = h0*h2;
1604          const double w_2 = h0*h1;          const double w_2 = h0*h1;
# Line 1440  void Brick::assembleIntegrate(vector<dou Line 1681  void Brick::assembleIntegrate(vector<dou
1681              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1682                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1683          } // end of parallel section          } // end of parallel section
1684        } // function space selector
     }  
1685  }  }
1686    
1687  //protected  //protected
# Line 1521  void Brick::dofToNodes(escript::Data& ou Line 1761  void Brick::dofToNodes(escript::Data& ou
1761                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1762          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1763      }      }
1764        Paso_Coupler_free(coupler);
1765  }  }
1766    
1767  //private  //private
# Line 1763  void Brick::createPattern() Line 2004  void Brick::createPattern()
2004                          offsetInShared.push_back(offsetInShared.back()+nDOF0);                          offsetInShared.push_back(offsetInShared.back()+nDOF0);
2005                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
2006                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2007                          const int firstNode=(i1+1)/2*m_N0*(m_N1-1)                          const int firstNode=left+(i1+1)/2*m_N0*(m_N1-1)
2008                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2009                          for (dim_t i=0; i<nDOF0; i++, numShared++) {                          for (dim_t i=0; i<nDOF0; i++, numShared++) {
2010                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
# Line 1780  void Brick::createPattern() Line 2021  void Brick::createPattern()
2021                          offsetInShared.push_back(offsetInShared.back()+nDOF1);                          offsetInShared.push_back(offsetInShared.back()+nDOF1);
2022                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2023                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2024                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=bottom*m_N0
2025                                                +(i0+1)/2*(m_N0-1)
2026                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2027                          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          for (dim_t i=0; i<nDOF1; i++, numShared++) {
2028                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
# Line 1797  void Brick::createPattern() Line 2039  void Brick::createPattern()
2039                          offsetInShared.push_back(offsetInShared.back()+nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF2);
2040                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2041                                             +(i1+1)/2*nDOF0*(nDOF1-1);                                             +(i1+1)/2*nDOF0*(nDOF1-1);
2042                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=front*m_N0*m_N1
2043                                                +(i0+1)/2*(m_N0-1)
2044                                              +(i1+1)/2*m_N0*(m_N1-1);                                              +(i1+1)/2*m_N0*(m_N1-1);
2045                          for (dim_t i=0; i<nDOF2; i++, numShared++) {                          for (dim_t i=0; i<nDOF2; i++, numShared++) {
2046                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
# Line 1956  void Brick::interpolateNodesOnElements(e Line 2199  void Brick::interpolateNodesOnElements(e
2199      if (reduced) {      if (reduced) {
2200          out.requireWrite();          out.requireWrite();
2201          const double c0 = .125;          const double c0 = .125;
2202  #pragma omp parallel for  #pragma omp parallel
2203          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2204              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2205                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2206                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2207                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2208                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2209                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2210                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2211                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2212                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
2213                      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) {
2214                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2215                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2216                          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));
2217                      } // end of component loop i                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2218                  } // end of k0 loop                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2219              } // end of k1 loop                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2220          } // end of k2 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2221                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2222                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2223                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2224                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2225                            for (index_t i=0; i < numComp; ++i) {
2226                                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]);
2227                            } // end of component loop i
2228                        } // end of k0 loop
2229                    } // end of k1 loop
2230                } // end of k2 loop
2231            } // end of parallel section
2232      } else {      } else {
2233          out.requireWrite();          out.requireWrite();
2234          const double c0 = .0094373878376559314545;          const double c0 = .0094373878376559314545;
2235          const double c1 = .035220810900864519624;          const double c1 = .035220810900864519624;
2236          const double c2 = .13144585576580214704;          const double c2 = .13144585576580214704;
2237          const double c3 = .49056261216234406855;          const double c3 = .49056261216234406855;
2238  #pragma omp parallel for  #pragma omp parallel
2239          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2240              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2241                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2242                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2243                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2244                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2245                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2246                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2247                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2248                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));  #pragma omp for
2249                      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) {
2250                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2251                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2252                          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));
2253                          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));
2254                          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));
2255                          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));
2256                          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));
2257                          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));
2258                          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));
2259                          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));
2260                      } // end of component loop i                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2261                  } // end of k0 loop                          for (index_t i=0; i < numComp; ++i) {
2262              } // 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]);
2263          } // 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]);
2264                                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]);
2265                                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]);
2266                                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]);
2267                                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]);
2268                                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]);
2269                                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]);
2270                            } // end of component loop i
2271                        } // end of k0 loop
2272                    } // end of k1 loop
2273                } // end of k2 loop
2274            } // end of parallel section
2275      }      }
2276  }  }
2277    
# Line 2020  void Brick::interpolateNodesOnFaces(escr Line 2285  void Brick::interpolateNodesOnFaces(escr
2285          const double c0 = .25;          const double c0 = .25;
2286  #pragma omp parallel  #pragma omp parallel
2287          {          {
2288                vector<double> f_000(numComp);
2289                vector<double> f_001(numComp);
2290                vector<double> f_010(numComp);
2291                vector<double> f_011(numComp);
2292                vector<double> f_100(numComp);
2293                vector<double> f_101(numComp);
2294                vector<double> f_110(numComp);
2295                vector<double> f_111(numComp);
2296              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2297  #pragma omp for nowait  #pragma omp for nowait
2298                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2299                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2300                          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));
2301                          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));
2302                          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));
2303                          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));
2304                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2305                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2306                              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 2312  void Brick::interpolateNodesOnFaces(escr
2312  #pragma omp for nowait  #pragma omp for nowait
2313                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2314                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2315                          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));
2316                          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));
2317                          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));
2318                          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));
2319                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2320                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2321                              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 2327  void Brick::interpolateNodesOnFaces(escr
2327  #pragma omp for nowait  #pragma omp for nowait
2328                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2329                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2330                          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));
2331                          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));
2332                          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));
2333                          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));
2334                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2335                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2336                              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 2342  void Brick::interpolateNodesOnFaces(escr
2342  #pragma omp for nowait  #pragma omp for nowait
2343                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2344                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2345                          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));
2346                          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));
2347                          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));
2348                          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));
2349                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2350                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2351                              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 2357  void Brick::interpolateNodesOnFaces(escr
2357  #pragma omp for nowait  #pragma omp for nowait
2358                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2359                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2360                          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));
2361                          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));
2362                          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));
2363                          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));
2364                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2365                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2366                              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 2372  void Brick::interpolateNodesOnFaces(escr
2372  #pragma omp for nowait  #pragma omp for nowait
2373                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2374                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2375                          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));
2376                          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));
2377                          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));
2378                          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));
2379                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2380                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2381                              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 2391  void Brick::interpolateNodesOnFaces(escr
2391          const double c2 = 0.62200846792814621559;          const double c2 = 0.62200846792814621559;
2392  #pragma omp parallel  #pragma omp parallel
2393          {          {
2394                vector<double> f_000(numComp);
2395                vector<double> f_001(numComp);
2396                vector<double> f_010(numComp);
2397                vector<double> f_011(numComp);
2398                vector<double> f_100(numComp);
2399                vector<double> f_101(numComp);
2400                vector<double> f_110(numComp);
2401                vector<double> f_111(numComp);
2402              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2403  #pragma omp for nowait  #pragma omp for nowait
2404                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2405                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2406                          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));
2407                          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));
2408                          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));
2409                          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));
2410                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2411                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2412                              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 2421  void Brick::interpolateNodesOnFaces(escr
2421  #pragma omp for nowait  #pragma omp for nowait
2422                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2423                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2424                          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));
2425                          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));
2426                          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));
2427                          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));
2428                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2429                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2430                              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 2439  void Brick::interpolateNodesOnFaces(escr
2439  #pragma omp for nowait  #pragma omp for nowait
2440                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2441                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2442                          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));
2443                          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));
2444                          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));
2445                          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));
2446                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2447                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2448                              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 2457  void Brick::interpolateNodesOnFaces(escr
2457  #pragma omp for nowait  #pragma omp for nowait
2458                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2459                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2460                          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));
2461                          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));
2462                          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));
2463                          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));
2464                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2465                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2466                              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 2475  void Brick::interpolateNodesOnFaces(escr
2475  #pragma omp for nowait  #pragma omp for nowait
2476                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2477                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2478                          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));
2479                          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));
2480                          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));
2481                          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));
2482                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2483                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2484                              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 2493  void Brick::interpolateNodesOnFaces(escr
2493  #pragma omp for nowait  #pragma omp for nowait
2494                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2495                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2496                          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));
2497                          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));
2498                          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));
2499                          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));
2500                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2501                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2502                              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.4010

  ViewVC Help
Powered by ViewVC 1.1.26