/[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 4010 by caltinay, Tue Oct 2 06:57:11 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2012 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  #include <ripley/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=max(1, (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 1027  void Rectangle::dofToNodes(escript::Data Line 1173  void Rectangle::dofToNodes(escript::Data
1173                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1174          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1175      }      }
1176        Paso_Coupler_free(coupler);
1177  }  }
1178    
1179  //private  //private
# Line 1309  void Rectangle::interpolateNodesOnElemen Line 1456  void Rectangle::interpolateNodesOnElemen
1456      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1457      if (reduced) {      if (reduced) {
1458          out.requireWrite();          out.requireWrite();
1459          const double c0 = .25;          const double c0 = 0.25;
1460  #pragma omp parallel for  #pragma omp parallel
1461          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1462              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1463                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_01(numComp);
1464                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_10(numComp);
1465                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));              vector<double> f_11(numComp);
1466                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));  #pragma omp for
1467                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1468                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1469                      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));
1470                  } /* end of component loop i */                      memcpy(&f_01[0], in.getSampleDataRO(INDEX2(k0,k1+1, m_N0)), numComp*sizeof(double));
1471              } /* end of k0 loop */                      memcpy(&f_10[0], in.getSampleDataRO(INDEX2(k0+1,k1, m_N0)), numComp*sizeof(double));
1472          } /* end of k1 loop */                      memcpy(&f_11[0], in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0)), numComp*sizeof(double));
1473                        double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1474                        for (index_t i=0; i < numComp; ++i) {
1475                            o[INDEX2(i,numComp,0)] = c0*(f_00[i] + f_01[i] + f_10[i] + f_11[i]);
1476                        } /* end of component loop i */
1477                    } /* end of k0 loop */
1478                } /* end of k1 loop */
1479            } /* end of parallel section */
1480      } else {      } else {
1481          out.requireWrite();          out.requireWrite();
1482          const double c0 = .16666666666666666667;          const double c0 = 0.16666666666666666667;
1483          const double c1 = .044658198738520451079;          const double c1 = 0.044658198738520451079;
1484          const double c2 = .62200846792814621559;          const double c2 = 0.62200846792814621559;
1485  #pragma omp parallel for  #pragma omp parallel
1486          for (index_t k1=0; k1 < m_NE1; ++k1) {          {
1487              for (index_t k0=0; k0 < m_NE0; ++k0) {              vector<double> f_00(numComp);
1488                  const double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));              vector<double> f_01(numComp);
1489                  const double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));              vector<double> f_10(numComp);
1490                  const double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));              vector<double> f_11(numComp);
1491                  const double* f_00 = in.getSampleDataRO(INDEX2(k0,k1, m_N0));  #pragma omp for
1492                  double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));              for (index_t k1=0; k1 < m_NE1; ++k1) {
1493                  for (index_t i=0; i < numComp; ++i) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1494                      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));
1495                      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));
1496                      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));
1497                      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));
1498                  } /* end of component loop i */                      double* o = out.getSampleDataRW(INDEX2(k0,k1,m_NE0));
1499              } /* end of k0 loop */                      for (index_t i=0; i < numComp; ++i) {
1500          } /* end of k1 loop */                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_10[i]) + c1*f_11[i] + c2*f_00[i];
1501                            o[INDEX2(i,numComp,1)] = c0*(f_00[i] + f_11[i]) + c1*f_01[i] + c2*f_10[i];
1502                            o[INDEX2(i,numComp,2)] = c0*(f_00[i] + f_11[i]) + c1*f_10[i] + c2*f_01[i];
1503                            o[INDEX2(i,numComp,3)] = c0*(f_01[i] + f_10[i]) + c1*f_00[i] + c2*f_11[i];
1504                        } /* end of component loop i */
1505                    } /* end of k0 loop */
1506                } /* end of k1 loop */
1507            } /* end of parallel section */
1508      }      }
1509  }  }
1510    
# Line 1354  void Rectangle::interpolateNodesOnFaces( Line 1515  void Rectangle::interpolateNodesOnFaces(
1515      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1516      if (reduced) {      if (reduced) {
1517          out.requireWrite();          out.requireWrite();
1518          const double c0 = .5;          const double c0 = 0.5;
1519  #pragma omp parallel  #pragma omp parallel
1520          {          {
1521                vector<double> f_00(numComp);
1522                vector<double> f_01(numComp);
1523                vector<double> f_10(numComp);
1524                vector<double> f_11(numComp);
1525              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1526  #pragma omp for nowait  #pragma omp for nowait
1527                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1528                      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));
1529                      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));
1530                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1531                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1532                          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 1536  void Rectangle::interpolateNodesOnFaces(
1536              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1537  #pragma omp for nowait  #pragma omp for nowait
1538                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1539                      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));
1540                      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));
1541                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1542                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1543                          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 1547  void Rectangle::interpolateNodesOnFaces(
1547              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1548  #pragma omp for nowait  #pragma omp for nowait
1549                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1550                      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));
1551                      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));
1552                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1553                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1554                          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 1558  void Rectangle::interpolateNodesOnFaces(
1558              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1559  #pragma omp for nowait  #pragma omp for nowait
1560                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1561                      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));
1562                      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));
1563                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1564                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1565                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);                          o[INDEX2(i,numComp,0)] = c0*(f_01[i] + f_11[i]);
1566                      } /* end of component loop i */                      } /* end of component loop i */
1567                  } /* end of k0 loop */                  } /* end of k0 loop */
1568              } /* end of face 3 */              } /* end of face 3 */
1569          } // end of parallel section          } /* end of parallel section */
1570      } else {      } else {
1571          out.requireWrite();          out.requireWrite();
1572          const double c0 = 0.21132486540518711775;          const double c0 = 0.21132486540518711775;
1573          const double c1 = 0.78867513459481288225;          const double c1 = 0.78867513459481288225;
1574  #pragma omp parallel  #pragma omp parallel
1575          {          {
1576                vector<double> f_00(numComp);
1577                vector<double> f_01(numComp);
1578                vector<double> f_10(numComp);
1579                vector<double> f_11(numComp);
1580              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
1581  #pragma omp for nowait  #pragma omp for nowait
1582                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1583                      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));
1584                      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));
1585                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1586                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1587                          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];
1588                          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];
1589                      } /* end of component loop i */                      } /* end of component loop i */
1590                  } /* end of k1 loop */                  } /* end of k1 loop */
1591              } /* end of face 0 */              } /* end of face 0 */
1592              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
1593  #pragma omp for nowait  #pragma omp for nowait
1594                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
1595                      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));
1596                      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));
1597                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1598                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1599                          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];
1600                          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];
1601                      } /* end of component loop i */                      } /* end of component loop i */
1602                  } /* end of k1 loop */                  } /* end of k1 loop */
1603              } /* end of face 1 */              } /* end of face 1 */
1604              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
1605  #pragma omp for nowait  #pragma omp for nowait
1606                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1607                      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));
1608                      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));
1609                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1610                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1611                          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];
1612                          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];
1613                      } /* end of component loop i */                      } /* end of component loop i */
1614                  } /* end of k0 loop */                  } /* end of k0 loop */
1615              } /* end of face 2 */              } /* end of face 2 */
1616              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
1617  #pragma omp for nowait  #pragma omp for nowait
1618                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
1619                      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));
1620                      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));
1621                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1622                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
1623                          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];
1624                          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];
1625                      } /* end of component loop i */                      } /* end of component loop i */
1626                  } /* end of k0 loop */                  } /* end of k0 loop */
1627              } /* end of face 3 */              } /* end of face 3 */
1628          } // end of parallel section          } /* end of parallel section */
1629      }      }
1630  }  }
1631    
# Line 3376  void Rectangle::assemblePDEBoundarySingl Line 3545  void Rectangle::assemblePDEBoundarySingl
3545      {      {
3546          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3547              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3548  #pragma omp for nowait  #pragma omp for
3549                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3550                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3551                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3438  void Rectangle::assemblePDEBoundarySingl Line 3607  void Rectangle::assemblePDEBoundarySingl
3607    
3608          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3609              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3610  #pragma omp for nowait  #pragma omp for
3611                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3612                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3613                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3500  void Rectangle::assemblePDEBoundarySingl Line 3669  void Rectangle::assemblePDEBoundarySingl
3669    
3670          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3671              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3672  #pragma omp for nowait  #pragma omp for
3673                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3674                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3675                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3562  void Rectangle::assemblePDEBoundarySingl Line 3731  void Rectangle::assemblePDEBoundarySingl
3731    
3732          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3733              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3734  #pragma omp for nowait  #pragma omp for
3735                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3736                      const index_t e = m_faceOffset[3]+k0;                      const index_t e = m_faceOffset[3]+k0;
3737                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
# Line 3641  void Rectangle::assemblePDEBoundarySingl Line 3810  void Rectangle::assemblePDEBoundarySingl
3810      {      {
3811          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3812              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3813  #pragma omp for nowait  #pragma omp for
3814                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3815                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3816                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3675  void Rectangle::assemblePDEBoundarySingl Line 3844  void Rectangle::assemblePDEBoundarySingl
3844    
3845          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
3846              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
3847  #pragma omp for nowait  #pragma omp for
3848                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3849                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3850                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3709  void Rectangle::assemblePDEBoundarySingl Line 3878  void Rectangle::assemblePDEBoundarySingl
3878    
3879          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
3880              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3881  #pragma omp for nowait  #pragma omp for
3882                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3883                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3884                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3742  void Rectangle::assemblePDEBoundarySingl Line 3911  void Rectangle::assemblePDEBoundarySingl
3911    
3912          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
3913              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
3914  #pragma omp for nowait  #pragma omp for
3915                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
3916                      vector<double> EM_S(4*4, 0);                      vector<double> EM_S(4*4, 0);
3917                      vector<double> EM_F(4, 0);                      vector<double> EM_F(4, 0);
# Line 3811  void Rectangle::assemblePDEBoundarySyste Line 3980  void Rectangle::assemblePDEBoundarySyste
3980      {      {
3981          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
3982              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
3983  #pragma omp for nowait  #pragma omp for
3984                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
3985                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
3986                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 3886  void Rectangle::assemblePDEBoundarySyste Line 4055  void Rectangle::assemblePDEBoundarySyste
4055    
4056          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4057              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4058  #pragma omp for nowait  #pragma omp for
4059                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4060                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4061                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 3961  void Rectangle::assemblePDEBoundarySyste Line 4130  void Rectangle::assemblePDEBoundarySyste
4130    
4131          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4132              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4133  #pragma omp for nowait  #pragma omp for
4134                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4135                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4136                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4036  void Rectangle::assemblePDEBoundarySyste Line 4205  void Rectangle::assemblePDEBoundarySyste
4205    
4206          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4207              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4208  #pragma omp for nowait  #pragma omp for
4209                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4210                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4211                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4135  void Rectangle::assemblePDEBoundarySyste Line 4304  void Rectangle::assemblePDEBoundarySyste
4304      {      {
4305          if (m_faceOffset[0] > -1) {          if (m_faceOffset[0] > -1) {
4306              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring
4307  #pragma omp for nowait  #pragma omp for
4308                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4309                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4310                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4176  void Rectangle::assemblePDEBoundarySyste Line 4345  void Rectangle::assemblePDEBoundarySyste
4345    
4346          if (m_faceOffset[1] > -1) {          if (m_faceOffset[1] > -1) {
4347              for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring                          for (index_t k1_0=0; k1_0<2; k1_0++) { // colouring            
4348  #pragma omp for nowait  #pragma omp for
4349                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {                  for (index_t k1=k1_0; k1<m_NE1; k1+=2) {
4350                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4351                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4217  void Rectangle::assemblePDEBoundarySyste Line 4386  void Rectangle::assemblePDEBoundarySyste
4386    
4387          if (m_faceOffset[2] > -1) {          if (m_faceOffset[2] > -1) {
4388              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4389  #pragma omp for nowait  #pragma omp for
4390                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4391                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4392                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);
# Line 4258  void Rectangle::assemblePDEBoundarySyste Line 4427  void Rectangle::assemblePDEBoundarySyste
4427    
4428          if (m_faceOffset[3] > -1) {          if (m_faceOffset[3] > -1) {
4429              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring              for (index_t k0_0=0; k0_0<2; k0_0++) { // colouring
4430  #pragma omp for nowait  #pragma omp for
4431                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {                  for (index_t k0 = k0_0; k0 < m_NE0; k0+=2) {
4432                      vector<double> EM_S(4*4*numEq*numComp, 0);                      vector<double> EM_S(4*4*numEq*numComp, 0);
4433                      vector<double> EM_F(4*numEq, 0);                      vector<double> EM_F(4*numEq, 0);

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

  ViewVC Help
Powered by ViewVC 1.1.26