/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp revision 3791 by caltinay, Wed Feb 1 05:10:22 2012 UTC trunk/ripley/src/Rectangle.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/Rectangle.h>  #include <ripley/Rectangle.h>
17  extern "C" {  extern "C" {
# Line 32  namespace ripley { Line 34  namespace ripley {
34  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,  Rectangle::Rectangle(int n0, int n1, double x0, double y0, double x1,
35                       double y1, int d0, int d1) :                       double y1, int d0, int d1) :
36      RipleyDomain(2),      RipleyDomain(2),
     m_gNE0(n0),  
     m_gNE1(n1),  
37      m_x0(x0),      m_x0(x0),
38      m_y0(y0),      m_y0(y0),
39      m_l0(x1-x0),      m_l0(x1-x0),
40      m_l1(y1-y0),      m_l1(y1-y0)
     m_NX(d0),  
     m_NY(d1)  
41  {  {
42        // ignore subdivision parameters for serial run
43        if (m_mpiInfo->size == 1) {
44            d0=1;
45            d1=1;
46        }
47    
48        bool warn=false;
49        // if number of subdivisions is non-positive, try to subdivide by the same
50        // ratio as the number of elements
51        if (d0<=0 && d1<=0) {
52            warn=true;
53            d0=(int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));
54            d1=m_mpiInfo->size/d0;
55            if (d0*d1 != m_mpiInfo->size) {
56                // ratios not the same so subdivide side with more elements only
57                if (n0>n1) {
58                    d0=0;
59                    d1=1;
60                } else {
61                    d0=1;
62                    d1=0;
63                }
64            }
65        }
66        if (d0<=0) {
67            // d1 is preset, determine d0 - throw further down if result is no good
68            d0=m_mpiInfo->size/d1;
69        } else if (d1<=0) {
70            // d0 is preset, determine d1 - throw further down if result is no good
71            d1=m_mpiInfo->size/d0;
72        }
73    
74        m_NX=d0;
75        m_NY=d1;
76    
77      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
78      // among number of ranks      // among number of ranks
79      if (m_NX*m_NY != m_mpiInfo->size)      if (m_NX*m_NY != m_mpiInfo->size)
80          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
81    
82      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0)      if (warn) {
83          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
84                << d1 << "). This may not be optimal!" << endl;
85        }
86    
87        if ((n0+1)%m_NX > 0) {
88            double Dx=m_l0/n0;
89            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
90            m_l0=Dx*n0;
91            cout << "Warning: Adjusted number of elements and length. N0="
92                << n0 << ", l0=" << m_l0 << endl;
93        }
94        if ((n1+1)%m_NY > 0) {
95            double Dy=m_l1/n1;
96            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
97            m_l1=Dy*n1;
98            cout << "Warning: Adjusted number of elements and length. N1="
99                << n1 << ", l1=" << m_l1 << endl;
100        }
101    
102        m_gNE0=n0;
103        m_gNE1=n1;
104    
105      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2))
106          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");
# Line 106  bool Rectangle::operator==(const Abstrac Line 159  bool Rectangle::operator==(const Abstrac
159      return false;      return false;
160  }  }
161    
162    void Rectangle::readBinaryGrid(escript::Data& out, string filename,
163                const vector<int>& first, const vector<int>& numValues) const
164    {
165        // check destination function space
166        int myN0, myN1;
167        if (out.getFunctionSpace().getTypeCode() == Nodes) {
168            myN0 = m_N0;
169            myN1 = m_N1;
170        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
171                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
172            myN0 = m_NE0;
173            myN1 = m_NE1;
174        } else
175            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
176    
177        // check file existence and size
178        ifstream f(filename.c_str(), ifstream::binary);
179        if (f.fail()) {
180            throw RipleyException("readBinaryGrid(): cannot open file");
181        }
182        f.seekg(0, ios::end);
183        const int numComp = out.getDataPointSize();
184        const int filesize = f.tellg();
185        const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(float);
186        if (filesize < reqsize) {
187            f.close();
188            throw RipleyException("readBinaryGrid(): not enough data in file");
189        }
190    
191        // check if this rank contributes anything
192        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
193                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1) {
194            f.close();
195            return;
196        }
197    
198        // now determine how much this rank has to write
199    
200        // first coordinates in data object to write to
201        const int first0 = max(0, first[0]-m_offset0);
202        const int first1 = max(0, first[1]-m_offset1);
203        // indices to first value in file
204        const int idx0 = max(0, m_offset0-first[0]);
205        const int idx1 = max(0, m_offset1-first[1]);
206        // number of values to write
207        const int num0 = min(numValues[0]-idx0, myN0-first0);
208        const int num1 = min(numValues[1]-idx1, myN1-first1);
209    
210        out.requireWrite();
211        vector<float> values(num0*numComp);
212        const int dpp = out.getNumDataPointsPerSample();
213    
214        for (index_t y=0; y<num1; y++) {
215            const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);
216            f.seekg(fileofs*sizeof(float));
217            f.read((char*)&values[0], num0*numComp*sizeof(float));
218            for (index_t x=0; x<num0; x++) {
219                double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0);
220                for (index_t c=0; c<numComp; c++) {
221                    for (index_t q=0; q<dpp; q++) {
222                        *dest++ = static_cast<double>(values[x*numComp+c]);
223                    }
224                }
225            }
226        }
227    
228        f.close();
229    }
230    
231  void Rectangle::dump(const string& fileName) const  void Rectangle::dump(const string& fileName) const
232  {  {
233  #if USE_SILO  #if USE_SILO
# Line 114  void Rectangle::dump(const string& fileN Line 236  void Rectangle::dump(const string& fileN
236          fn+=".silo";          fn+=".silo";
237      }      }
238    
     const int NUM_SILO_FILES = 1;  
     const char* blockDirFmt = "/block%04d";  
