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

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

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

branches/ripleygmg_from_3668/ripley/src/Brick.cpp revision 3791 by caltinay, Wed Feb 1 05:10:22 2012 UTC trunk/ripley/src/Brick.cpp revision 3981 by jfenwick, Fri Sep 21 02:47:54 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            d1=(int)(d0*n1/(float)n0);
58            d2=m_mpiInfo->size/(d0*d1);
59            if (d0*d1*d2 != m_mpiInfo->size) {
60                // ratios not the same so leave "smallest" side undivided and try
61                // dividing 2 sides only
62                if (n0>=n1) {
63                    if (n1>=n2) {
64                        d0=d1=0;
65                        d2=1;
66                    } else {
67                        d0=d2=0;
68                        d1=1;
69                    }
70                } else {
71                    if (n0>=n2) {
72                        d0=d1=0;
73                        d2=1;
74                    } else {
75                        d0=1;
76                        d1=d2=0;
77                    }
78                }
79            }
80        }
81        if (d0<=0 && d1<=0) {
82            warn=true;
83            d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));
84            d1=m_mpiInfo->size/d0;
85            if (d0*d1*d2 != m_mpiInfo->size) {
86                // ratios not the same so subdivide side with more elements only
87                if (n0>n1) {
88                    d0=0;
89                    d1=1;
90                } else {
91                    d0=1;
92                    d1=0;
93                }
94            }
95        } else if (d0<=0 && d2<=0) {
96            warn=true;
97            d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1)));
98            d2=m_mpiInfo->size/d0;
99            if (d0*d1*d2 != m_mpiInfo->size) {
100                // ratios not the same so subdivide side with more elements only
101                if (n0>n2) {
102                    d0=0;
103                    d2=1;
104                } else {
105                    d0=1;
106                    d2=0;
107                }
108            }
109        } else if (d1<=0 && d2<=0) {
110            warn=true;
111            d1=int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1)));
112            d2=m_mpiInfo->size/d1;
113            if (d0*d1*d2 != m_mpiInfo->size) {
114                // ratios not the same so subdivide side with more elements only
115                if (n1>n2) {
116                    d1=0;
117                    d2=1;
118                } else {
119                    d1=1;
120                    d2=0;
121                }
122            }
123        }
124        if (d0<=0) {
125            // d1,d2 are preset, determine d0
126            d0=m_mpiInfo->size/(d1*d2);
127        } else if (d1<=0) {
128            // d0,d2 are preset, determine d1
129            d1=m_mpiInfo->size/(d0*d2);
130        } else if (d2<=0) {
131            // d0,d1 are preset, determine d2
132            d2=m_mpiInfo->size/(d0*d1);
133        }
134    
135        m_NX=d0;
136        m_NY=d1;
137        m_NZ=d2;
138    
139      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
140      // among number of ranks      // among number of ranks
141      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
142          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
143    
144      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0 || (n2+1)%m_NZ > 0)      if (warn) {
145          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
146                << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
147        }
148    
149        if ((n0+1)%m_NX > 0) {
150            double Dx=m_l0/n0;
151            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
152            m_l0=Dx*n0;
153            cout << "Warning: Adjusted number of elements and length. N0="
154                << n0 << ", l0=" << m_l0 << endl;
155        }
156        if ((n1+1)%m_NY > 0) {
157            double Dy=m_l1/n1;
158            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
159            m_l1=Dy*n1;
160            cout << "Warning: Adjusted number of elements and length. N1="
161                << n1 << ", l1=" << m_l1 << endl;
162        }
163        if ((n2+1)%m_NZ > 0) {
164            double Dz=m_l2/n2;
165            n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;
166            m_l2=Dz*n2;
167            cout << "Warning: Adjusted number of elements and length. N2="
168                << n2 << ", l2=" << m_l2 << endl;
169        }
170    
171        m_gNE0=n0;
172        m_gNE1=n1;
173        m_gNE2=n2;
174    
175      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))
176          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 240  bool Brick::operator==(const AbstractDom
240      return false;      return false;
241  }  }
242    
243    void Brick::readBinaryGrid(escript::Data& out, string filename,
244                               const vector<int>& first,
245                               const vector<int>& numValues) const
246    {
247        // check destination function space
248        int myN0, myN1, myN2;
249        if (out.getFunctionSpace().getTypeCode() == Nodes) {
250            myN0 = m_N0;
251            myN1 = m_N1;
252            myN2 = m_N2;
253        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
254                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
255            myN0 = m_NE0;
256            myN1 = m_NE1;
257            myN2 = m_NE2;
258        } else
259            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
260    
261        // check file existence and size
262        ifstream f(filename.c_str(), ifstream::binary);
263        if (f.fail()) {
264            throw RipleyException("readBinaryGrid(): cannot open file");
265        }
266        f.seekg(0, ios::end);
267        const int numComp = out.getDataPointSize();
268        const int filesize = f.tellg();
269        const int reqsize = numValues[0]*numValues[1]*numValues[2]*numComp*sizeof(float);
270        if (filesize < reqsize) {
271            f.close();
272            throw RipleyException("readBinaryGrid(): not enough data in file");
273        }
274    
275        // check if this rank contributes anything
276        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
277                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
278                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
279            f.close();
280            return;
281        }
282    
283        // now determine how much this rank has to write
284    
285        // first coordinates in data object to write to
286        const int first0 = max(0, first[0]-m_offset0);
287        const int first1 = max(0, first[1]-m_offset1);
288        const int first2 = max(0, first[2]-m_offset2);
289        // indices to first value in file
290        const int idx0 = max(0, m_offset0-first[0]);
291        const int idx1 = max(0, m_offset1-first[1]);
292        const int idx2 = max(0, m_offset2-first[2]);
293        // number of values to write
294        const int num0 = min(numValues[0]-idx0, myN0-first0);
295        const int num1 = min(numValues[1]-idx1, myN1-first1);
296        const int num2 = min(numValues[2]-idx2, myN2-first2);
297    
298        out.requireWrite();
299        vector<float> values(num0*numComp);
300        const int dpp = out.getNumDataPointsPerSample();
301    
302        for (index_t z=0; z<num2; z++) {
303            for (index_t y=0; y<num1; y++) {
304                const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
305                f.seekg(fileofs*sizeof(float));
306                f.read((char*)&values[0], num0*numComp*sizeof(float));
307                for (index_t x=0; x<num0; x++) {
308                    double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0+(first2+z)*myN0*myN1);
309                    for (index_t c=0; c<numComp; c++) {
310                        for (index_t q=0; q<dpp; q++) {
311                            *dest++ = static_cast<double>(values[x*numComp+c]);
312                        }
313                    }
314                }
315            }
316        }
317    
318        f.close();
319    }
320    
321  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
322  {  {
323  #if USE_SILO  #if USE_SILO
# Line 129  void Brick::dump(const string& fileName) Line 326  void Brick::dump(const string& fileName)
326          fn+=".silo";          fn+=".silo";
327      }      }
328    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
329      int driver=DB_HDF5;          int driver=DB_HDF5;    
330      string siloPath;      string siloPath;
331      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
332    
333  #ifdef ESYS_MPI  #ifdef ESYS_MPI
334      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
335        const int NUM_SILO_FILES = 1;
336        const char* blockDirFmt = "/block%04d";
337  #endif  #endif
338    
339      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 528  void Brick::setToSize(escript::Data& out Line 725  void Brick::setToSize(escript::Data& out
725          const double xSize=getFirstCoordAndSpacing(0).second;          const double xSize=getFirstCoordAndSpacing(0).second;
726          const double ySize=getFirstCoordAndSpacing(1).second;          const double ySize=getFirstCoordAndSpacing(1).second;
727          const double zSize=getFirstCoordAndSpacing(2).second;          const double zSize=getFirstCoordAndSpacing(2).second;
728          const double size=min(min(xSize,ySize),zSize);          const double size=sqrt(xSize*xSize+ySize*ySize+zSize*zSize);
729  #pragma omp parallel for  #pragma omp parallel for
730          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
731              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 759  void Brick::assembleGradient(escript::Da Line 956  void Brick::assembleGradient(escript::Da
956      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
957      const double h0 = m_l0/m_gNE0;      const double h0 = m_l0/m_gNE0;
958      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
959      const double h2 = m_l1/m_gNE2;      const double h2 = m_l2/m_gNE2;
960      const double C0 = .044658198738520451079;      const double C0 = .044658198738520451079;
961      const double C1 = .16666666666666666667;      const double C1 = .16666666666666666667;
962      const double C2 = .21132486540518711775;      const double C2 = .21132486540518711775;
# Line 770  void Brick::assembleGradient(escript::Da Line 967  void Brick::assembleGradient(escript::Da
967    
968      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
969          out.requireWrite();          out.requireWrite();
970  #pragma omp parallel for  #pragma omp parallel
971          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
972              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
973                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
974                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
975                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
976                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
977                      const double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
978                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
979                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_111(numComp);
980                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
981                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              for (index_t k2=0; k2 < m_NE2; ++k2) {
982                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
983                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
984                          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));
985                          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));
986                          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));
987                          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));
988                          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));
989                          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));
990                          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));
991                          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));
992                          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));
993                          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) {
994                          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;
995                          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;
996                          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;
997                          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;
998                          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;
999                          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;
1000                          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;
1001                          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;
1002                          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;
1003                          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;
1004                          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;
1005                          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;
1006                          o[INDEX3(i,1,3,numComp,3)] = V5;                              o[INDEX3(i,0,0,numComp,3)] = V0;
1007                          o[INDEX3(i,2,3,numComp,3)] = V11;                              o[INDEX3(i,1,0,numComp,3)] = V4;
1008                          o[INDEX3(i,0,4,numComp,3)] = V2;                              o[INDEX3(i,2,0,numComp,3)] = V8;
1009                          o[INDEX3(i,1,4,numComp,3)] = V6;                              o[INDEX3(i,0,1,numComp,3)] = V0;
1010                          o[INDEX3(i,2,4,numComp,3)] = V8;                              o[INDEX3(i,1,1,numComp,3)] = V5;
1011                          o[INDEX3(i,0,5,numComp,3)] = V2;                              o[INDEX3(i,2,1,numComp,3)] = V9;
1012                          o[INDEX3(i,1,5,numComp,3)] = V7;                              o[INDEX3(i,0,2,numComp,3)] = V1;
1013                          o[INDEX3(i,2,5,numComp,3)] = V9;                              o[INDEX3(i,1,2,numComp,3)] = V4;
1014                          o[INDEX3(i,0,6,numComp,3)] = V3;                              o[INDEX3(i,2,2,numComp,3)] = V10;
1015                          o[INDEX3(i,1,6,numComp,3)] = V6;                              o[INDEX3(i,0,3,numComp,3)] = V1;
1016                          o[INDEX3(i,2,6,numComp,3)] = V10;                              o[INDEX3(i,1,3,numComp,3)] = V5;
1017                          o[INDEX3(i,0,7,numComp,3)] = V3;                              o[INDEX3(i,2,3,numComp,3)] = V11;
1018                          o[INDEX3(i,1,7,numComp,3)] = V7;                              o[INDEX3(i,0,4,numComp,3)] = V2;
1019                          o[INDEX3(i,2,7,numComp,3)] = V11;                              o[INDEX3(i,1,4,numComp,3)] = V6;
1020                      } // end of component loop i                              o[INDEX3(i,2,4,numComp,3)] = V8;
1021                  } // end of k0 loop                              o[INDEX3(i,0,5,numComp,3)] = V2;
1022              } // end of k1 loop                              o[INDEX3(i,1,5,numComp,3)] = V7;
1023          } // end of k2 loop                              o[INDEX3(i,2,5,numComp,3)] = V9;
1024                                o[INDEX3(i,0,6,numComp,3)] = V3;
1025                                o[INDEX3(i,1,6,numComp,3)] = V6;
1026                                o[INDEX3(i,2,6,numComp,3)] = V10;
1027                                o[INDEX3(i,0,7,numComp,3)] = V3;
1028                                o[INDEX3(i,1,7,numComp,3)] = V7;
1029                                o[INDEX3(i,2,7,numComp,3)] = V11;
1030                            } // end of component loop i
1031                        } // end of k0 loop
1032                    } // end of k1 loop
1033                } // end of k2 loop
1034            } // end of parallel section
1035      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
1036          out.requireWrite();          out.requireWrite();
1037  #pragma omp parallel for  #pragma omp parallel
1038          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
1039              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
1040                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
1041                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
1042                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
1043                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
1044                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
1045                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
1046                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
1047                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
1048                      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) {
1049                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1050                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1051                          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));
1052                          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));
1053                          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));
1054                      } // end of component loop i                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1055                  } // end of k0 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
1056              } // end of k1 loop                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1057          } // end of k2 loop                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
1058                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
1059                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1060                            for (index_t i=0; i < numComp; ++i) {
1061                                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;
1062                                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;
1063                                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;
1064                            } // end of component loop i
1065                        } // end of k0 loop
1066                    } // end of k1 loop
1067                } // end of k2 loop
1068            } // end of parallel section
1069      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1070          out.requireWrite();          out.requireWrite();
1071  #pragma omp parallel  #pragma omp parallel
1072          {          {
1073                vector<double> f_000(numComp);
1074                vector<double> f_001(numComp);
1075                vector<double> f_010(numComp);
1076                vector<double> f_011(numComp);
1077                vector<double> f_100(numComp);
1078                vector<double> f_101(numComp);
1079                vector<double> f_110(numComp);
1080                vector<double> f_111(numComp);
1081              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1082  #pragma omp for nowait  #pragma omp for nowait
1083                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1084                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1085                          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));
1086                          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));
1087                          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));
1088                          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));
1089                          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));
1090                          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));
1091                          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));
1092                          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));
1093                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1094                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1095                              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 1116  void Brick::assembleGradient(escript::Da
1116  #pragma omp for nowait  #pragma omp for nowait
1117                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1118                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1119                          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));
1120                          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));
1121                          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));
1122                          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));
1123                          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));
1124                          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));
1125                          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));
1126                          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));
1127                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1128                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1129                              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 1150  void Brick::assembleGradient(escript::Da
1150  #pragma omp for nowait  #pragma omp for nowait
1151                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1152                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1153                          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));
1154                          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));
1155                          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));
1156                          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));
1157                          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));
1158                          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));
1159                          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));
1160                          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));
1161                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1162                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1163                              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 1183  void Brick::assembleGradient(escript::Da
1183  #pragma omp for nowait  #pragma omp for nowait
1184                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1185                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1186                          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));
1187                          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));
1188                          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));
1189                          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));
1190                          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));
1191                          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));
1192                          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));
1193                          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));
1194                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1195                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1196                              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 1217  void Brick::assembleGradient(escript::Da
1217  #pragma omp for nowait  #pragma omp for nowait
1218                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1219                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1220                          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));
1221                          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));
1222                          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));
1223                          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));
1224                          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));
1225                          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));
1226                          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));
1227                          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));
1228                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1229                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1230                              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 1251  void Brick::assembleGradient(escript::Da
1251  #pragma omp for nowait  #pragma omp for nowait
1252                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1253                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1254                          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));
1255                          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));
1256                          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));
1257                          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));
1258                          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));
1259                          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));
1260                          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));
1261                          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));
1262                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1263                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1264                              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 1286  void Brick::assembleGradient(escript::Da
1286          out.requireWrite();          out.requireWrite();
1287  #pragma omp parallel  #pragma omp parallel
1288          {          {
1289                vector<double> f_000(numComp);
1290                vector<double> f_001(numComp);
1291                vector<double> f_010(numComp);
1292                vector<double> f_011(numComp);
1293                vector<double> f_100(numComp);
1294                vector<double> f_101(numComp);
1295                vector<double> f_110(numComp);
1296                vector<double> f_111(numComp);
1297              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1298  #pragma omp for nowait  #pragma omp for nowait
1299                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1300                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1301                          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));
1302                          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));
1303                          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));
1304                          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));
1305                          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));
1306                          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));
1307                          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));
1308                          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));
1309                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1310                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1311                              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 1319  void Brick::assembleGradient(escript::Da
1319  #pragma omp for nowait  #pragma omp for nowait
1320                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1321                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
1322                          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));
1323                          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));
1324                          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));
1325                          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));
1326                          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));
1327                          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));
1328                          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));
1329                          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));
1330                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1331                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1332                              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 1340  void Brick::assembleGradient(escript::Da
1340  #pragma omp for nowait  #pragma omp for nowait
1341                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1342                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1343                          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));
1344                          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));
1345                          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));
1346                          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));
1347                          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));
1348                          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));
1349                          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));
1350                          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));
1351                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1352                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1353                              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 1361  void Brick::assembleGradient(escript::Da
1361  #pragma omp for nowait  #pragma omp for nowait
1362                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
1363                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1364                          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));
1365                          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));
1366                          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));
1367                          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));
1368                          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));
1369                          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));
1370                          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));
1371                          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));
1372                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1373                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1374                              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 1382  void Brick::assembleGradient(escript::Da
1382  #pragma omp for nowait  #pragma omp for nowait
1383                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1384                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1385                          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));
1386                          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));
1387                          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));
1388                          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));
1389                          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));
1390                          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));
1391                          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));
1392                          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));
1393                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1394                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1395                              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 1403  void Brick::assembleGradient(escript::Da
1403  #pragma omp for nowait  #pragma omp for nowait
1404                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1405                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
1406                          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));
1407                          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));
1408                          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));
1409                          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));
1410                          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));
1411                          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));
1412                          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));
1413                          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));
1414                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1415                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1416                              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 1434  void Brick::assembleIntegrate(vector<dou
1434      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
1435      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
1436      const index_t front = (m_offset2==0 ? 0 : 1);      const index_t front = (m_offset2==0 ? 0 : 1);
1437      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs = arg.getFunctionSpace().getTypeCode();
1438        if (fs == Elements && arg.actsExpanded()) {
1439          const double w_0 = h0*h1*h2/8.;          const double w_0 = h0*h1*h2/8.;
1440  #pragma omp parallel  #pragma omp parallel
1441          {          {
# Line 1228  void Brick::assembleIntegrate(vector<dou Line 1464  void Brick::assembleIntegrate(vector<dou
1464              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1465                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1466          } // end of parallel section          } // end of parallel section
1467      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
1468        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
1469          const double w_0 = h0*h1*h2;          const double w_0 = h0*h1*h2;
1470  #pragma omp parallel  #pragma omp parallel
1471          {          {
# Line 1249  void Brick::assembleIntegrate(vector<dou Line 1486  void Brick::assembleIntegrate(vector<dou
1486              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1487                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1488          } // end of parallel section          } // end of parallel section
1489      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1490        } else if (fs == FaceElements && arg.actsExpanded()) {
1491          const double w_0 = h1*h2/4.;          const double w_0 = h1*h2/4.;
1492          const double w_1 = h0*h2/4.;          const double w_1 = h0*h2/4.;
1493          const double w_2 = h0*h1/4.;          const double w_2 = h0*h1/4.;
# Line 1357  void Brick::assembleIntegrate(vector<dou Line 1595  void Brick::assembleIntegrate(vector<dou
1595                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1596          } // end of parallel section          } // end of parallel section
1597    
1598      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1599          const double w_0 = h1*h2;          const double w_0 = h1*h2;
1600          const double w_1 = h0*h2;          const double w_1 = h0*h2;
1601          const double w_2 = h0*h1;          const double w_2 = h0*h1;
# Line 1440  void Brick::assembleIntegrate(vector<dou Line 1678  void Brick::assembleIntegrate(vector<dou
1678              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1679                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1680          } // end of parallel section          } // end of parallel section
1681        } // function space selector
     }  
