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

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

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

revision 4988 by caltinay, Wed Jun 4 01:15:10 2014 UTC revision 5084 by caltinay, Sun Jun 29 23:29:51 2014 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 #include <limits>  
   
17  #include <ripley/Brick.h>  #include <ripley/Brick.h>
 #include <esysUtils/esysFileWriter.h>  
18  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
19  #include <ripley/WaveAssembler3D.h>  #include <ripley/WaveAssembler3D.h>
20  #include <ripley/LameAssembler3D.h>  #include <ripley/LameAssembler3D.h>
21    #include <ripley/blocktools.h>
22  #include <ripley/domainhelpers.h>  #include <ripley/domainhelpers.h>
23    #include <esysUtils/esysFileWriter.h>
24    #include <esysUtils/EsysRandom.h>
25    #include <paso/SystemMatrix.h>
26    
27  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
28    
29  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 36  Line 38 
38  #endif  #endif
39    
40  #include <iomanip>  #include <iomanip>
41    #include <limits>
 #include "esysUtils/EsysRandom.h"  
 #include "blocktools.h"  
42    
43    
44  using namespace std;  using namespace std;
# Line 531  void Brick::readBinaryGridImpl(escript:: Line 531  void Brick::readBinaryGridImpl(escript::
531                                  if (params.byteOrder != BYTEORDER_NATIVE) {                                  if (params.byteOrder != BYTEORDER_NATIVE) {
532                                      char* cval = reinterpret_cast<char*>(&val);                                      char* cval = reinterpret_cast<char*>(&val);
533                                      // this will alter val!!                                      // this will alter val!!
534                                      byte_swap32(cval);                                      if (sizeof(ValueType)>4) {
535                                            byte_swap64(cval);
536                                        } else {
537                                            byte_swap32(cval);
538                                        }
539                                  }                                  }
540                                  if (!std::isnan(val)) {                                  if (!isnan(val)) {
541                                      for (int q=0; q<dpp; q++) {                                      for (int q=0; q<dpp; q++) {
542                                          *dest++ = static_cast<double>(val);                                          *dest++ = static_cast<double>(val);
543                                      }                                      }
# Line 653  void Brick::readBinaryGridZippedImpl(esc Line 657  void Brick::readBinaryGridZippedImpl(esc
657                                      // this will alter val!!                                      // this will alter val!!
658                                      byte_swap32(cval);                                      byte_swap32(cval);
659                                  }                                  }
660                                  if (!std::isnan(val)) {                                  if (!isnan(val)) {
661                                      for (int q=0; q<dpp; q++) {                                      for (int q=0; q<dpp; q++) {
662                                          *dest++ = static_cast<double>(val);                                          *dest++ = static_cast<double>(val);
663                                      }                                      }
# Line 738  void Brick::writeBinaryGridImpl(const es Line 742  void Brick::writeBinaryGridImpl(const es
742                      oss.write((char*)&fvalue, sizeof(fvalue));                      oss.write((char*)&fvalue, sizeof(fvalue));
743                  } else {                  } else {
744                      char* value = reinterpret_cast<char*>(&fvalue);                      char* value = reinterpret_cast<char*>(&fvalue);
745                      oss.write(byte_swap32(value), sizeof(fvalue));                      if (sizeof(fvalue)>4) {
746                            byte_swap64(value);
747                        } else {
748                            byte_swap32(value);
749                        }
750                        oss.write(value, sizeof(fvalue));
751                  }                  }
752              }              }
753              fw.writeAt(oss, fileofs);              fw.writeAt(oss, fileofs);
# Line 811  void Brick::dump(const string& fileName) Line 820  void Brick::dump(const string& fileName)
820      boost::scoped_ptr<double> y(new double[m_NN[1]]);      boost::scoped_ptr<double> y(new double[m_NN[1]]);
821      boost::scoped_ptr<double> z(new double[m_NN[2]]);      boost::scoped_ptr<double> z(new double[m_NN[2]]);
822      double* coords[3] = { x.get(), y.get(), z.get() };      double* coords[3] = { x.get(), y.get(), z.get() };
823        const int NN0 = m_NN[0];
824        const int NN1 = m_NN[1];
825        const int NN2 = m_NN[2];
826    
827  #pragma omp parallel  #pragma omp parallel
828      {      {
829  #pragma omp for  #pragma omp for
830          for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {          for (dim_t i0 = 0; i0 < NN0; i0++) {
831              coords[0][i0]=getLocalCoordinate(i0, 0);              coords[0][i0]=getLocalCoordinate(i0, 0);
832          }          }
833  #pragma omp for  #pragma omp for
834          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {          for (dim_t i1 = 0; i1 < NN1; i1++) {
835              coords[1][i1]=getLocalCoordinate(i1, 1);              coords[1][i1]=getLocalCoordinate(i1, 1);
836          }          }
837  #pragma omp for  #pragma omp for
838          for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {          for (dim_t i2 = 0; i2 < NN2; i2++) {
839              coords[2][i2]=getLocalCoordinate(i2, 2);              coords[2][i2]=getLocalCoordinate(i2, 2);
840          }          }
841      }      }
# Line 973  bool Brick::ownSample(int fsType, index_ Line 986  bool Brick::ownSample(int fsType, index_
986    
987  void Brick::setToNormal(escript::Data& out) const  void Brick::setToNormal(escript::Data& out) const
988  {  {
989        const int NE0 = m_NE[0];
990        const int NE1 = m_NE[1];
991        const int NE2 = m_NE[2];
992    
993      if (out.getFunctionSpace().getTypeCode() == FaceElements) {      if (out.getFunctionSpace().getTypeCode() == FaceElements) {
994          out.requireWrite();          out.requireWrite();
995  #pragma omp parallel  #pragma omp parallel
996          {          {
997              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
998  #pragma omp for nowait  #pragma omp for nowait
999                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1000                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1001                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1002                          // set vector at four quadrature points                          // set vector at four quadrature points
1003                          *o++ = -1.; *o++ = 0.; *o++ = 0.;                          *o++ = -1.; *o++ = 0.; *o++ = 0.;
# Line 993  void Brick::setToNormal(escript::Data& o Line 1010  void Brick::setToNormal(escript::Data& o
1010    
1011              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1012  #pragma omp for nowait  #pragma omp for nowait
1013                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1014                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1015                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1016                          // set vector at four quadrature points                          // set vector at four quadrature points
1017                          *o++ = 1.; *o++ = 0.; *o++ = 0.;                          *o++ = 1.; *o++ = 0.; *o++ = 0.;
# Line 1007  void Brick::setToNormal(escript::Data& o Line 1024  void Brick::setToNormal(escript::Data& o
1024    
1025              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1026  #pragma omp for nowait  #pragma omp for nowait
1027                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1028                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1029                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1030                          // set vector at four quadrature points                          // set vector at four quadrature points
1031                          *o++ = 0.; *o++ = -1.; *o++ = 0.;                          *o++ = 0.; *o++ = -1.; *o++ = 0.;
# Line 1021  void Brick::setToNormal(escript::Data& o Line 1038  void Brick::setToNormal(escript::Data& o
1038    
1039              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1040  #pragma omp for nowait  #pragma omp for nowait
1041                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1042                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1043                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1044                          // set vector at four quadrature points                          // set vector at four quadrature points
1045                          *o++ = 0.; *o++ = 1.; *o++ = 0.;                          *o++ = 0.; *o++ = 1.; *o++ = 0.;
# Line 1035  void Brick::setToNormal(escript::Data& o Line 1052  void Brick::setToNormal(escript::Data& o
1052    
1053              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1054  #pragma omp for nowait  #pragma omp for nowait
1055                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1056                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1057                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1058                          // set vector at four quadrature points                          // set vector at four quadrature points
1059                          *o++ = 0.; *o++ = 0.; *o++ = -1.;                          *o++ = 0.; *o++ = 0.; *o++ = -1.;
# Line 1049  void Brick::setToNormal(escript::Data& o Line 1066  void Brick::setToNormal(escript::Data& o
1066    
1067              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1068  #pragma omp for nowait  #pragma omp for nowait
1069                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1070                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1071                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1072                          // set vector at four quadrature points                          // set vector at four quadrature points
1073                          *o++ = 0.; *o++ = 0.; *o++ = 1.;                          *o++ = 0.; *o++ = 0.; *o++ = 1.;
# Line 1067  void Brick::setToNormal(escript::Data& o Line 1084  void Brick::setToNormal(escript::Data& o
1084          {          {
1085              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1086  #pragma omp for nowait  #pragma omp for nowait
1087                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1088                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1089                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1090                          *o++ = -1.;                          *o++ = -1.;
1091                          *o++ = 0.;                          *o++ = 0.;
# Line 1079  void Brick::setToNormal(escript::Data& o Line 1096  void Brick::setToNormal(escript::Data& o
1096    
1097              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1098  #pragma omp for nowait  #pragma omp for nowait
1099                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1100                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1101                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1102                          *o++ = 1.;                          *o++ = 1.;
1103                          *o++ = 0.;                          *o++ = 0.;
# Line 1091  void Brick::setToNormal(escript::Data& o Line 1108  void Brick::setToNormal(escript::Data& o
1108    
1109              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1110  #pragma omp for nowait  #pragma omp for nowait
1111                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1112                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1113                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1114                          *o++ = 0.;                          *o++ = 0.;
1115                          *o++ = -1.;                          *o++ = -1.;
# Line 1103  void Brick::setToNormal(escript::Data& o Line 1120  void Brick::setToNormal(escript::Data& o
1120    
1121              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1122  #pragma omp for nowait  #pragma omp for nowait
1123                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1124                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1125                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1126                          *o++ = 0.;                          *o++ = 0.;
1127                          *o++ = 1.;                          *o++ = 1.;
# Line 1115  void Brick::setToNormal(escript::Data& o Line 1132  void Brick::setToNormal(escript::Data& o
1132    
1133              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1134  #pragma omp for nowait  #pragma omp for nowait
1135                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1136                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1137                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1138                          *o++ = 0.;                          *o++ = 0.;
1139                          *o++ = 0.;                          *o++ = 0.;
# Line 1127  void Brick::setToNormal(escript::Data& o Line 1144  void Brick::setToNormal(escript::Data& o
1144    
1145              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1146  #pragma omp for nowait  #pragma omp for nowait
1147                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1148                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1149                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1150                          *o++ = 0.;                          *o++ = 0.;
1151                          *o++ = 0.;                          *o++ = 0.;
# Line 1153  void Brick::setToSize(escript::Data& out Line 1170  void Brick::setToSize(escript::Data& out
1170          out.requireWrite();          out.requireWrite();
1171          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
1172          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);          const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);
1173            const int NE = getNumElements();
1174  #pragma omp parallel for  #pragma omp parallel for
1175          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < NE; ++k) {
1176              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
1177              fill(o, o+numQuad, size);              fill(o, o+numQuad, size);
1178          }          }
# Line 1162  void Brick::setToSize(escript::Data& out Line 1180  void Brick::setToSize(escript::Data& out
1180              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {              || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1181          out.requireWrite();          out.requireWrite();
1182          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
1183            const int NE0 = m_NE[0];
1184            const int NE1 = m_NE[1];
1185            const int NE2 = m_NE[2];
1186  #pragma omp parallel  #pragma omp parallel
1187          {          {
1188              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1189                  const double size=min(m_dx[1],m_dx[2]);                  const double size=min(m_dx[1],m_dx[2]);
1190  #pragma omp for nowait  #pragma omp for nowait
1191                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1192                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1193                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
1194                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1195                      }                      }
# Line 1178  void Brick::setToSize(escript::Data& out Line 1199  void Brick::setToSize(escript::Data& out
1199              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1200                  const double size=min(m_dx[1],m_dx[2]);                  const double size=min(m_dx[1],m_dx[2]);
1201  #pragma omp for nowait  #pragma omp for nowait
1202                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1203                      for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                      for (index_t k1 = 0; k1 < NE1; ++k1) {
1204                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
1205                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1206                      }                      }
# Line 1189  void Brick::setToSize(escript::Data& out Line 1210  void Brick::setToSize(escript::Data& out
1210              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1211                  const double size=min(m_dx[0],m_dx[2]);                  const double size=min(m_dx[0],m_dx[2]);
1212  #pragma omp for nowait  #pragma omp for nowait
1213                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1214                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1215                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
1216                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1217                      }                      }
# Line 1200  void Brick::setToSize(escript::Data& out Line 1221  void Brick::setToSize(escript::Data& out
1221              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1222                  const double size=min(m_dx[0],m_dx[2]);                  const double size=min(m_dx[0],m_dx[2]);
1223  #pragma omp for nowait  #pragma omp for nowait
1224                  for (index_t k2 = 0; k2 < m_NE[2]; ++k2) {                  for (index_t k2 = 0; k2 < NE2; ++k2) {
1225                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1226                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
1227                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1228                      }                      }
# Line 1211  void Brick::setToSize(escript::Data& out Line 1232  void Brick::setToSize(escript::Data& out
1232              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1233                  const double size=min(m_dx[0],m_dx[1]);                  const double size=min(m_dx[0],m_dx[1]);
1234  #pragma omp for nowait  #pragma omp for nowait
1235                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1236                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1237                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
1238                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1239                      }                      }
# Line 1222  void Brick::setToSize(escript::Data& out Line 1243  void Brick::setToSize(escript::Data& out
1243              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1244                  const double size=min(m_dx[0],m_dx[1]);                  const double size=min(m_dx[0],m_dx[1]);
1245  #pragma omp for nowait  #pragma omp for nowait
1246                  for (index_t k1 = 0; k1 < m_NE[1]; ++k1) {                  for (index_t k1 = 0; k1 < NE1; ++k1) {
1247                      for (index_t k0 = 0; k0 < m_NE[0]; ++k0) {                      for (index_t k0 = 0; k0 < NE0; ++k0) {
1248                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
1249                          fill(o, o+numQuad, size);                          fill(o, o+numQuad, size);
1250                      }                      }
# Line 1266  void Brick::assembleCoordinates(escript: Line 1287  void Brick::assembleCoordinates(escript:
1287      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
1288          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
1289    
1290        const int NN0 = m_NN[0];
1291        const int NN1 = m_NN[1];
1292        const int NN2 = m_NN[2];
1293      arg.requireWrite();      arg.requireWrite();
1294  #pragma omp parallel for  #pragma omp parallel for
1295      for (dim_t i2 = 0; i2 < m_NN[2]; i2++) {      for (dim_t i2 = 0; i2 < NN2; i2++) {
1296          for (dim_t i1 = 0; i1 < m_NN[1]; i1++) {          for (dim_t i1 = 0; i1 < NN1; i1++) {
1297              for (dim_t i0 = 0; i0 < m_NN[0]; i0++) {              for (dim_t i0 = 0; i0 < NN0; i0++) {
1298                  double* point = arg.getSampleDataRW(i0+m_NN[0]*i1+m_NN[0]*m_NN[1]*i2);                  double* point = arg.getSampleDataRW(i0+NN0*i1+NN0*NN1*i2);
1299                  point[0] = getLocalCoordinate(i0, 0);                  point[0] = getLocalCoordinate(i0, 0);
1300                  point[1] = getLocalCoordinate(i1, 1);                  point[1] = getLocalCoordinate(i1, 1);
1301                  point[2] = getLocalCoordinate(i2, 2);                  point[2] = getLocalCoordinate(i2, 2);
# Line 1291  void Brick::assembleGradient(escript::Da Line 1315  void Brick::assembleGradient(escript::Da
1315      const double C4 = .5;      const double C4 = .5;
1316      const double C5 = .62200846792814621559;      const double C5 = .62200846792814621559;
1317      const double C6 = .78867513459481288225;      const double C6 = .78867513459481288225;
1318        const int NE0 = m_NE[0];
1319        const int NE1 = m_NE[1];
1320        const int NE2 = m_NE[2];
1321    
1322      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
1323          out.requireWrite();          out.requireWrite();
# Line 1305  void Brick::assembleGradient(escript::Da Line 1332  void Brick::assembleGradient(escript::Da
1332              vector<double> f_110(numComp);              vector<double> f_110(numComp);
1333              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1334  #pragma omp for  #pragma omp for
1335              for (index_t k2=0; k2 < m_NE[2]; ++k2) {              for (index_t k2=0; k2 < NE2; ++k2) {
1336                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1337                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1338                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1339                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1340                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1316  void Brick::assembleGradient(escript::Da Line 1343  void Brick::assembleGradient(escript::Da
1343                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1344                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1345                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1346                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,NE0,NE1));
1347                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1348                              const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / m_dx[0];
1349                              const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];                              const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / m_dx[0];
# Line 1372  void Brick::assembleGradient(escript::Da Line 1399  void Brick::assembleGradient(escript::Da
1399              vector<double> f_110(numComp);              vector<double> f_110(numComp);
1400              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1401  #pragma omp for  #pragma omp for
1402              for (index_t k2=0; k2 < m_NE[2]; ++k2) {              for (index_t k2=0; k2 < NE2; ++k2) {
1403                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1404                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1405                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1406                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1407                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1383  void Brick::assembleGradient(escript::Da Line 1410  void Brick::assembleGradient(escript::Da
1410                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1411                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1412                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1413                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE[0],m_NE[1]));                          double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,NE0,NE1));
1414                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1415                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1416                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1407  void Brick::assembleGradient(escript::Da Line 1434  void Brick::assembleGradient(escript::Da
1434              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1435              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1436  #pragma omp for nowait  #pragma omp for nowait
1437                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1438                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1439                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1440                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1441                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1417  void Brick::assembleGradient(escript::Da Line 1444  void Brick::assembleGradient(escript::Da
1444                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1445                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1446                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1447                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,NE1));
1448                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1449                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];                              const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / m_dx[1];
1450                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];                              const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / m_dx[1];
# Line 1441  void Brick::assembleGradient(escript::Da Line 1468  void Brick::assembleGradient(escript::Da
1468              } // end of face 0              } // end of face 0
1469              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1470  #pragma omp for nowait  #pragma omp for nowait
1471                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1472                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1473                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1474                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1475                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1451  void Brick::assembleGradient(escript::Da Line 1478  void Brick::assembleGradient(escript::Da
1478                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1479                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1480                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1481                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,NE1));
1482                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1483                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];                              const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / m_dx[1];
1484                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];                              const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / m_dx[1];
# Line 1475  void Brick::assembleGradient(escript::Da Line 1502  void Brick::assembleGradient(escript::Da
1502              } // end of face 1              } // end of face 1
1503              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1504  #pragma omp for nowait  #pragma omp for nowait
1505                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1506                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1507                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1508                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1509                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1485  void Brick::assembleGradient(escript::Da Line 1512  void Brick::assembleGradient(escript::Da
1512                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1513                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1514                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1515                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,NE0));
1516                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1517                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / m_dx[0];
1518                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];                              const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / m_dx[2];
# Line 1508  void Brick::assembleGradient(escript::Da Line 1535  void Brick::assembleGradient(escript::Da
1535              } // end of face 2              } // end of face 2
1536              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1537  #pragma omp for nowait  #pragma omp for nowait
1538                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1539                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1540                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1541                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1542                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1518  void Brick::assembleGradient(escript::Da Line 1545  void Brick::assembleGradient(escript::Da
1545                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1546                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1547                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1548                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,NE0));
1549                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1550                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];                              const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1551                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];                              const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
# Line 1542  void Brick::assembleGradient(escript::Da Line 1569  void Brick::assembleGradient(escript::Da
1569              } // end of face 3              } // end of face 3
1570              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1571  #pragma omp for nowait  #pragma omp for nowait
1572                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1573                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1574                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1575                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1576                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1552  void Brick::assembleGradient(escript::Da Line 1579  void Brick::assembleGradient(escript::Da
1579                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1580                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1581                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1582                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,NE0));
1583                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1584                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];                              const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / m_dx[0];
1585                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];                              const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / m_dx[0];
# Line 1576  void Brick::assembleGradient(escript::Da Line 1603  void Brick::assembleGradient(escript::Da
1603              } // end of face 4              } // end of face 4
1604              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1605  #pragma omp for nowait  #pragma omp for nowait
1606                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1607                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1608                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1609                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1610                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1586  void Brick::assembleGradient(escript::Da Line 1613  void Brick::assembleGradient(escript::Da
1613                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1614                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1615                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1616                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,NE0));
1617                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1618                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];                              const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / m_dx[0];
1619                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];                              const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / m_dx[0];
# Line 1623  void Brick::assembleGradient(escript::Da Line 1650  void Brick::assembleGradient(escript::Da
1650              vector<double> f_111(numComp);              vector<double> f_111(numComp);
1651              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1652  #pragma omp for nowait  #pragma omp for nowait
1653                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1654                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1655                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(0,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1656                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(0,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1657                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(0,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1633  void Brick::assembleGradient(escript::Da Line 1660  void Brick::assembleGradient(escript::Da
1660                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1661                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1662                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1663                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,NE1));
1664                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1665                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1666                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / m_dx[1];
# Line 1644  void Brick::assembleGradient(escript::Da Line 1671  void Brick::assembleGradient(escript::Da
1671              } // end of face 0              } // end of face 0
1672              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1673  #pragma omp for nowait  #pragma omp for nowait
1674                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1675                      for (index_t k1=0; k1 < m_NE[1]; ++k1) {                      for (index_t k1=0; k1 < NE1; ++k1) {
1676                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1677                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1678                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(m_NN[0]-2,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1654  void Brick::assembleGradient(escript::Da Line 1681  void Brick::assembleGradient(escript::Da
1681                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1682                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1683                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(m_NN[0]-1,k1+1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1684                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));                          double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,NE1));
1685                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1686                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / m_dx[0];
1687                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / m_dx[1];
# Line 1665  void Brick::assembleGradient(escript::Da Line 1692  void Brick::assembleGradient(escript::Da
1692              } // end of face 1              } // end of face 1
1693              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1694  #pragma omp for nowait  #pragma omp for nowait
1695                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1696                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1697                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,0,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1698                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1699                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1675  void Brick::assembleGradient(escript::Da Line 1702  void Brick::assembleGradient(escript::Da
1702                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1703                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1704                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1705                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,NE0));
1706                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1707                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / m_dx[0];
1708                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1686  void Brick::assembleGradient(escript::Da Line 1713  void Brick::assembleGradient(escript::Da
1713              } // end of face 2              } // end of face 2
1714              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1715  #pragma omp for nowait  #pragma omp for nowait
1716                  for (index_t k2=0; k2 < m_NE[2]; ++k2) {                  for (index_t k2=0; k2 < NE2; ++k2) {
1717                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1718                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1719                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1720                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1696  void Brick::assembleGradient(escript::Da Line 1723  void Brick::assembleGradient(escript::Da
1723                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-2,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1724                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1725                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,m_NN[1]-1,k2+1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1726                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,NE0));
1727                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1728                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / m_dx[0];
1729                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / m_dx[1];
# Line 1707  void Brick::assembleGradient(escript::Da Line 1734  void Brick::assembleGradient(escript::Da
1734              } // end of face 3              } // end of face 3
1735              if (m_faceOffset[4] > -1) {              if (m_faceOffset[4] > -1) {
1736  #pragma omp for nowait  #pragma omp for nowait
1737                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1738                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1739                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1740                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1741                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1717  void Brick::assembleGradient(escript::Da Line 1744  void Brick::assembleGradient(escript::Da
1744                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1745                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_NN[0],m_NN[1])), numComp*sizeof(double));
1746                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1747                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,NE0));
1748                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1749                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / m_dx[0];
1750                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / m_dx[1];
# Line 1728  void Brick::assembleGradient(escript::Da Line 1755  void Brick::assembleGradient(escript::Da
1755              } // end of face 4              } // end of face 4
1756              if (m_faceOffset[5] > -1) {              if (m_faceOffset[5] > -1) {
1757  #pragma omp for nowait  #pragma omp for nowait
1758                  for (index_t k1=0; k1 < m_NE[1]; ++k1) {                  for (index_t k1=0; k1 < NE1; ++k1) {
1759                      for (index_t k0=0; k0 < m_NE[0]; ++k0) {                      for (index_t k0=0; k0 < NE0; ++k0) {
1760                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_000[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1761                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_001[0], in.getSampleDataRO(INDEX3(k0,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1762                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_010[0], in.getSampleDataRO(INDEX3(k0,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
# Line 1738  void Brick::assembleGradient(escript::Da Line 1765  void Brick::assembleGradient(escript::Da
1765                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_101[0], in.getSampleDataRO(INDEX3(k0+1,k1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1766                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_110[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-2, m_NN[0],m_NN[1])), numComp*sizeof(double));
1767                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));                          memcpy(&f_111[0], in.getSampleDataRO(INDEX3(k0+1,k1+1,m_NN[2]-1, m_NN[0],m_NN[1])), numComp*sizeof(double));
1768                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));                          double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,NE0));
1769                          for (index_t i=0; i < numComp; ++i) {                          for (index_t i=0; i < numComp; ++i) {
1770                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];                              o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / m_dx[0];
1771                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];                              o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / m_dx[1];
# Line 2056  void Brick::dofToNodes(escript::Data& ou Line 2083  void Brick::dofToNodes(escript::Data& ou
2083      coupler->startCollect(in.getDataRO());      coupler->startCollect(in.getDataRO());
2084    
2085      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
2086        const dim_t numNodes = getNumNodes();
2087      out.requireWrite();      out.requireWrite();
2088      const double* buffer = coupler->finishCollect();      const double* buffer = coupler->finishCollect();
2089    
2090  #pragma omp parallel for  #pragma omp parallel for
2091      for (index_t i=0; i<getNumNodes(); i++) {      for (index_t i=0; i<numNodes; i++) {
2092          const double* src=(m_dofMap[i]<numDOF ?          const double* src=(m_dofMap[i]<numDOF ?
2093                  in.getSampleDataRO(m_dofMap[i])                  in.getSampleDataRO(m_dofMap[i])
2094                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
# Line 2127  void Brick::populateSampleIds() Line 2155  void Brick::populateSampleIds()
2155      else      else
2156          m_faceCount[5]=0;          m_faceCount[5]=0;
2157    
2158      m_faceId.resize(getNumFaceElements());      const dim_t NFE = getNumFaceElements();
2159        m_faceId.resize(NFE);
2160    
2161      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
2162      const index_t bottom = (m_offset[1]==0 ? 0 : 1);      const index_t bottom = (m_offset[1]==0 ? 0 : 1);
# Line 2135  void Brick::populateSampleIds() Line 2164  void Brick::populateSampleIds()
2164      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];      const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
2165      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];      const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
2166      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];      const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
2167        const dim_t NN0 = m_NN[0];
2168        const dim_t NN1 = m_NN[1];
2169        const dim_t NN2 = m_NN[2];
2170        const dim_t NE0 = m_NE[0];
2171        const dim_t NE1 = m_NE[1];
2172        const dim_t NE2 = m_NE[2];
2173    
2174      // the following is a compromise between efficiency and code length to      // the following is a compromise between efficiency and code length to
2175      // set the node id's according to the order mentioned above.      // set the node id's according to the order mentioned above.
# Line 2153  void Brick::populateSampleIds() Line 2188  void Brick::populateSampleIds()
2188          // set edge id's          // set edge id's
2189          // edges in x-direction, including corners          // edges in x-direction, including corners
2190  #pragma omp for nowait  #pragma omp for nowait
2191          for (dim_t i=0; i<m_NN[0]; i++) {          for (dim_t i=0; i<NN0; i++) {
2192              m_nodeId[i] = globalNodeId(i, 0, 0); // LF              m_nodeId[i] = globalNodeId(i, 0, 0); // LF
2193              m_nodeId[m_NN[0]*(m_NN[1]-1)+i] = globalNodeId(i, m_NN[1]-1, 0); // UF              m_nodeId[NN0*(NN1-1)+i] = globalNodeId(i, NN1-1, 0); // UF
2194              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+i] = globalNodeId(i, 0, m_NN[2]-1); // LB              m_nodeId[NN0*NN1*(NN2-1)+i] = globalNodeId(i, 0, NN2-1); // LB
2195              m_nodeId[m_NN[0]*m_NN[1]*m_NN[2]-m_NN[0]+i] = globalNodeId(i, m_NN[1]-1, m_NN[2]-1); // UB              m_nodeId[NN0*NN1*NN2-NN0+i] = globalNodeId(i, NN1-1, NN2-1); // UB
2196          }          }
2197          // edges in y-direction, without corners          // edges in y-direction, without corners
2198  #pragma omp for nowait  #pragma omp for nowait
2199          for (dim_t i=1; i<m_NN[1]-1; i++) {          for (dim_t i=1; i<NN1-1; i++) {
2200              m_nodeId[m_NN[0]*i] = globalNodeId(0, i, 0); // FL              m_nodeId[NN0*i] = globalNodeId(0, i, 0); // FL
2201              m_nodeId[m_NN[0]*(i+1)-1] = globalNodeId(m_NN[0]-1, i, 0); // FR              m_nodeId[NN0*(i+1)-1] = globalNodeId(NN0-1, i, 0); // FR
2202              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+m_NN[0]*i] = globalNodeId(0, i, m_NN[2]-1); // BL              m_nodeId[NN0*NN1*(NN2-1)+NN0*i] = globalNodeId(0, i, NN2-1); // BL
2203              m_nodeId[m_NN[0]*m_NN[1]*(m_NN[2]-1)+m_NN[0]*(i+1)-1] = globalNodeId(m_NN[0]-1, i, m_NN[2]-1); // BR              m_nodeId[NN0*NN1*(NN2-1)+NN0*(i+1)-1] = globalNodeId(NN0-1, i, NN2-1); // BR
2204          }          }
2205          // edges in z-direction, without corners          // edges in z-direction, without corners
2206  #pragma omp for  #pragma omp for
2207          for (dim_t i=1; i<m_NN[2]-1; i++) {          for (dim_t i=1; i<NN2-1; i++) {
2208              m_nodeId[m_NN[0]*m_NN[1]*i] = globalNodeId(0, 0, i); // LL              m_nodeId[NN0*NN1*i] = globalNodeId(0, 0, i); // LL
2209              m_nodeId[m_NN[0]*m_NN[1]*i+m_NN[0]-1] = globalNodeId(m_NN[0]-1, 0, i); // LR              m_nodeId[NN0*NN1*i+NN0-1] = globalNodeId(NN0-1, 0, i); // LR
2210              m_nodeId[m_NN[0]*m_NN[1]*(i+1)-m_NN[0]] = globalNodeId(0, m_NN[1]-1, i); // UL              m_nodeId[NN0*NN1*(i+1)-NN0] = globalNodeId(0, NN1-1, i); // UL
2211              m_nodeId[m_NN[0]*m_NN[1]*(i+1)-1] = globalNodeId(m_NN[0]-1, m_NN[1]-1, i); // UR              m_nodeId[NN0*NN1*(i+1)-1] = globalNodeId(NN0-1, NN1-1, i); // UR
2212          }          }
2213          // implicit barrier here because some node IDs will be overwritten          // implicit barrier here because some node IDs will be overwritten
2214          // below          // below
# Line 2183  void Brick::populateSampleIds() Line 2218  void Brick::populateSampleIds()
2218          for (dim_t i=0; i<nDOF2; i++) {          for (dim_t i=0; i<nDOF2; i++) {
2219              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2220                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2221                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=k+left+(j+bottom)*NN0+(i+front)*NN0*NN1;
2222                      const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;                      const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;
2223                      m_dofId[dofIdx] = m_nodeId[nodeIdx]                      m_dofId[dofIdx] = m_nodeId[nodeIdx]
2224                          = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;                          = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
# Line 2196  void Brick::populateSampleIds() Line 2231  void Brick::populateSampleIds()
2231  #pragma omp for nowait  #pragma omp for nowait
2232              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2233                  for (dim_t j=0; j<nDOF1; j++) {                  for (dim_t j=0; j<nDOF1; j++) {
2234                      const index_t nodeIdx=(j+bottom)*m_NN[0]+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=(j+bottom)*NN0+(i+front)*NN0*NN1;
2235                      const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;                      const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;
2236                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2237                          = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
# Line 2207  void Brick::populateSampleIds() Line 2242  void Brick::populateSampleIds()
2242  #pragma omp for nowait  #pragma omp for nowait
2243              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2244                  for (dim_t j=0; j<nDOF1; j++) {                  for (dim_t j=0; j<nDOF1; j++) {
2245                      const index_t nodeIdx=(j+bottom+1)*m_NN[0]-1+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=(j+bottom+1)*NN0-1+(i+front)*NN0*NN1;
2246                      const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;                      const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;
2247                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2248                          = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
# Line 2218  void Brick::populateSampleIds() Line 2253  void Brick::populateSampleIds()
2253  #pragma omp for nowait  #pragma omp for nowait
2254              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2255                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2256                      const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1];                      const index_t nodeIdx=k+left+(i+front)*NN0*NN1;
2257                      const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;                      const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;
2258                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2259                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
# Line 2229  void Brick::populateSampleIds() Line 2264  void Brick::populateSampleIds()
2264  #pragma omp for nowait  #pragma omp for nowait
2265              for (dim_t i=0; i<nDOF2; i++) {              for (dim_t i=0; i<nDOF2; i++) {
2266                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2267                      const index_t nodeIdx=k+left+(i+front)*m_NN[0]*m_NN[1]+m_NN[0]*(m_NN[1]-1);                      const index_t nodeIdx=k+left+(i+front)*NN0*NN1+NN0*(NN1-1);
2268                      const index_t dofId=k+i*nDOF0*nDOF1;                      const index_t dofId=k+i*nDOF0*nDOF1;
2269                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2270                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
# Line 2240  void Brick::populateSampleIds() Line 2275  void Brick::populateSampleIds()
2275  #pragma omp for nowait  #pragma omp for nowait
2276              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2277                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2278                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0];                      const index_t nodeIdx=k+left+(j+bottom)*NN0;
2279                      const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);                      const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);
2280                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2281                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;
# Line 2251  void Brick::populateSampleIds() Line 2286  void Brick::populateSampleIds()
2286  #pragma omp for nowait  #pragma omp for nowait
2287              for (dim_t j=0; j<nDOF1; j++) {              for (dim_t j=0; j<nDOF1; j++) {
2288                  for (dim_t k=0; k<nDOF0; k++) {                  for (dim_t k=0; k<nDOF0; k++) {
2289                      const index_t nodeIdx=k+left+(j+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1);                      const index_t nodeIdx=k+left+(j+bottom)*NN0+NN0*NN1*(NN2-1);
2290                      const index_t dofId=k+j*nDOF0;                      const index_t dofId=k+j*nDOF0;
2291                      m_nodeId[nodeIdx]                      m_nodeId[nodeIdx]
2292                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;                          = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;
# Line 2261  void Brick::populateSampleIds() Line 2296  void Brick::populateSampleIds()
2296    
2297          // populate element id's          // populate element id's
2298  #pragma omp for nowait  #pragma omp for nowait
2299          for (dim_t i2=0; i2<m_NE[2]; i2++) {          for (dim_t i2=0; i2<NE2; i2++) {
2300              for (dim_t i1=0; i1<m_NE[1]; i1++) {              for (dim_t i1=0; i1<NE1; i1++) {
2301                  for (dim_t i0=0; i0<m_NE[0]; i0++) {                  for (dim_t i0=0; i0<NE0; i0++) {
2302                      m_elementId[i0+i1*m_NE[0]+i2*m_NE[0]*m_NE[1]] =                      m_elementId[i0+i1*NE0+i2*NE0*NE1] =
2303                          (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]                          (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]
2304                          +(m_offset[1]+i1)*m_gNE[0]                          +(m_offset[1]+i1)*m_gNE[0]
2305                          +m_offset[0]+i0;                          +m_offset[0]+i0;
# Line 2274  void Brick::populateSampleIds() Line 2309  void Brick::populateSampleIds()
2309    
2310          // face elements          // face elements
2311  #pragma omp for  #pragma omp for
2312          for (dim_t k=0; k<getNumFaceElements(); k++)          for (dim_t k=0; k<NFE; k++)
2313              m_faceId[k]=k;              m_faceId[k]=k;
2314      } // end parallel section      } // end parallel section
2315    

Legend:
Removed from v.4988  
changed lines
  Added in v.5084

  ViewVC Help
Powered by ViewVC 1.1.26