239      int driver=DB_HDF5;          int driver=DB_HDF5;    
240      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
241        const char* blockDirFmt = "/block%04d";
242    
243  #ifdef ESYS_MPI  #ifdef ESYS_MPI
244      PMPIO_baton_t* baton = NULL;      PMPIO_baton_t* baton = NULL;
245        const int NUM_SILO_FILES = 1;
246  #endif  #endif
247    
248      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 434  void Rectangle::setToSize(escript::Data& Line 556  void Rectangle::setToSize(escript::Data&
556          const dim_t numQuad=out.getNumDataPointsPerSample();          const dim_t numQuad=out.getNumDataPointsPerSample();
557          const double hSize=getFirstCoordAndSpacing(0).second;          const double hSize=getFirstCoordAndSpacing(0).second;
558          const double vSize=getFirstCoordAndSpacing(1).second;          const double vSize=getFirstCoordAndSpacing(1).second;
559          const double size=min(hSize,vSize);          const double size=sqrt(hSize*hSize+vSize*vSize);
560  #pragma omp parallel for  #pragma omp parallel for
561          for (index_t k = 0; k < getNumElements(); ++k) {          for (index_t k = 0; k < getNumElements(); ++k) {
562              double* o = out.getSampleDataRW(k);              double* o = out.getSampleDataRW(k);
# Line 632  void Rectangle::assembleGradient(escript Line 754  void Rectangle::assembleGradient(escript
754    
755      if (out.getFunctionSpace().getTypeCode() == Elements) {      if (out.getFunctionSpace().getTypeCode() == Elements) {
756          out.requireWrite();          out.requireWrite();
757  #pragma omp parallel for  #pragma omp parallel
758          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
759              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
760                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
761                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
762                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
763                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
764                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
765                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
766                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
767                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
768                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
769                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
770                      o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
771                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      for (index_t i=0; i < numComp; ++i) {
772                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
773                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
774                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
775              } /* end of k0 loop */                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
776          } /* end of k1 loop */                          o[INDEX3(i,0,2,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
777                            o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
778                            o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
779                            o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
780                        } // end of component loop i
781                    } // end of k0 loop
782                } // end of k1 loop
783            } // end of parallel section
784      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
785          out.requireWrite();          out.requireWrite();
786  #pragma omp parallel for  #pragma omp parallel
787          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
788              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
789                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
790                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
791                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
792                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
793                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
794                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
795                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
796                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
797                  } /* end of component loop i */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
798              } /* end of k0 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
799          } /* end of k1 loop */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
800                        for (index_t i=0; i < numComp; ++i) {
801                            o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
802                            o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
803                        } // end of component loop i
804                    } // end of k0 loop
805                } // end of k1 loop
806            } // end of parallel section
807      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
808          out.requireWrite();          out.requireWrite();
809  #pragma omp parallel  #pragma omp parallel
810          {          {
811                vector<double> f_00(numComp);
812                vector<double> f_01(numComp);
813                vector<double> f_10(numComp);
814                vector<double> f_11(numComp);
815              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
816  #pragma omp for nowait  #pragma omp for nowait
817                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
818                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
819                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
820                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
821                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
822                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
823                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
824                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
825                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
826                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
827                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
828                      } /* end of component loop i */                      } // end of component loop i
829                  } /* end of k1 loop */                  } // end of k1 loop
830              } /* end of face 0 */              } // end of face 0
831              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
832  #pragma omp for nowait  #pragma omp for nowait
833                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
834                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
835                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
836                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
837                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
838                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
839                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
840                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx1 + f_01[i]*cx3 + f_10[i]*cx6 + f_11[i]*cx4;
841                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
842                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
843                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
844                      } /* end of component loop i */                      } // end of component loop i
845                  } /* end of k1 loop */                  } // end of k1 loop
846              } /* end of face 1 */              } // end of face 1
847              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
848  #pragma omp for nowait  #pragma omp for nowait
849                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
850                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
851                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
852                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
853                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
854                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
855                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
856                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
857                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
858                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
859                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
860                      } /* end of component loop i */                      } // end of component loop i
861                  } /* end of k0 loop */                  } // end of k0 loop
862              } /* end of face 2 */              } // end of face 2
863              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
864  #pragma omp for nowait  #pragma omp for nowait
865                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
866                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
867                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
868                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
869                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
870                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
871                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
872                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
873                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
874                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
875                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
876                      } /* end of component loop i */                      } // end of component loop i
877                  } /* end of k0 loop */                  } // end of k0 loop
878              } /* end of face 3 */              } // end of face 3
879          } // end of parallel section          } // end of parallel section
880    
881      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
882          out.requireWrite();          out.requireWrite();
883  #pragma omp parallel  #pragma omp parallel
884          {          {
885                vector<double> f_00(numComp);
886                vector<double> f_01(numComp);
887                vector<double> f_10(numComp);
888                vector<double> f_11(numComp);
889              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
890  #pragma omp for nowait  #pragma omp for nowait
891                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
892                      const double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
893                      const double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
894                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(1,k1, m_N0)), numComp*sizeof(double));
895                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(1,k1+1, m_N0)), numComp*sizeof(double));
896                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
897                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
898                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
899                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
900                      } /* end of component loop i */                      } // end of component loop i
901                  } /* end of k1 loop */                  } // end of k1 loop
902              } /* end of face 0 */              } // end of face 0
903              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
904  #pragma omp for nowait  #pragma omp for nowait
905                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
906                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0)), numComp*sizeof(double));
907                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0)), numComp*sizeof(double));
908                      const double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
909                      const double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
910                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
911                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
912                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
913                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
914                      } /* end of component loop i */                      } // end of component loop i
915                  } /* end of k1 loop */                  } // end of k1 loop
916              } /* end of face 1 */              } // end of face 1
917              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
918  #pragma omp for nowait  #pragma omp for nowait
919                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
920                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
921                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,1, m_N0)), numComp*sizeof(double));
922                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
923                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,1, m_N0)), numComp*sizeof(double));
924                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
925                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
926                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
927                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
928                      } /* end of component loop i */                      } // end of component loop i
929                  } /* end of k0 loop */                  } // end of k0 loop
930              } /* end of face 2 */              } // end of face 2
931              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
932  #pragma omp for nowait  #pragma omp for nowait
933                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
934                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0)), numComp*sizeof(double));
935                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
936                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0)), numComp*sizeof(double));
937                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
938                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
939                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
940                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
941                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
942                      } /* end of component loop i */                      } // end of component loop i
943                  } /* end of k0 loop */                  } // end of k0 loop
944              } /* end of face 3 */              } // end of face 3
945          } // end of parallel section          } // end of parallel section
946      }      }
947  }  }
# Line 809  void Rectangle::assembleIntegrate(vector Line 954  void Rectangle::assembleIntegrate(vector
954      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
955      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
956      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
957      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs=arg.getFunctionSpace().getTypeCode();
958          const double w = h0*h1/4.;      if (fs == Elements && arg.actsExpanded()) {
959  #pragma omp parallel  #pragma omp parallel
960          {          {
961              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
962                const double w = h0*h1/4.;
963  #pragma omp for nowait  #pragma omp for nowait
964              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
965                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
# Line 824  void Rectangle::assembleIntegrate(vector Line 970  void Rectangle::assembleIntegrate(vector
970                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
971                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
972                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
973                      }  /* end of component loop i */                      }  // end of component loop i
974                  } /* end of k0 loop */                  } // end of k0 loop
975              } /* end of k1 loop */              } // end of k1 loop
   