1682  }  }
1683    
1684  //protected  //protected
# Line 1763  void Brick::createPattern() Line 2000  void Brick::createPattern()
2000                          offsetInShared.push_back(offsetInShared.back()+nDOF0);                          offsetInShared.push_back(offsetInShared.back()+nDOF0);
2001                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)                          const int firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
2002                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2003                          const int firstNode=(i1+1)/2*m_N0*(m_N1-1)                          const int firstNode=left+(i1+1)/2*m_N0*(m_N1-1)
2004                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2005                          for (dim_t i=0; i<nDOF0; i++, numShared++) {                          for (dim_t i=0; i<nDOF0; i++, numShared++) {
2006                              sendShared.push_back(firstDOF+i);                              sendShared.push_back(firstDOF+i);
# Line 1780  void Brick::createPattern() Line 2017  void Brick::createPattern()
2017                          offsetInShared.push_back(offsetInShared.back()+nDOF1);                          offsetInShared.push_back(offsetInShared.back()+nDOF1);
2018                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2019                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);                                             +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
2020                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=bottom*m_N0
2021                                                +(i0+1)/2*(m_N0-1)
2022                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);                                              +(i2+1)/2*m_N0*m_N1*(m_N2-1);
2023                          for (dim_t i=0; i<nDOF1; i++, numShared++) {                          for (dim_t i=0; i<nDOF1; i++, numShared++) {
2024                              sendShared.push_back(firstDOF+i*nDOF0);                              sendShared.push_back(firstDOF+i*nDOF0);
# Line 1797  void Brick::createPattern() Line 2035  void Brick::createPattern()
2035                          offsetInShared.push_back(offsetInShared.back()+nDOF2);                          offsetInShared.push_back(offsetInShared.back()+nDOF2);
2036                          const int firstDOF=(i0+1)/2*(nDOF0-1)                          const int firstDOF=(i0+1)/2*(nDOF0-1)
2037                                             +(i1+1)/2*nDOF0*(nDOF1-1);                                             +(i1+1)/2*nDOF0*(nDOF1-1);
2038                          const int firstNode=(i0+1)/2*(m_N0-1)                          const int firstNode=front*m_N0*m_N1
2039                                                +(i0+1)/2*(m_N0-1)
2040                                              +(i1+1)/2*m_N0*(m_N1-1);                                              +(i1+1)/2*m_N0*(m_N1-1);
2041                          for (dim_t i=0; i<nDOF2; i++, numShared++) {                          for (dim_t i=0; i<nDOF2; i++, numShared++) {
2042                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);                              sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
# Line 1956  void Brick::interpolateNodesOnElements(e Line 2195  void Brick::interpolateNodesOnElements(e
2195      if (reduced) {      if (reduced) {
2196          out.requireWrite();          out.requireWrite();
2197          const double c0 = .125;          const double c0 = .125;
2198  #pragma omp parallel for  #pragma omp parallel
2199          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2200              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2201                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2202                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2203                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2204                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2205                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2206                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2207                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2208                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));  #pragma omp for
2209                      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) {
2210                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2211                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2212                          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));
2213                      } // end of component loop i                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2214                  } // end of k0 loop                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2215              } // end of k1 loop                          memcpy(&f_011[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2216          } // end of k2 loop                          memcpy(&f_100[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1)), numComp*sizeof(double));
2217                            memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2218                            memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1)), numComp*sizeof(double));
2219                            memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1)), numComp*sizeof(double));
2220                            double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2221                            for (index_t i=0; i < numComp; ++i) {
2222                                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]);
2223                            } // end of component loop i
2224                        } // end of k0 loop
2225                    } // end of k1 loop
2226                } // end of k2 loop
2227            } // end of parallel section
2228      } else {      } else {
2229          out.requireWrite();          out.requireWrite();
2230          const double c0 = .0094373878376559314545;          const double c0 = .0094373878376559314545;
2231          const double c1 = .035220810900864519624;          const double c1 = .035220810900864519624;
2232          const double c2 = .13144585576580214704;          const double c2 = .13144585576580214704;
2233          const double c3 = .49056261216234406855;          const double c3 = .49056261216234406855;
2234  #pragma omp parallel for  #pragma omp parallel
2235          for (index_t k2=0; k2 < m_NE2; ++k2) {          {
2236              for (index_t k1=0; k1 < m_NE1; ++k1) {              vector<double> f_000(numComp);
2237                  for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_001(numComp);
2238                      const double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));              vector<double> f_010(numComp);
2239                      const double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));              vector<double> f_011(numComp);
2240                      const double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));              vector<double> f_100(numComp);
2241                      const double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));              vector<double> f_101(numComp);
2242                      const double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));              vector<double> f_110(numComp);
2243                      const double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));              vector<double> f_111(numComp);
2244                      const double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));  #pragma omp for
2245                      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) {
2246                      double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2247                      for (index_t i=0; i < numComp; ++i) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2248                          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));
2249                          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));
2250                          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));
2251                          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));
2252                          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));
2253                          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));
2254                          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));
2255                          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));
2256                      } // end of component loop i                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
2257                  } // end of k0 loop                          for (index_t i=0; i < numComp; ++i) {
2258              } // 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]);
2259          } // 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]);
2260                                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]);
2261                                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]);
2262                                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]);
2263                                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]);
2264                                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]);
2265                                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]);
2266                            } // end of component loop i
2267                        } // end of k0 loop
2268                    } // end of k1 loop
2269                } // end of k2 loop
2270            } // end of parallel section
2271      }      }
2272  }  }
2273    
# Line 2020  void Brick::interpolateNodesOnFaces(escr Line 2281  void Brick::interpolateNodesOnFaces(escr
2281          const double c0 = .25;          const double c0 = .25;
2282  #pragma omp parallel  #pragma omp parallel
2283          {          {
2284                vector<double> f_000(numComp);
2285                vector<double> f_001(numComp);
2286                vector<double> f_010(numComp);
2287                vector<double> f_011(numComp);
2288                vector<double> f_100(numComp);
2289                vector<double> f_101(numComp);
2290                vector<double> f_110(numComp);
2291                vector<double> f_111(numComp);
2292              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2293  #pragma omp for nowait  #pragma omp for nowait
2294                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2295                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2296                          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));
2297                          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));
2298                          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));
2299                          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));
2300                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2301                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2302                              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 2308  void Brick::interpolateNodesOnFaces(escr
2308  #pragma omp for nowait  #pragma omp for nowait
2309                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2310                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2311                          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));
2312                          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));
2313                          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));
2314                          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));
2315                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2316                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2317                              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 2323  void Brick::interpolateNodesOnFaces(escr
2323  #pragma omp for nowait  #pragma omp for nowait
2324                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2325                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2326                          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));
2327                          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));
2328                          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));
2329                          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));
2330                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2331                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2332                              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 2338  void Brick::interpolateNodesOnFaces(escr
2338  #pragma omp for nowait  #pragma omp for nowait
2339                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2340                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2341                          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));
2342                          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));
2343                          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));
2344                          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));
2345                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2346                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2347                              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 2353  void Brick::interpolateNodesOnFaces(escr
2353  #pragma omp for nowait  #pragma omp for nowait
2354                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2355                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2356                          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));
2357                          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));
2358                          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));
2359                          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));
2360                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2361                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2362                              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 2368  void Brick::interpolateNodesOnFaces(escr
2368  #pragma omp for nowait  #pragma omp for nowait
2369                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2370                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2371                          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));
2372                          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));
2373                          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));
2374                          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));
2375                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2376                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2377                              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 2387  void Brick::interpolateNodesOnFaces(escr
2387          const double c2 = 0.62200846792814621559;          const double c2 = 0.62200846792814621559;
2388  #pragma omp parallel  #pragma omp parallel
2389          {          {
2390                vector<double> f_000(numComp);
2391                vector<double> f_001(numComp);
2392                vector<double> f_010(numComp);
2393                vector<double> f_011(numComp);
2394                vector<double> f_100(numComp);
2395                vector<double> f_101(numComp);
2396                vector<double> f_110(numComp);
2397                vector<double> f_111(numComp);
2398              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
2399  #pragma omp for nowait  #pragma omp for nowait
2400                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2401                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2402                          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));
2403                          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));
2404                          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));
2405                          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));
2406                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
2407                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2408                              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 2417  void Brick::interpolateNodesOnFaces(escr
2417  #pragma omp for nowait  #pragma omp for nowait
2418                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2419                      for (index_t k1=0; k1 < m_NE1; ++k1) {                      for (index_t k1=0; k1 < m_NE1; ++k1) {
2420                          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));
2421                          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));
2422                          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));
2423                          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));
2424                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
2425                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2426                              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 2435  void Brick::interpolateNodesOnFaces(escr
2435  #pragma omp for nowait  #pragma omp for nowait
2436                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2437                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2438                          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));
2439                          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));
2440                          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));
2441                          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));
2442                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
2443                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2444                              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 2453  void Brick::interpolateNodesOnFaces(escr
2453  #pragma omp for nowait  #pragma omp for nowait
2454                  for (index_t k2=0; k2 < m_NE2; ++k2) {                  for (index_t k2=0; k2 < m_NE2; ++k2) {
2455                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2456                          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));
2457                          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));
2458                          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));
2459                          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));
2460                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
2461                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2462                              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 2471  void Brick::interpolateNodesOnFaces(escr
2471  #pragma omp for nowait  #pragma omp for nowait
2472                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2473                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2474                          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));
2475                          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));
2476                          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));
2477                          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));
2478                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
2479                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2480                              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 2489  void Brick::interpolateNodesOnFaces(escr
2489  #pragma omp for nowait  #pragma omp for nowait
2490                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
2491                      for (index_t k0=0; k0 < m_NE0; ++k0) {                      for (index_t k0=0; k0 < m_NE0; ++k0) {
2492                          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));
2493                          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));
2494                          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));
2495                          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));
2496                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
2497                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
2498                              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]);
# Line 9438  void Brick::assemblePDEBoundarySingle(Pa Line 9715  void Brick::assemblePDEBoundarySingle(Pa
9715      {      {
9716          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
9717              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9718  #pragma omp for nowait  #pragma omp for
9719                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9720                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
9721                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 9562  void Brick::assemblePDEBoundarySingle(Pa Line 9839  void Brick::assemblePDEBoundarySingle(Pa
9839    
9840          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
9841              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9842  #pragma omp for nowait  #pragma omp for
9843                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9844                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
9845                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 9686  void Brick::assemblePDEBoundarySingle(Pa Line 9963  void Brick::assemblePDEBoundarySingle(Pa
9963    
9964          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
9965              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
9966  #pragma omp for nowait  #pragma omp for
9967                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
9968                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
9969                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 9810  void Brick::assemblePDEBoundarySingle(Pa Line 10087  void Brick::assemblePDEBoundarySingle(Pa
10087    
10088          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
10089              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10090  #pragma omp for nowait  #pragma omp for
10091                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10092                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10093                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 9934  void Brick::assemblePDEBoundarySingle(Pa Line 10211  void Brick::assemblePDEBoundarySingle(Pa
10211    
10212          if (m_faceOffset[4] > -1) {          if (m_faceOffset[4] > -1) {
10213              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10214  #pragma omp for nowait  #pragma omp for
10215                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10216                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10217                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10058  void Brick::assemblePDEBoundarySingle(Pa Line 10335  void Brick::assemblePDEBoundarySingle(Pa
10335    
10336          if (m_faceOffset[5] > -1) {          if (m_faceOffset[5] > -1) {
10337              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10338  #pragma omp for nowait  #pragma omp for
10339                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10340                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10341                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10202  void Brick::assemblePDEBoundarySingleRed Line 10479  void Brick::assemblePDEBoundarySingleRed
10479      {      {
10480          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
10481              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10482  #pragma omp for nowait  #pragma omp for
10483                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10484                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
10485                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10254  void Brick::assemblePDEBoundarySingleRed Line 10531  void Brick::assemblePDEBoundarySingleRed
10531    
10532          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
10533              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10534  #pragma omp for nowait  #pragma omp for
10535                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10536                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
10537                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10306  void Brick::assemblePDEBoundarySingleRed Line 10583  void Brick::assemblePDEBoundarySingleRed
10583    
10584          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
10585              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10586  #pragma omp for nowait  #pragma omp for
10587                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10588                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10589                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10357  void Brick::assemblePDEBoundarySingleRed Line 10634  void Brick::assemblePDEBoundarySingleRed
10634    
10635          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
10636              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10637  #pragma omp for nowait  #pragma omp for
10638                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10639                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10640                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10409  void Brick::assemblePDEBoundarySingleRed Line 10686  void Brick::assemblePDEBoundarySingleRed
10686    
10687          if (m_faceOffset[4] > -1) {          if (m_faceOffset[4] > -1) {
10688              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10689  #pragma omp for nowait  #pragma omp for
10690                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10691                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10692                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10461  void Brick::assemblePDEBoundarySingleRed Line 10738  void Brick::assemblePDEBoundarySingleRed
10738    
10739          if (m_faceOffset[5] > -1) {          if (m_faceOffset[5] > -1) {
10740              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
10741  #pragma omp for nowait  #pragma omp for
10742                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
10743                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
10744                          vector<double> EM_S(8*8, 0);                          vector<double> EM_S(8*8, 0);
# Line 10570  void Brick::assemblePDEBoundarySystem(Pa Line 10847  void Brick::assemblePDEBoundarySystem(Pa
10847      {      {
10848          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
10849              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10850  #pragma omp for nowait  #pragma omp for
10851                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10852                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
10853                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 10708  void Brick::assemblePDEBoundarySystem(Pa Line 10985  void Brick::assemblePDEBoundarySystem(Pa
10985    
10986          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
10987              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
10988  #pragma omp for nowait  #pragma omp for
10989                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
10990                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
10991                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 10846  void Brick::assemblePDEBoundarySystem(Pa Line 11123  void Brick::assemblePDEBoundarySystem(Pa
11123    
11124          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
11125              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11126  #pragma omp for nowait  #pragma omp for
11127                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11128                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11129                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 10984  void Brick::assemblePDEBoundarySystem(Pa Line 11261  void Brick::assemblePDEBoundarySystem(Pa
11261    
11262          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
11263              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11264  #pragma omp for nowait  #pragma omp for
11265                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11266                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11267                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11122  void Brick::assemblePDEBoundarySystem(Pa Line 11399  void Brick::assemblePDEBoundarySystem(Pa
11399    
11400          if (m_faceOffset[4] > -1) {          if (m_faceOffset[4] > -1) {
11401              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11402  #pragma omp for nowait  #pragma omp for
11403                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11404                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11405                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11260  void Brick::assemblePDEBoundarySystem(Pa Line 11537  void Brick::assemblePDEBoundarySystem(Pa
11537    
11538          if (m_faceOffset[5] > -1) {          if (m_faceOffset[5] > -1) {
11539              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11540  #pragma omp for nowait  #pragma omp for
11541                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11542                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11543                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11425  void Brick::assemblePDEBoundarySystemRed Line 11702  void Brick::assemblePDEBoundarySystemRed
11702      {      {
11703          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
11704              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11705  #pragma omp for nowait  #pragma omp for
11706                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11707                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
11708                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11483  void Brick::assemblePDEBoundarySystemRed Line 11760  void Brick::assemblePDEBoundarySystemRed
11760    
11761          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
11762              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11763  #pragma omp for nowait  #pragma omp for
11764                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11765                      for (index_t k1=0; k1<m_NE1; ++k1) {                      for (index_t k1=0; k1<m_NE1; ++k1) {
11766                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11541  void Brick::assemblePDEBoundarySystemRed Line 11818  void Brick::assemblePDEBoundarySystemRed
11818    
11819          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
11820              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11821  #pragma omp for nowait  #pragma omp for
11822                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11823                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11824                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11599  void Brick::assemblePDEBoundarySystemRed Line 11876  void Brick::assemblePDEBoundarySystemRed
11876    
11877          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
11878              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring              for (index_t k2_0=0; k2_0<2; k2_0++) { // colouring
11879  #pragma omp for nowait  #pragma omp for
11880                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {                  for (index_t k2=k2_0; k2<m_NE2; k2+=2) {
11881                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11882                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11657  void Brick::assemblePDEBoundarySystemRed Line 11934  void Brick::assemblePDEBoundarySystemRed
11934    
11935          if (m_faceOffset[4] > -1) {          if (m_faceOffset[4] > -1) {
11936              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11937  #pragma omp for nowait  #pragma omp for
11938                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11939                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11940                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);
# Line 11715  void Brick::assemblePDEBoundarySystemRed Line 11992  void Brick::assemblePDEBoundarySystemRed
11992    
11993          if (m_faceOffset[5] > -1) {          if (m_faceOffset[5] > -1) {
11994              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
11995  #pragma omp for nowait  #pragma omp for
11996                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
11997                      for (index_t k0=0; k0<m_NE0; ++k0) {                      for (index_t k0=0; k0<m_NE0; ++k0) {
11998                          vector<double> EM_S(8*8*numEq*numComp, 0);                          vector<double> EM_S(8*8*numEq*numComp, 0);

Legend:
Removed from v.3791  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26