976  #pragma omp critical  #pragma omp critical
977              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
978                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
979          } // end of parallel section          } // end of parallel section
980      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
981        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
982          const double w = h0*h1;          const double w = h0*h1;
983  #pragma omp parallel  #pragma omp parallel
984          {          {
# Line 843  void Rectangle::assembleIntegrate(vector Line 989  void Rectangle::assembleIntegrate(vector
989                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
990                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
991                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
992                      }  /* end of component loop i */                      }
993                  } /* end of k0 loop */                  }
994              } /* end of k1 loop */              }
   
995  #pragma omp critical  #pragma omp critical
996              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
997                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
998          } // end of parallel section          } // end of parallel section
999      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
1000          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
1001  #pragma omp parallel  #pragma omp parallel
1002          {          {
1003              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
1004                const double w0 = h0/2.;
1005                const double w1 = h1/2.;
1006              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1007  #pragma omp for nowait  #pragma omp for nowait
1008                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
# Line 865  void Rectangle::assembleIntegrate(vector Line 1011  void Rectangle::assembleIntegrate(vector
1011                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1012                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1013                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1014                      }  /* end of component loop i */                      }  // end of component loop i
1015                  } /* end of k1 loop */                  } // end of k1 loop
1016              }              }
1017    
1018              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 877  void Rectangle::assembleIntegrate(vector Line 1023  void Rectangle::assembleIntegrate(vector
1023                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1024                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1025                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
1026                      }  /* end of component loop i */                      }  // end of component loop i
1027                  } /* end of k1 loop */                  } // end of k1 loop
1028              }              }
1029    
1030              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 889  void Rectangle::assembleIntegrate(vector Line 1035  void Rectangle::assembleIntegrate(vector
1035                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1036                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1037                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1038                      }  /* end of component loop i */                      }  // end of component loop i
1039                  } /* end of k0 loop */                  } // end of k0 loop
1040              }              }
1041    
1042              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 901  void Rectangle::assembleIntegrate(vector Line 1047  void Rectangle::assembleIntegrate(vector
1047                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
1048                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
1049                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
1050                      }  /* end of component loop i */                      }  // end of component loop i
1051                  } /* end of k0 loop */                  } // end of k0 loop
1052              }              }
   
1053  #pragma omp critical  #pragma omp critical
1054              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1055                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1056          } // end of parallel section          } // end of parallel section
1057      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
1058        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
1059  #pragma omp parallel  #pragma omp parallel
1060          {          {
1061              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 919  void Rectangle::assembleIntegrate(vector Line 1065  void Rectangle::assembleIntegrate(vector
1065                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
1066                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1067                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
1068                      }  /* end of component loop i */                      }
1069                  } /* end of k1 loop */                  }
1070              }              }
1071    
1072              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 929  void Rectangle::assembleIntegrate(vector Line 1075  void Rectangle::assembleIntegrate(vector
1075                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
1076                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1077                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
1078                      }  /* end of component loop i */                      }
1079                  } /* end of k1 loop */                  }
1080              }              }
1081    
1082              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 939  void Rectangle::assembleIntegrate(vector Line 1085  void Rectangle::assembleIntegrate(vector
1085                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
1086                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1087                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
1088                      }  /* end of component loop i */                      }
1089                  } /* end of k0 loop */                  }
1090              }              }
1091    
1092              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 949  void Rectangle::assembleIntegrate(vector Line 1095  void Rectangle::assembleIntegrate(vector
1095                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
1096                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1097                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
1098                      }  /* end of component loop i */                      }
1099                  } /* end of k0 loop */                  }
1100              }              }
1101    
1102  #pragma omp critical  #pragma omp critical
1103              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
1104                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
1105          } // end of parallel section          } // end of parallel section
1106      }      } // function space selector
1107  }  }
1108    
1109  //protected  //protected
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1455  void Rectangle::interpolateNodesOnElemen
1455      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1456      if (reduced) {      if (reduced) {
1457          out.requireWrite();          out.requireWrite();
1458          const double c0 = .25;          const double c0 = 0.25;
1459  #pragma omp parallel for  #pragma omp parallel
1460          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1461              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1462                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1463                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1464                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1465                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1466                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1467                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1468                      o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
1469                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1470              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1471          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1472                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1473                        for (index_t i=0; i < numComp; ++i) {
1474                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1475                        } /* end of component loop i */
1476                    } /* end of k0 loop */
1477                } /* end of k1 loop */
1478            } /* end of parallel section */
1479      } else {      } else {
1480          out.requireWrite();          out.requireWrite();
1481          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1482          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1483          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1484  #pragma omp parallel for  #pragma omp parallel
1485          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1486              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1487                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1488                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1489                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1490                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1491                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1492                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1493                      o[INDEX2(i,numComp,0)] = f_00[i]*c2 + f_11[i]*c1 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,k1, m_N0)), numComp*sizeof(double));
1494                      o[INDEX2(i,numComp,1)] = f_01[i]*c1 + f_10[i]*c2 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1495                      o[INDEX2(i,numComp,2)] = f_01[i]*c2 + f_10[i]*c1 + c0*(f_00[i] + f_11[i]);                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1496                      o[INDEX2(i,numComp,3)] = f_00[i]*c1 + f_11[i]*c2 + c0*(f_01[i] + f_10[i]);                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1497                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1498              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1499          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1500                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1501                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1502                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1503                        } /* end of component loop i */
1504                    } /* end of k0 loop */
1505                } /* end of k1 loop */
1506            } /* end of parallel section */
1507      }      }
1508  }  }
1509    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1514  void Rectangle::interpolateNodesOnFaces(
1514      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1515      if (reduced) {      if (reduced) {
1516          out.requireWrite();          out.requireWrite();
1517          const double c0 = .5;          const double c0 = 0.5;
1518  #pragma omp parallel  #pragma omp parallel
1519          {          {
1520                vector<double> f_00(numComp);
1521                vector<double> f_01(numComp);
1522                vector<double> f_10(numComp);
1523                vector<double> f_11(numComp);
1524              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1525  #pragma omp for nowait  #pragma omp for nowait
1526                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1527                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1528                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1529                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1530                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1531                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i]);
# Line 1371  void Rectangle::interpolateNodesOnFaces( Line 1535  void Rectangle::interpolateNodesOnFaces(
1535              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1536  #pragma omp for nowait  #pragma omp for nowait
1537                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1538                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1539                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1540                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1541                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1542                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_10[i] + f_11[i]);
# Line 1382  void Rectangle::interpolateNodesOnFaces( Line 1546  void Rectangle::interpolateNodesOnFaces(
1546              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1547  #pragma omp for nowait  #pragma omp for nowait
1548                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1549                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1550                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1551                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1552                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1553                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_10[i]);
# Line 1393  void Rectangle::interpolateNodesOnFaces( Line 1557  void Rectangle::interpolateNodesOnFaces(
1557              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1558  #pragma omp for nowait  #pragma omp for nowait
1559                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1560                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1561                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1562                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1563                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1564                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1565                      } /* end of component loop i */                      } /* end of component loop i */
1566                  } /* end of k0 loop */                  } /* end of k0 loop */
1567              } /* end of face 3 */              } /* end of face 3 */
1568          } // end of parallel section          } /* end of parallel section */
1569      } else {      } else {
1570          out.requireWrite();          out.requireWrite();
1571          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1572          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1573  #pragma omp parallel  #pragma omp parallel
1574          {          {
1575                vector<double> f_00(numComp);
1576                vector<double> f_01(numComp);
1577                vector<double> f_10(numComp);
1578                vector<double> f_11(numComp);
1579              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1580  #pragma omp for nowait  #pragma omp for nowait
1581                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1582                      const double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(0,k1, m_N0)), numComp*sizeof(double));
1583                      const double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(0,k1+1, m_N0)), numComp*sizeof(double));
1584                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1585                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1586                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_01[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_01[i] + c1*f_00[i];
1587                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_01[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_01[i];
1588                      } /* end of component loop i */                      } /* end of component loop i */
1589                  } /* end of k1 loop */                  } /* end of k1 loop */
1590              } /* end of face 0 */              } /* end of face 0 */
1591              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1592  #pragma omp for nowait  #pragma omp for nowait
1593                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1594                      const double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0)), numComp*sizeof(double));
1595                      const double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0)), numComp*sizeof(double));
1596                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1597                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1598                          o[INDEX2(i,numComp,0)] = f_10[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c1*f_10[i] + c0*f_11[i];
1599                          o[INDEX2(i,numComp,1)] = f_10[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c1*f_11[i] + c0*f_10[i];
1600                      } /* end of component loop i */                      } /* end of component loop i */
1601                  } /* end of k1 loop */                  } /* end of k1 loop */
1602              } /* end of face 1 */              } /* end of face 1 */
1603              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1604  #pragma omp for nowait  #pragma omp for nowait
1605                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1606                      const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                      memcpy(&f_00[0], in.getSampleDataRO(INDEX2(k0,0, m_N0)), numComp*sizeof(double));
1607                      const double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,0, m_N0)), numComp*sizeof(double));
1608                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1609                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1610                          o[INDEX2(i,numComp,0)] = f_00[i]*c1 + f_10[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_10[i] + c1*f_00[i];
1611                          o[INDEX2(i,numComp,1)] = f_00[i]*c0 + f_10[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_00[i] + c1*f_10[i];
1612                      } /* end of component loop i */                      } /* end of component loop i */
1613                  } /* end of k0 loop */                  } /* end of k0 loop */
1614              } /* end of face 2 */              } /* end of face 2 */
1615              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1616  #pragma omp for nowait  #pragma omp for nowait
1617                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1618                      const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0)), numComp*sizeof(double));
1619                      const double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0)), numComp*sizeof(double));
1620                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1621                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1622                          o[INDEX2(i,numComp,0)] = f_01[i]*c1 + f_11[i]*c0;                          o[INDEX2(i,numComp,0)] = c0*f_11[i] + c1*f_01[i];
1623                          o[INDEX2(i,numComp,1)] = f_01[i]*c0 + f_11[i]*c1;                          o[INDEX2(i,numComp,1)] = c0*f_01[i] + c1*f_11[i];
1624                      } /* end of component loop i */                      } /* end of component loop i */
1625                  } /* end of k0 loop */                  } /* end of k0 loop */
1626              } /* end of face 3 */              } /* end of face 3 */
1627          } // end of parallel section          } /* end of parallel section */
1628      }      }
1629  }  }
1630    
# Line 3376  void Rectangle::assemblePDEBoundarySingl Line 3544  void Rectangle::assemblePDEBoundarySingl
3544      {      {
3545          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3546              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3547  #pragma omp for nowait  #pragma omp for
3548                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3549                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3550                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3438  void Rectangle::assemblePDEBoundarySingl Line 3606  void Rectangle::assemblePDEBoundarySingl
3606    
3607          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3608              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3609  #pragma omp for nowait  #pragma omp for
3610                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3611                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3612                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3500  void Rectangle::assemblePDEBoundarySingl Line 3668  void Rectangle::assemblePDEBoundarySingl
3668    
3669          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3670              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3671  #pragma omp for nowait  #pragma omp for
3672                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3673                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3674                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3562  void Rectangle::assemblePDEBoundarySingl Line 3730  void Rectangle::assemblePDEBoundarySingl
3730    
3731          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3732              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3733  #pragma omp for nowait  #pragma omp for
3734                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3735                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
3736                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
# Line 3641  void Rectangle::assemblePDEBoundarySingl Line 3809  void Rectangle::assemblePDEBoundarySingl
3809      {      {
3810          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3811              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3812  #pragma omp for nowait  #pragma omp for
3813                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3814                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3815                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3675  void Rectangle::assemblePDEBoundarySingl Line 3843  void Rectangle::assemblePDEBoundarySingl
3843    
3844          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3845              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3846  #pragma omp for nowait  #pragma omp for
3847                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3848                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3849                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3709  void Rectangle::assemblePDEBoundarySingl Line 3877  void Rectangle::assemblePDEBoundarySingl
3877    
3878          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3879              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3880  #pragma omp for nowait  #pragma omp for
3881                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3882                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3883                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3742  void Rectangle::assemblePDEBoundarySingl Line 3910  void Rectangle::assemblePDEBoundarySingl
3910    
3911          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3912              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3913  #pragma omp for nowait  #pragma omp for
3914                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3915                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3916                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3811  void Rectangle::assemblePDEBoundarySyste Line 3979  void Rectangle::assemblePDEBoundarySyste
3979      {      {
3980          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3981              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3982  #pragma omp for nowait  #pragma omp for
3983                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3984                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3985                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 3886  void Rectangle::assemblePDEBoundarySyste Line 4054  void Rectangle::assemblePDEBoundarySyste
4054    
4055          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4056              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4057  #pragma omp for nowait  #pragma omp for
4058                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4059                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4060                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 3961  void Rectangle::assemblePDEBoundarySyste Line 4129  void Rectangle::assemblePDEBoundarySyste
4129    
4130          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4131              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4132  #pragma omp for nowait  #pragma omp for
4133                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4134                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4135                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4036  void Rectangle::assemblePDEBoundarySyste Line 4204  void Rectangle::assemblePDEBoundarySyste
4204    
4205          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4206              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4207  #pragma omp for nowait  #pragma omp for
4208                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4209                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4210                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4135  void Rectangle::assemblePDEBoundarySyste Line 4303  void Rectangle::assemblePDEBoundarySyste
4303      {      {
4304          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
4305              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4306  #pragma omp for nowait  #pragma omp for
4307                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4308                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4309                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4176  void Rectangle::assemblePDEBoundarySyste Line 4344  void Rectangle::assemblePDEBoundarySyste
4344    
4345          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4346              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4347  #pragma omp for nowait  #pragma omp for
4348                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4349                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4350                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4217  void Rectangle::assemblePDEBoundarySyste Line 4385  void Rectangle::assemblePDEBoundarySyste
4385    
4386          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4387              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4388  #pragma omp for nowait  #pragma omp for
4389                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4390                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4391                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4258  void Rectangle::assemblePDEBoundarySyste Line 4426  void Rectangle::assemblePDEBoundarySyste
4426    
4427          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4428              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4429  #pragma omp for nowait  #pragma omp for
4430                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4431                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4432                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);

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

  ViewVC Help
Powered by ViewVC 1.1.26