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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3744 - (hide annotations)
Tue Dec 13 06:41:54 2011 UTC (7 years, 3 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 94091 byte(s)
Implemented (Face)Elements->Reduced(Face)Elements interpolation and started
separating DOF and nodes.
Also, implemented operator==() so that a==b for different domain objects a and
b which are in the same state.

1 caltinay 3691
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2011 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     #include <ripley/Brick.h>
15     extern "C" {
16     #include "paso/SystemMatrixPattern.h"
17     }
18    
19     #if USE_SILO
20     #include <silo.h>
21     #ifdef ESYS_MPI
22     #include <pmpio.h>
23     #endif
24     #endif
25    
26     #include <iomanip>
27    
28     using namespace std;
29    
30     namespace ripley {
31    
32     Brick::Brick(int n0, int n1, int n2, double l0, double l1, double l2, int d0,
33     int d1, int d2) :
34     RipleyDomain(3),
35     m_gNE0(n0),
36     m_gNE1(n1),
37     m_gNE2(n2),
38     m_l0(l0),
39     m_l1(l1),
40     m_l2(l2),
41     m_NX(d0),
42     m_NY(d1),
43     m_NZ(d2)
44     {
45     // ensure number of subdivisions is valid and nodes can be distributed
46     // among number of ranks
47     if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
48     throw RipleyException("Invalid number of spatial subdivisions");
49    
50     if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)
51     throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
52    
53     // local number of elements
54     m_NE0 = n0/m_NX;
55     m_NE1 = n1/m_NY;
56     m_NE2 = n2/m_NZ;
57     // local number of nodes (not necessarily owned)
58     m_N0 = m_NE0+1;
59     m_N1 = m_NE1+1;
60     m_N2 = m_NE2+1;
61     // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
62     m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
63     m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);
64     m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));
65     populateSampleIds();
66     }
67    
68    
69     Brick::~Brick()
70     {
71     }
72    
73     string Brick::getDescription() const
74     {
75     return "ripley::Brick";
76     }
77    
78     bool Brick::operator==(const AbstractDomain& other) const
79     {
80 caltinay 3744 const Brick* o=dynamic_cast<const Brick*>(&other);
81     if (o) {
82     return (RipleyDomain::operator==(other) &&
83     m_gNE0==o->m_gNE0 && m_gNE1==o->m_gNE1 && m_gNE2==o->m_gNE2
84     && m_l0==o->m_l0 && m_l1==o->m_l1 && m_l2==o->m_l2
85     && m_NX==o->m_NX && m_NY==o->m_NY && m_NZ==o->m_NZ);
86     }
87 caltinay 3691
88     return false;
89     }
90    
91     void Brick::dump(const string& fileName) const
92     {
93     #if USE_SILO
94     string fn(fileName);
95     if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
96     fn+=".silo";
97     }
98    
99     const int NUM_SILO_FILES = 1;
100     const char* blockDirFmt = "/block%04d";
101     int driver=DB_HDF5;
102     string siloPath;
103     DBfile* dbfile = NULL;
104    
105     #ifdef ESYS_MPI
106     PMPIO_baton_t* baton = NULL;
107     #endif
108    
109     if (m_mpiInfo->size > 1) {
110     #ifdef ESYS_MPI
111     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
112     0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
113     PMPIO_DefaultClose, (void*)&driver);
114     // try the fallback driver in case of error
115     if (!baton && driver != DB_PDB) {
116     driver = DB_PDB;
117     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
118     0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
119     PMPIO_DefaultClose, (void*)&driver);
120     }
121     if (baton) {
122     char str[64];
123     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
124     siloPath = str;
125     dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
126     }
127     #endif
128     } else {
129     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
130     getDescription().c_str(), driver);
131     // try the fallback driver in case of error
132     if (!dbfile && driver != DB_PDB) {
133     driver = DB_PDB;
134     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
135     getDescription().c_str(), driver);
136     }
137     }
138    
139     if (!dbfile)
140     throw RipleyException("dump: Could not create Silo file");
141    
142     /*
143     if (driver==DB_HDF5) {
144     // gzip level 1 already provides good compression with minimal
145     // performance penalty. Some tests showed that gzip levels >3 performed
146     // rather badly on escript data both in terms of time and space
147     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
148     }
149     */
150    
151     boost::scoped_ptr<double> x(new double[m_N0]);
152     boost::scoped_ptr<double> y(new double[m_N1]);
153     boost::scoped_ptr<double> z(new double[m_N2]);
154     double* coords[3] = { x.get(), y.get(), z.get() };
155 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
156     pair<double,double> ydy = getFirstCoordAndSpacing(1);
157     pair<double,double> zdz = getFirstCoordAndSpacing(2);
158 caltinay 3691 #pragma omp parallel
159     {
160     #pragma omp for
161     for (dim_t i0 = 0; i0 < m_N0; i0++) {
162 caltinay 3698 coords[0][i0]=xdx.first+i0*xdx.second;
163 caltinay 3691 }
164     #pragma omp for
165     for (dim_t i1 = 0; i1 < m_N1; i1++) {
166 caltinay 3698 coords[1][i1]=ydy.first+i1*ydy.second;
167 caltinay 3691 }
168     #pragma omp for
169     for (dim_t i2 = 0; i2 < m_N2; i2++) {
170 caltinay 3698 coords[2][i2]=zdz.first+i2*zdz.second;
171 caltinay 3691 }
172     }
173 caltinay 3698 IndexVector dims = getNumNodesPerDim();
174     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,
175 caltinay 3691 DB_COLLINEAR, NULL);
176    
177 caltinay 3698 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
178     NULL, 0, DB_INT, DB_NODECENT, NULL);
179 caltinay 3691
180 caltinay 3698 // write element ids
181     dims = getNumElementsPerDim();
182     DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
183     &dims[0], 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
184    
185     // rank 0 writes multimesh and multivar
186 caltinay 3691 if (m_mpiInfo->rank == 0) {
187     vector<string> tempstrings;
188     vector<char*> names;
189     for (dim_t i=0; i<m_mpiInfo->size; i++) {
190     stringstream path;
191     path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
192     tempstrings.push_back(path.str());
193     names.push_back((char*)tempstrings.back().c_str());
194     }
195     vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
196     DBSetDir(dbfile, "/");
197     DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
198     &types[0], NULL);
199     tempstrings.clear();
200     names.clear();
201     for (dim_t i=0; i<m_mpiInfo->size; i++) {
202     stringstream path;
203     path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
204     tempstrings.push_back(path.str());
205     names.push_back((char*)tempstrings.back().c_str());
206     }
207     types.assign(m_mpiInfo->size, DB_QUADVAR);
208     DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
209     &types[0], NULL);
210 caltinay 3698 tempstrings.clear();
211     names.clear();
212     for (dim_t i=0; i<m_mpiInfo->size; i++) {
213     stringstream path;
214     path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
215     tempstrings.push_back(path.str());
216     names.push_back((char*)tempstrings.back().c_str());
217     }
218     DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
219     &types[0], NULL);
220 caltinay 3691 }
221    
222     if (m_mpiInfo->size > 1) {
223     #ifdef ESYS_MPI
224     PMPIO_HandOffBaton(baton, dbfile);
225     PMPIO_Finish(baton);
226     #endif
227     } else {
228     DBClose(dbfile);
229     }
230    
231     #else // USE_SILO
232     throw RipleyException("dump(): no Silo support");
233     #endif
234     }
235    
236     const int* Brick::borrowSampleReferenceIDs(int fsType) const
237     {
238 caltinay 3697 switch (fsType) {
239     case Nodes:
240 caltinay 3733 case DegreesOfFreedom: //FIXME
241     case ReducedDegreesOfFreedom: //FIXME
242 caltinay 3697 return &m_nodeId[0];
243     case Elements:
244 caltinay 3733 case ReducedElements:
245 caltinay 3697 return &m_elementId[0];
246 caltinay 3733 case ReducedFaceElements:
247 caltinay 3697 return &m_faceId[0];
248     default:
249     break;
250     }
251 caltinay 3691
252 caltinay 3697 stringstream msg;
253     msg << "borrowSampleReferenceIDs() not implemented for function space type "
254     << fsType;
255     throw RipleyException(msg.str());
256 caltinay 3691 }
257    
258     bool Brick::ownSample(int fsCode, index_t id) const
259     {
260     #ifdef ESYS_MPI
261     if (fsCode == Nodes) {
262 caltinay 3698 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
263     const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
264 caltinay 3691 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
265     } else
266     throw RipleyException("ownSample() only implemented for Nodes");
267     #else
268     return true;
269     #endif
270     }
271    
272 caltinay 3703 void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const
273     {
274     escript::Data& in = *const_cast<escript::Data*>(&cIn);
275     const dim_t numComp = in.getDataPointSize();
276     const double h0 = m_l0/m_gNE0;
277     const double h1 = m_l1/m_gNE1;
278     const double h2 = m_l1/m_gNE2;
279 caltinay 3731 const double C0 = .044658198738520451079;
280     const double C1 = .16666666666666666667;
281     const double C2 = .21132486540518711775;
282     const double C3 = .25;
283     const double C4 = .5;
284     const double C5 = .62200846792814621559;
285     const double C6 = .78867513459481288225;
286    
287 caltinay 3703 if (out.getFunctionSpace().getTypeCode() == Elements) {
288 caltinay 3731 /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
289 caltinay 3703 #pragma omp parallel for
290 caltinay 3707 for (index_t k2=0; k2 < m_NE2; ++k2) {
291     for (index_t k1=0; k1 < m_NE1; ++k1) {
292     for (index_t k0=0; k0 < m_NE0; ++k0) {
293 caltinay 3703 const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
294     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
295     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
296     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
297     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
298     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
299     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
300     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
301     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
302     for (index_t i=0; i < numComp; ++i) {
303 caltinay 3731 const double V0=((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
304     const double V1=((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
305     const double V2=((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
306     const double V3=((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
307     const double V4=((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
308     const double V5=((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
309     const double V6=((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
310     const double V7=((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
311     const double V8=((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
312     const double V9=((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
313     const double V10=((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
314     const double V11=((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
315     o[INDEX3(i,0,0,numComp,3)] = V0;
316     o[INDEX3(i,1,0,numComp,3)] = V4;
317     o[INDEX3(i,2,0,numComp,3)] = V8;
318     o[INDEX3(i,0,1,numComp,3)] = V0;
319     o[INDEX3(i,1,1,numComp,3)] = V5;
320     o[INDEX3(i,2,1,numComp,3)] = V9;
321     o[INDEX3(i,0,2,numComp,3)] = V1;
322     o[INDEX3(i,1,2,numComp,3)] = V4;
323     o[INDEX3(i,2,2,numComp,3)] = V10;
324     o[INDEX3(i,0,3,numComp,3)] = V1;
325     o[INDEX3(i,1,3,numComp,3)] = V5;
326     o[INDEX3(i,2,3,numComp,3)] = V11;
327     o[INDEX3(i,0,4,numComp,3)] = V2;
328     o[INDEX3(i,1,4,numComp,3)] = V6;
329     o[INDEX3(i,2,4,numComp,3)] = V8;
330     o[INDEX3(i,0,5,numComp,3)] = V2;
331     o[INDEX3(i,1,5,numComp,3)] = V7;
332     o[INDEX3(i,2,5,numComp,3)] = V9;
333     o[INDEX3(i,0,6,numComp,3)] = V3;
334     o[INDEX3(i,1,6,numComp,3)] = V6;
335     o[INDEX3(i,2,6,numComp,3)] = V10;
336     o[INDEX3(i,0,7,numComp,3)] = V3;
337     o[INDEX3(i,1,7,numComp,3)] = V7;
338     o[INDEX3(i,2,7,numComp,3)] = V11;
339 caltinay 3703 } /* end of component loop i */
340     } /* end of k0 loop */
341     } /* end of k1 loop */
342     } /* end of k2 loop */
343     /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
344 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
345 caltinay 3731 /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
346 caltinay 3711 #pragma omp parallel for
347     for (index_t k2=0; k2 < m_NE2; ++k2) {
348     for (index_t k1=0; k1 < m_NE1; ++k1) {
349     for (index_t k0=0; k0 < m_NE0; ++k0) {
350     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
351     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
352     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
353     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
354     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
355     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
356     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
357     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
358     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
359     for (index_t i=0; i < numComp; ++i) {
360 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
361     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
362     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
363 caltinay 3711 } /* end of component loop i */
364     } /* end of k0 loop */
365     } /* end of k1 loop */
366     } /* end of k2 loop */
367     /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
368 caltinay 3707 } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
369 caltinay 3731 /*** GENERATOR SNIP_GRAD_FACES TOP */
370 caltinay 3722 #pragma omp parallel
371     {
372     if (m_faceOffset[0] > -1) {
373     #pragma omp for nowait
374     for (index_t k2=0; k2 < m_NE2; ++k2) {
375     for (index_t k1=0; k1 < m_NE1; ++k1) {
376     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
377     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
378     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
379     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
380     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
381     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
382     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
383     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
384     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
385     for (index_t i=0; i < numComp; ++i) {
386 caltinay 3731 const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;
387     const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / h1;
388     const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / h2;
389     const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / h2;
390     o[INDEX3(i,0,0,numComp,3)] = ((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
391     o[INDEX3(i,1,0,numComp,3)] = V0;
392     o[INDEX3(i,2,0,numComp,3)] = V2;
393     o[INDEX3(i,0,1,numComp,3)] = ((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
394     o[INDEX3(i,1,1,numComp,3)] = V0;
395     o[INDEX3(i,2,1,numComp,3)] = V3;
396     o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
397     o[INDEX3(i,1,2,numComp,3)] = V1;
398     o[INDEX3(i,2,2,numComp,3)] = V2;
399     o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
400     o[INDEX3(i,1,3,numComp,3)] = V1;
401     o[INDEX3(i,2,3,numComp,3)] = V3;
402 caltinay 3722 } /* end of component loop i */
403     } /* end of k1 loop */
404     } /* end of k2 loop */
405     } /* end of face 0 */
406     if (m_faceOffset[1] > -1) {
407     #pragma omp for nowait
408     for (index_t k2=0; k2 < m_NE2; ++k2) {
409     for (index_t k1=0; k1 < m_NE1; ++k1) {
410     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
411     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
412     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
413     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
414     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
415     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
416     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
417     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
418     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
419     for (index_t i=0; i < numComp; ++i) {
420 caltinay 3731 const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
421     const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
422     const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
423     const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
424     o[INDEX3(i,0,0,numComp,3)] = ((f_100[i]-f_000[i])*C5 + (f_111[i]-f_011[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
425     o[INDEX3(i,1,0,numComp,3)] = V0;
426     o[INDEX3(i,2,0,numComp,3)] = V2;
427     o[INDEX3(i,0,1,numComp,3)] = ((f_110[i]-f_010[i])*C5 + (f_101[i]-f_001[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
428     o[INDEX3(i,1,1,numComp,3)] = V0;
429     o[INDEX3(i,2,1,numComp,3)] = V3;
430     o[INDEX3(i,0,2,numComp,3)] = ((f_101[i]-f_001[i])*C5 + (f_110[i]-f_010[i])*C0 + (f_100[i]+f_111[i]-f_000[i]-f_011[i])*C1) / h0;
431     o[INDEX3(i,1,2,numComp,3)] = V1;
432     o[INDEX3(i,2,2,numComp,3)] = V2;
433     o[INDEX3(i,0,3,numComp,3)] = ((f_111[i]-f_011[i])*C5 + (f_100[i]-f_000[i])*C0 + (f_101[i]+f_110[i]-f_001[i]-f_010[i])*C1) / h0;
434     o[INDEX3(i,1,3,numComp,3)] = V1;
435     o[INDEX3(i,2,3,numComp,3)] = V3;
436 caltinay 3722 } /* end of component loop i */
437     } /* end of k1 loop */
438     } /* end of k2 loop */
439     } /* end of face 1 */
440     if (m_faceOffset[2] > -1) {
441     #pragma omp for nowait
442     for (index_t k2=0; k2 < m_NE2; ++k2) {
443     for (index_t k0=0; k0 < m_NE0; ++k0) {
444     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
445     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
446     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
447     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
448     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
449     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
450     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
451     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
452     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
453     for (index_t i=0; i < numComp; ++i) {
454 caltinay 3731 const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;
455     const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / h2;
456     const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / h2;
457     o[INDEX3(i,0,0,numComp,3)] = V0;
458     o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
459     o[INDEX3(i,2,0,numComp,3)] = V1;
460     o[INDEX3(i,0,1,numComp,3)] = V0;
461     o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
462     o[INDEX3(i,2,1,numComp,3)] = V2;
463     o[INDEX3(i,0,2,numComp,3)] = V0;
464     o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
465     o[INDEX3(i,2,2,numComp,3)] = V1;
466     o[INDEX3(i,0,3,numComp,3)] = V0;
467     o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
468     o[INDEX3(i,2,3,numComp,3)] = V2;
469 caltinay 3722 } /* end of component loop i */
470     } /* end of k0 loop */
471     } /* end of k2 loop */
472     } /* end of face 2 */
473     if (m_faceOffset[3] > -1) {
474     #pragma omp for nowait
475     for (index_t k2=0; k2 < m_NE2; ++k2) {
476     for (index_t k0=0; k0 < m_NE0; ++k0) {
477     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
478     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
479     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
480     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
481     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
482     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
483     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
484     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
485     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
486     for (index_t i=0; i < numComp; ++i) {
487 caltinay 3731 const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
488     const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
489     const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
490     const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
491     o[INDEX3(i,0,0,numComp,3)] = V0;
492     o[INDEX3(i,1,0,numComp,3)] = ((f_010[i]-f_000[i])*C5 + (f_111[i]-f_101[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
493     o[INDEX3(i,2,0,numComp,3)] = V2;
494     o[INDEX3(i,0,1,numComp,3)] = V0;
495     o[INDEX3(i,1,1,numComp,3)] = ((f_110[i]-f_100[i])*C5 + (f_011[i]-f_001[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
496     o[INDEX3(i,2,1,numComp,3)] = V3;
497     o[INDEX3(i,0,2,numComp,3)] = V1;
498     o[INDEX3(i,1,2,numComp,3)] = ((f_011[i]-f_001[i])*C5 + (f_110[i]-f_100[i])*C0 + (f_010[i]+f_111[i]-f_000[i]-f_101[i])*C1) / h1;
499     o[INDEX3(i,2,2,numComp,3)] = V2;
500     o[INDEX3(i,0,3,numComp,3)] = V1;
501     o[INDEX3(i,1,3,numComp,3)] = ((f_111[i]-f_101[i])*C5 + (f_010[i]-f_000[i])*C0 + (f_011[i]+f_110[i]-f_001[i]-f_100[i])*C1) / h1;
502     o[INDEX3(i,2,3,numComp,3)] = V3;
503 caltinay 3722 } /* end of component loop i */
504     } /* end of k0 loop */
505     } /* end of k2 loop */
506     } /* end of face 3 */
507     if (m_faceOffset[4] > -1) {
508     #pragma omp for nowait
509 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
510 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
511     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
512     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
513     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
514     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
515     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
516     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
517     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
518     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
519     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
520     for (index_t i=0; i < numComp; ++i) {
521 caltinay 3731 const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;
522     const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / h0;
523     const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / h1;
524     const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / h1;
525     o[INDEX3(i,0,0,numComp,3)] = V0;
526     o[INDEX3(i,1,0,numComp,3)] = V2;
527     o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
528     o[INDEX3(i,0,1,numComp,3)] = V0;
529     o[INDEX3(i,1,1,numComp,3)] = V3;
530     o[INDEX3(i,2,1,numComp,3)] = ((f_101[i]-f_100[i])*C5 + (f_011[i]-f_010[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
531     o[INDEX3(i,0,2,numComp,3)] = V1;
532     o[INDEX3(i,1,2,numComp,3)] = V2;
533     o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
534     o[INDEX3(i,0,3,numComp,3)] = V1;
535     o[INDEX3(i,1,3,numComp,3)] = V3;
536     o[INDEX3(i,2,3,numComp,3)] = ((f_111[i]-f_110[i])*C5 + (f_001[i]-f_000[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
537 caltinay 3722 } /* end of component loop i */
538     } /* end of k0 loop */
539 caltinay 3707 } /* end of k1 loop */
540 caltinay 3722 } /* end of face 4 */
541     if (m_faceOffset[5] > -1) {
542     #pragma omp for nowait
543 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
544 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
545     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
546     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
547     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
548     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
549     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
550     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
551     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
552     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
553     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
554     for (index_t i=0; i < numComp; ++i) {
555 caltinay 3731 const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
556     const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
557     const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
558     const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
559     o[INDEX3(i,0,0,numComp,3)] = V0;
560     o[INDEX3(i,1,0,numComp,3)] = V2;
561     o[INDEX3(i,2,0,numComp,3)] = ((f_001[i]-f_000[i])*C5 + (f_111[i]-f_110[i])*C0 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
562     o[INDEX3(i,0,1,numComp,3)] = V0;
563     o[INDEX3(i,1,1,numComp,3)] = V3;
564     o[INDEX3(i,2,1,numComp,3)] = ((f_011[i]-f_010[i])*C0 + (f_101[i]-f_100[i])*C5 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
565     o[INDEX3(i,0,2,numComp,3)] = V1;
566     o[INDEX3(i,1,2,numComp,3)] = V2;
567     o[INDEX3(i,2,2,numComp,3)] = ((f_011[i]-f_010[i])*C5 + (f_101[i]-f_100[i])*C0 + (f_001[i]+f_111[i]-f_000[i]-f_110[i])*C1) / h2;
568     o[INDEX3(i,0,3,numComp,3)] = V1;
569     o[INDEX3(i,1,3,numComp,3)] = V3;
570     o[INDEX3(i,2,3,numComp,3)] = ((f_001[i]-f_000[i])*C0 + (f_111[i]-f_110[i])*C5 + (f_011[i]+f_101[i]-f_010[i]-f_100[i])*C1) / h2;
571 caltinay 3722 } /* end of component loop i */
572     } /* end of k0 loop */
573 caltinay 3707 } /* end of k1 loop */
574 caltinay 3722 } /* end of face 5 */
575     } // end of parallel section
576 caltinay 3707 /* GENERATOR SNIP_GRAD_FACES BOTTOM */
577 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
578 caltinay 3731 /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
579 caltinay 3722 #pragma omp parallel
580     {
581     if (m_faceOffset[0] > -1) {
582     #pragma omp for nowait
583     for (index_t k2=0; k2 < m_NE2; ++k2) {
584     for (index_t k1=0; k1 < m_NE1; ++k1) {
585     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
586     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
587     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
588     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
589     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
590     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
591     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
592     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
593     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
594     for (index_t i=0; i < numComp; ++i) {
595 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
596     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;
597     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;
598 caltinay 3722 } /* end of component loop i */
599     } /* end of k1 loop */
600     } /* end of k2 loop */
601     } /* end of face 0 */
602     if (m_faceOffset[1] > -1) {
603     #pragma omp for nowait
604     for (index_t k2=0; k2 < m_NE2; ++k2) {
605     for (index_t k1=0; k1 < m_NE1; ++k1) {
606     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
607     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
608     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
609     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
610     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
611     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
612     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
613     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
614     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
615     for (index_t i=0; i < numComp; ++i) {
616 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_010[i]-f_011[i])*C3 / h0;
617     o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;
618     o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;
619 caltinay 3722 } /* end of component loop i */
620     } /* end of k1 loop */
621     } /* end of k2 loop */
622     } /* end of face 1 */
623     if (m_faceOffset[2] > -1) {
624     #pragma omp for nowait
625     for (index_t k2=0; k2 < m_NE2; ++k2) {
626     for (index_t k0=0; k0 < m_NE0; ++k0) {
627     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
628     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
629     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
630     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
631     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
632     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
633     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
634     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
635     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
636     for (index_t i=0; i < numComp; ++i) {
637 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
638     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
639     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;
640 caltinay 3722 } /* end of component loop i */
641     } /* end of k0 loop */
642     } /* end of k2 loop */
643     } /* end of face 2 */
644     if (m_faceOffset[3] > -1) {
645     #pragma omp for nowait
646     for (index_t k2=0; k2 < m_NE2; ++k2) {
647     for (index_t k0=0; k0 < m_NE0; ++k0) {
648     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
649     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
650     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
651     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
652     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
653     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
654     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
655     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
656     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
657     for (index_t i=0; i < numComp; ++i) {
658 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
659     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]+f_110[i]+f_111[i]-f_000[i]-f_001[i]-f_100[i]-f_101[i])*C3 / h1;
660     o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;
661 caltinay 3722 } /* end of component loop i */
662     } /* end of k0 loop */
663     } /* end of k2 loop */
664     } /* end of face 3 */
665     if (m_faceOffset[4] > -1) {
666     #pragma omp for nowait
667 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
668 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
669     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
670     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
671     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
672     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
673     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
674     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
675     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
676     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
677     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
678     for (index_t i=0; i < numComp; ++i) {
679 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
680     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;
681     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C4 / h2;
682 caltinay 3722 } /* end of component loop i */
683     } /* end of k0 loop */
684 caltinay 3711 } /* end of k1 loop */
685 caltinay 3722 } /* end of face 4 */
686     if (m_faceOffset[5] > -1) {
687     #pragma omp for nowait
688 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
689 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
690     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
691     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
692     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
693     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
694     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
695     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
696     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
697     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
698     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
699     for (index_t i=0; i < numComp; ++i) {
700 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
701     o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;
702     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]+f_101[i]+f_111[i]-f_000[i]-f_010[i]-f_100[i]-f_110[i])*C3 / h2;
703 caltinay 3722 } /* end of component loop i */
704     } /* end of k0 loop */
705 caltinay 3711 } /* end of k1 loop */
706 caltinay 3722 } /* end of face 5 */
707     } // end of parallel section
708 caltinay 3711 /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
709 caltinay 3703 } else {
710 caltinay 3707 stringstream msg;
711     msg << "setToGradient() not implemented for "
712     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
713     throw RipleyException(msg.str());
714 caltinay 3703 }
715     }
716    
717 caltinay 3713 void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
718     {
719     escript::Data& in = *const_cast<escript::Data*>(&arg);
720     const dim_t numComp = in.getDataPointSize();
721     const double h0 = m_l0/m_gNE0;
722     const double h1 = m_l1/m_gNE1;
723     const double h2 = m_l2/m_gNE2;
724     if (arg.getFunctionSpace().getTypeCode() == Elements) {
725     const double w_0 = h0*h1*h2/8.;
726     #pragma omp parallel
727     {
728     vector<double> int_local(numComp, 0);
729 caltinay 3722 #pragma omp for nowait
730 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
731     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
732     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
733     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
734     for (index_t i=0; i < numComp; ++i) {
735     const register double f_0 = f[INDEX2(i,0,numComp)];
736     const register double f_1 = f[INDEX2(i,1,numComp)];
737     const register double f_2 = f[INDEX2(i,2,numComp)];
738     const register double f_3 = f[INDEX2(i,3,numComp)];
739     const register double f_4 = f[INDEX2(i,4,numComp)];
740     const register double f_5 = f[INDEX2(i,5,numComp)];
741     const register double f_6 = f[INDEX2(i,6,numComp)];
742     const register double f_7 = f[INDEX2(i,7,numComp)];
743     int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
744     } /* end of component loop i */
745     } /* end of k0 loop */
746     } /* end of k1 loop */
747     } /* end of k2 loop */
748    
749     #pragma omp critical
750     for (index_t i=0; i<numComp; i++)
751     integrals[i]+=int_local[i];
752 caltinay 3722 } // end of parallel section
753 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
754     const double w_0 = h0*h1*h2;
755     #pragma omp parallel
756     {
757     vector<double> int_local(numComp, 0);
758 caltinay 3722 #pragma omp for nowait
759 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
760     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
761     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
762     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
763     for (index_t i=0; i < numComp; ++i) {
764     int_local[i]+=f[i]*w_0;
765     } /* end of component loop i */
766     } /* end of k0 loop */
767     } /* end of k1 loop */
768     } /* end of k2 loop */
769    
770     #pragma omp critical
771     for (index_t i=0; i<numComp; i++)
772     integrals[i]+=int_local[i];
773 caltinay 3722 } // end of parallel section
774 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
775     const double w_0 = h1*h2/4.;
776     const double w_1 = h0*h2/4.;
777     const double w_2 = h0*h1/4.;
778     #pragma omp parallel
779     {
780     vector<double> int_local(numComp, 0);
781     if (m_faceOffset[0] > -1) {
782 caltinay 3722 #pragma omp for nowait
783 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
784     for (index_t k1=0; k1 < m_NE1; ++k1) {
785     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
786     for (index_t i=0; i < numComp; ++i) {
787     const register double f_0 = f[INDEX2(i,0,numComp)];
788     const register double f_1 = f[INDEX2(i,1,numComp)];
789     const register double f_2 = f[INDEX2(i,2,numComp)];
790     const register double f_3 = f[INDEX2(i,3,numComp)];
791     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
792     } /* end of component loop i */
793     } /* end of k1 loop */
794     } /* end of k2 loop */
795     }
796    
797     if (m_faceOffset[1] > -1) {
798 caltinay 3722 #pragma omp for nowait
799 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
800     for (index_t k1=0; k1 < m_NE1; ++k1) {
801     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
802     for (index_t i=0; i < numComp; ++i) {
803     const register double f_0 = f[INDEX2(i,0,numComp)];
804     const register double f_1 = f[INDEX2(i,1,numComp)];
805     const register double f_2 = f[INDEX2(i,2,numComp)];
806     const register double f_3 = f[INDEX2(i,3,numComp)];
807     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
808     } /* end of component loop i */
809     } /* end of k1 loop */
810     } /* end of k2 loop */
811     }
812    
813     if (m_faceOffset[2] > -1) {
814 caltinay 3722 #pragma omp for nowait
815 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
816     for (index_t k0=0; k0 < m_NE0; ++k0) {
817     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
818     for (index_t i=0; i < numComp; ++i) {
819     const register double f_0 = f[INDEX2(i,0,numComp)];
820     const register double f_1 = f[INDEX2(i,1,numComp)];
821     const register double f_2 = f[INDEX2(i,2,numComp)];
822     const register double f_3 = f[INDEX2(i,3,numComp)];
823     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
824     } /* end of component loop i */
825     } /* end of k1 loop */
826     } /* end of k2 loop */
827     }
828    
829     if (m_faceOffset[3] > -1) {
830 caltinay 3722 #pragma omp for nowait
831 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
832     for (index_t k0=0; k0 < m_NE0; ++k0) {
833     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
834     for (index_t i=0; i < numComp; ++i) {
835     const register double f_0 = f[INDEX2(i,0,numComp)];
836     const register double f_1 = f[INDEX2(i,1,numComp)];
837     const register double f_2 = f[INDEX2(i,2,numComp)];
838     const register double f_3 = f[INDEX2(i,3,numComp)];
839     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
840     } /* end of component loop i */
841     } /* end of k1 loop */
842     } /* end of k2 loop */
843     }
844    
845     if (m_faceOffset[4] > -1) {
846 caltinay 3722 #pragma omp for nowait
847 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
848     for (index_t k0=0; k0 < m_NE0; ++k0) {
849     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
850     for (index_t i=0; i < numComp; ++i) {
851     const register double f_0 = f[INDEX2(i,0,numComp)];
852     const register double f_1 = f[INDEX2(i,1,numComp)];
853     const register double f_2 = f[INDEX2(i,2,numComp)];
854     const register double f_3 = f[INDEX2(i,3,numComp)];
855     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
856     } /* end of component loop i */
857     } /* end of k1 loop */
858     } /* end of k2 loop */
859     }
860    
861     if (m_faceOffset[5] > -1) {
862 caltinay 3722 #pragma omp for nowait
863 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
864     for (index_t k0=0; k0 < m_NE0; ++k0) {
865     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
866     for (index_t i=0; i < numComp; ++i) {
867     const register double f_0 = f[INDEX2(i,0,numComp)];
868     const register double f_1 = f[INDEX2(i,1,numComp)];
869     const register double f_2 = f[INDEX2(i,2,numComp)];
870     const register double f_3 = f[INDEX2(i,3,numComp)];
871     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
872     } /* end of component loop i */
873     } /* end of k1 loop */
874     } /* end of k2 loop */
875     }
876    
877     #pragma omp critical
878     for (index_t i=0; i<numComp; i++)
879     integrals[i]+=int_local[i];
880 caltinay 3722 } // end of parallel section
881 caltinay 3713
882     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
883     const double w_0 = h1*h2;
884     const double w_1 = h0*h2;
885     const double w_2 = h0*h1;
886     #pragma omp parallel
887     {
888     vector<double> int_local(numComp, 0);
889     if (m_faceOffset[0] > -1) {
890 caltinay 3722 #pragma omp for nowait
891 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
892     for (index_t k1=0; k1 < m_NE1; ++k1) {
893     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
894     for (index_t i=0; i < numComp; ++i) {
895     int_local[i]+=f[i]*w_0;
896     } /* end of component loop i */
897     } /* end of k1 loop */
898     } /* end of k2 loop */
899     }
900    
901     if (m_faceOffset[1] > -1) {
902 caltinay 3722 #pragma omp for nowait
903 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
904     for (index_t k1=0; k1 < m_NE1; ++k1) {
905     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
906     for (index_t i=0; i < numComp; ++i) {
907     int_local[i]+=f[i]*w_0;
908     } /* end of component loop i */
909     } /* end of k1 loop */
910     } /* end of k2 loop */
911     }
912    
913     if (m_faceOffset[2] > -1) {
914 caltinay 3722 #pragma omp for nowait
915 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
916     for (index_t k0=0; k0 < m_NE0; ++k0) {
917     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
918     for (index_t i=0; i < numComp; ++i) {
919     int_local[i]+=f[i]*w_1;
920     } /* end of component loop i */
921     } /* end of k1 loop */
922     } /* end of k2 loop */
923     }
924    
925     if (m_faceOffset[3] > -1) {
926 caltinay 3722 #pragma omp for nowait
927 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
928     for (index_t k0=0; k0 < m_NE0; ++k0) {
929     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
930     for (index_t i=0; i < numComp; ++i) {
931     int_local[i]+=f[i]*w_1;
932     } /* end of component loop i */
933     } /* end of k1 loop */
934     } /* end of k2 loop */
935     }
936    
937     if (m_faceOffset[4] > -1) {
938 caltinay 3722 #pragma omp for nowait
939 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
940     for (index_t k0=0; k0 < m_NE0; ++k0) {
941     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
942     for (index_t i=0; i < numComp; ++i) {
943     int_local[i]+=f[i]*w_2;
944     } /* end of component loop i */
945     } /* end of k1 loop */
946     } /* end of k2 loop */
947     }
948    
949     if (m_faceOffset[5] > -1) {
950 caltinay 3722 #pragma omp for nowait
951 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
952     for (index_t k0=0; k0 < m_NE0; ++k0) {
953     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
954     for (index_t i=0; i < numComp; ++i) {
955     int_local[i]+=f[i]*w_2;
956     } /* end of component loop i */
957     } /* end of k1 loop */
958     } /* end of k2 loop */
959     }
960    
961     #pragma omp critical
962     for (index_t i=0; i<numComp; i++)
963     integrals[i]+=int_local[i];
964 caltinay 3722 } // end of parallel section
965 caltinay 3713
966     } else {
967     stringstream msg;
968     msg << "setToIntegrals() not implemented for "
969     << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
970     throw RipleyException(msg.str());
971     }
972     }
973    
974 caltinay 3722 void Brick::setToNormal(escript::Data& out) const
975     {
976     if (out.getFunctionSpace().getTypeCode() == FaceElements) {
977     #pragma omp parallel
978     {
979     if (m_faceOffset[0] > -1) {
980     #pragma omp for nowait
981     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
982     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
983     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
984     // set vector at four quadrature points
985     *o++ = -1.; *o++ = 0.; *o++ = 0.;
986     *o++ = -1.; *o++ = 0.; *o++ = 0.;
987     *o++ = -1.; *o++ = 0.; *o++ = 0.;
988     *o++ = -1.; *o++ = 0.; *o = 0.;
989     }
990     }
991     }
992    
993     if (m_faceOffset[1] > -1) {
994     #pragma omp for nowait
995     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
996     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
997     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
998     // set vector at four quadrature points
999     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1000     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1001     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1002     *o++ = 1.; *o++ = 0.; *o = 0.;
1003     }
1004     }
1005     }
1006    
1007     if (m_faceOffset[2] > -1) {
1008     #pragma omp for nowait
1009     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1010     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1011     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1012     // set vector at four quadrature points
1013     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1014     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1015     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1016     *o++ = 0.; *o++ = -1.; *o = 0.;
1017     }
1018     }
1019     }
1020    
1021     if (m_faceOffset[3] > -1) {
1022     #pragma omp for nowait
1023     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1024     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1025     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1026     // set vector at four quadrature points
1027     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1028     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1029     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1030     *o++ = 0.; *o++ = 1.; *o = 0.;
1031     }
1032     }
1033     }
1034    
1035     if (m_faceOffset[4] > -1) {
1036     #pragma omp for nowait
1037     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1038     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1039     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1040     // set vector at four quadrature points
1041     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1042     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1043     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1044     *o++ = 0.; *o++ = 0.; *o = -1.;
1045     }
1046     }
1047     }
1048    
1049     if (m_faceOffset[5] > -1) {
1050     #pragma omp for nowait
1051     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1052     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1053     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1054     // set vector at four quadrature points
1055     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1056     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1057     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1058     *o++ = 0.; *o++ = 0.; *o = 1.;
1059     }
1060     }
1061     }
1062     } // end of parallel section
1063     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1064     #pragma omp parallel
1065     {
1066     if (m_faceOffset[0] > -1) {
1067     #pragma omp for nowait
1068     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1069     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1070     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1071     *o++ = -1.;
1072     *o++ = 0.;
1073     *o = 0.;
1074     }
1075     }
1076     }
1077    
1078     if (m_faceOffset[1] > -1) {
1079     #pragma omp for nowait
1080     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1081     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1082     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1083     *o++ = 1.;
1084     *o++ = 0.;
1085     *o = 0.;
1086     }
1087     }
1088     }
1089    
1090     if (m_faceOffset[2] > -1) {
1091     #pragma omp for nowait
1092     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1093     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1094     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1095     *o++ = 0.;
1096     *o++ = -1.;
1097     *o = 0.;
1098     }
1099     }
1100     }
1101    
1102     if (m_faceOffset[3] > -1) {
1103     #pragma omp for nowait
1104     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1105     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1106     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1107     *o++ = 0.;
1108     *o++ = 1.;
1109     *o = 0.;
1110     }
1111     }
1112     }
1113    
1114     if (m_faceOffset[4] > -1) {
1115     #pragma omp for nowait
1116     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1117     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1118     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1119     *o++ = 0.;
1120     *o++ = 0.;
1121     *o = -1.;
1122     }
1123     }
1124     }
1125    
1126     if (m_faceOffset[5] > -1) {
1127     #pragma omp for nowait
1128     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1129     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1130     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1131     *o++ = 0.;
1132     *o++ = 0.;
1133     *o = 1.;
1134     }
1135     }
1136     }
1137     } // end of parallel section
1138    
1139     } else {
1140     stringstream msg;
1141     msg << "setToNormal() not implemented for "
1142     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
1143     throw RipleyException(msg.str());
1144     }
1145     }
1146    
1147 caltinay 3691 Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
1148     bool reducedColOrder) const
1149     {
1150     if (reducedRowOrder || reducedColOrder)
1151     throw RipleyException("getPattern() not implemented for reduced order");
1152    
1153     throw RipleyException("getPattern() not implemented");
1154     }
1155    
1156     void Brick::Print_Mesh_Info(const bool full) const
1157     {
1158     RipleyDomain::Print_Mesh_Info(full);
1159     if (full) {
1160     cout << " Id Coordinates" << endl;
1161     cout.precision(15);
1162     cout.setf(ios::scientific, ios::floatfield);
1163 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
1164     pair<double,double> ydy = getFirstCoordAndSpacing(1);
1165     pair<double,double> zdz = getFirstCoordAndSpacing(2);
1166 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
1167     cout << " " << setw(5) << m_nodeId[i]
1168 caltinay 3698 << " " << xdx.first+(i%m_N0)*xdx.second
1169     << " " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
1170     << " " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
1171 caltinay 3691 }
1172     }
1173     }
1174    
1175 caltinay 3698 IndexVector Brick::getNumNodesPerDim() const
1176     {
1177     IndexVector ret;
1178     ret.push_back(m_N0);
1179     ret.push_back(m_N1);
1180     ret.push_back(m_N2);
1181     return ret;
1182     }
1183    
1184     IndexVector Brick::getNumElementsPerDim() const
1185     {
1186     IndexVector ret;
1187     ret.push_back(m_NE0);
1188     ret.push_back(m_NE1);
1189     ret.push_back(m_NE2);
1190     return ret;
1191     }
1192    
1193     IndexVector Brick::getNumFacesPerBoundary() const
1194     {
1195     IndexVector ret(6, 0);
1196     //left
1197     if (m_offset0==0)
1198     ret[0]=m_NE1*m_NE2;
1199     //right
1200     if (m_mpiInfo->rank%m_NX==m_NX-1)
1201     ret[1]=m_NE1*m_NE2;
1202     //bottom
1203     if (m_offset1==0)
1204     ret[2]=m_NE0*m_NE2;
1205     //top
1206     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1207     ret[3]=m_NE0*m_NE2;
1208     //front
1209     if (m_offset2==0)
1210     ret[4]=m_NE0*m_NE1;
1211     //back
1212     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1213     ret[5]=m_NE0*m_NE1;
1214     return ret;
1215     }
1216    
1217     pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
1218     {
1219     if (dim==0)
1220     return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
1221     else if (dim==1)
1222     return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
1223     else if (dim==2)
1224     return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
1225    
1226     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
1227     }
1228    
1229 caltinay 3744 //protected
1230     dim_t Brick::getNumDOF() const
1231     {
1232     return m_nodeDistribution[m_mpiInfo->rank+1]
1233     -m_nodeDistribution[m_mpiInfo->rank];
1234     }
1235 caltinay 3698
1236 caltinay 3691 //protected
1237     dim_t Brick::getNumFaceElements() const
1238     {
1239     dim_t n=0;
1240     //left
1241     if (m_offset0==0)
1242     n+=m_NE1*m_NE2;
1243     //right
1244     if (m_mpiInfo->rank%m_NX==m_NX-1)
1245     n+=m_NE1*m_NE2;
1246     //bottom
1247     if (m_offset1==0)
1248     n+=m_NE0*m_NE2;
1249     //top
1250     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1251     n+=m_NE0*m_NE2;
1252     //front
1253     if (m_offset2==0)
1254     n+=m_NE0*m_NE1;
1255     //back
1256     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1257     n+=m_NE0*m_NE1;
1258    
1259     return n;
1260     }
1261    
1262     //protected
1263     void Brick::assembleCoordinates(escript::Data& arg) const
1264     {
1265     escriptDataC x = arg.getDataC();
1266     int numDim = m_numDim;
1267     if (!isDataPointShapeEqual(&x, 1, &numDim))
1268     throw RipleyException("setToX: Invalid Data object shape");
1269     if (!numSamplesEqual(&x, 1, getNumNodes()))
1270     throw RipleyException("setToX: Illegal number of samples in Data object");
1271    
1272 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
1273     pair<double,double> ydy = getFirstCoordAndSpacing(1);
1274     pair<double,double> zdz = getFirstCoordAndSpacing(2);
1275 caltinay 3691 arg.requireWrite();
1276     #pragma omp parallel for
1277     for (dim_t i2 = 0; i2 < m_N2; i2++) {
1278     for (dim_t i1 = 0; i1 < m_N1; i1++) {
1279     for (dim_t i0 = 0; i0 < m_N0; i0++) {
1280     double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
1281 caltinay 3698 point[0] = xdx.first+i0*xdx.second;
1282     point[1] = ydy.first+i1*ydy.second;
1283     point[2] = zdz.first+i2*zdz.second;
1284 caltinay 3691 }
1285     }
1286     }
1287     }
1288    
1289     //private
1290     void Brick::populateSampleIds()
1291     {
1292 caltinay 3697 // identifiers are ordered from left to right, bottom to top, front to back
1293     // on each rank, except for the shared nodes which are owned by the rank
1294     // below / to the left / to the front of the current rank
1295 caltinay 3698
1296     // build node distribution vector first.
1297     // m_nodeDistribution[i] is the first node id on rank i, that is
1298     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1299     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1300     m_nodeDistribution[1]=getNumNodes();
1301     for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
1302     const index_t x = k%m_NX;
1303     const index_t y = k%(m_NX*m_NY)/m_NX;
1304     const index_t z = k/(m_NX*m_NY);
1305     index_t numNodes=getNumNodes();
1306     if (x>0)
1307     numNodes-=m_N1*m_N2;
1308     if (y>0)
1309     numNodes-=m_N0*m_N2;
1310     if (z>0)
1311     numNodes-=m_N0*m_N1;
1312     // if an edge was subtracted twice add it back
1313     if (x>0 && y>0)
1314     numNodes+=m_N2;
1315     if (x>0 && z>0)
1316     numNodes+=m_N1;
1317     if (y>0 && z>0)
1318     numNodes+=m_N0;
1319     // the corner node was removed 3x and added back 3x, so subtract it
1320     if (x>0 && y>0 && z>0)
1321     numNodes--;
1322     m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
1323     }
1324     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
1325    
1326 caltinay 3691 m_nodeId.resize(getNumNodes());
1327 caltinay 3698
1328     // the bottom, left and front planes are not owned by this rank so the
1329     // identifiers need to be computed accordingly
1330     const index_t left = (m_offset0==0 ? 0 : 1);
1331     const index_t bottom = (m_offset1==0 ? 0 : 1);
1332     const index_t front = (m_offset2==0 ? 0 : 1);
1333    
1334     // case 1: all nodes on left plane are owned by rank on the left
1335     if (left>0) {
1336     const int neighbour=m_mpiInfo->rank-1;
1337     const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1338     const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1339 caltinay 3691 #pragma omp parallel for
1340 caltinay 3698 for (dim_t i2=front; i2<m_N2; i2++) {
1341     for (dim_t i1=bottom; i1<m_N1; i1++) {
1342     m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1343     + (i1-bottom+1)*leftN0
1344     + (i2-front)*leftN0*leftN1 - 1;
1345     }
1346     }
1347 caltinay 3691 }
1348 caltinay 3698 // case 2: all nodes on bottom plane are owned by rank below
1349     if (bottom>0) {
1350     const int neighbour=m_mpiInfo->rank-m_NX;
1351     const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1352     const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1353     #pragma omp parallel for
1354     for (dim_t i2=front; i2<m_N2; i2++) {
1355     for (dim_t i0=left; i0<m_N0; i0++) {
1356     m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1357     + bottomN0*(bottomN1-1)
1358     + (i2-front)*bottomN0*bottomN1 + i0-left;
1359     }
1360     }
1361     }
1362     // case 3: all nodes on front plane are owned by rank in front
1363     if (front>0) {
1364     const int neighbour=m_mpiInfo->rank-m_NX*m_NY;
1365     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1366     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1367     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1368     #pragma omp parallel for
1369     for (dim_t i1=bottom; i1<m_N1; i1++) {
1370     for (dim_t i0=left; i0<m_N0; i0++) {
1371     m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]
1372     + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;
1373     }
1374     }
1375     }
1376     // case 4: nodes on front bottom edge are owned by the corresponding rank
1377     if (front>0 && bottom>0) {
1378     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);
1379     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1380     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1381     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1382     #pragma omp parallel for
1383     for (dim_t i0=left; i0<m_N0; i0++) {
1384     m_nodeId[i0]=m_nodeDistribution[neighbour]
1385     + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;
1386     }
1387     }
1388     // case 5: nodes on left bottom edge are owned by the corresponding rank
1389     if (left>0 && bottom>0) {
1390     const int neighbour=m_mpiInfo->rank-m_NX-1;
1391     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1392     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1393     #pragma omp parallel for
1394     for (dim_t i2=front; i2<m_N2; i2++) {
1395     m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1396     + (1+i2-front)*N0*N1-1;
1397     }
1398     }
1399     // case 6: nodes on left front edge are owned by the corresponding rank
1400     if (left>0 && front>0) {
1401     const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;
1402     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1403     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1404     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1405     #pragma omp parallel for
1406     for (dim_t i1=bottom; i1<m_N1; i1++) {
1407     m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
1408     + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;
1409     }
1410     }
1411     // case 7: bottom-left-front corner node owned by corresponding rank
1412     if (left>0 && bottom>0 && front>0) {
1413     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;
1414     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1415     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1416     const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);
1417     m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;
1418     }
1419 caltinay 3697
1420 caltinay 3698 // the rest of the id's are contiguous
1421     const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
1422     #pragma omp parallel for
1423     for (dim_t i2=front; i2<m_N2; i2++) {
1424     for (dim_t i1=bottom; i1<m_N1; i1++) {
1425     for (dim_t i0=left; i0<m_N0; i0++) {
1426     m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
1427     +(i1-bottom)*(m_N0-left)
1428     +(i2-front)*(m_N0-left)*(m_N1-bottom);
1429     }
1430     }
1431     }
1432 caltinay 3735 m_nodeTags.assign(getNumNodes(), 0);
1433     updateTagsInUse(Nodes);
1434 caltinay 3698
1435 caltinay 3697 // elements
1436     m_elementId.resize(getNumElements());
1437     #pragma omp parallel for
1438     for (dim_t k=0; k<getNumElements(); k++) {
1439     m_elementId[k]=k;
1440     }
1441 caltinay 3735 m_elementTags.assign(getNumElements(), 0);
1442     updateTagsInUse(Elements);
1443 caltinay 3697
1444     // face elements
1445     m_faceId.resize(getNumFaceElements());
1446     #pragma omp parallel for
1447     for (dim_t k=0; k<getNumFaceElements(); k++) {
1448     m_faceId[k]=k;
1449     }
1450 caltinay 3704
1451 caltinay 3722 // generate face offset vector and set face tags
1452 caltinay 3704 const IndexVector facesPerEdge = getNumFacesPerBoundary();
1453 caltinay 3722 const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
1454     const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
1455 caltinay 3704 m_faceOffset.assign(facesPerEdge.size(), -1);
1456 caltinay 3722 m_faceTags.clear();
1457 caltinay 3704 index_t offset=0;
1458     for (size_t i=0; i<facesPerEdge.size(); i++) {
1459     if (facesPerEdge[i]>0) {
1460     m_faceOffset[i]=offset;
1461     offset+=facesPerEdge[i];
1462 caltinay 3722 m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1463 caltinay 3704 }
1464     }
1465 caltinay 3722 setTagMap("left", LEFT);
1466     setTagMap("right", RIGHT);
1467     setTagMap("bottom", BOTTOM);
1468     setTagMap("top", TOP);
1469     setTagMap("front", FRONT);
1470     setTagMap("back", BACK);
1471     updateTagsInUse(FaceElements);
1472 caltinay 3691 }
1473    
1474 caltinay 3703 //protected
1475 caltinay 3711 void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,
1476     bool reduced) const
1477 caltinay 3703 {
1478     const dim_t numComp = in.getDataPointSize();
1479 caltinay 3711 if (reduced) {
1480 caltinay 3731 /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1481     const double c0 = .125;
1482 caltinay 3703 #pragma omp parallel for
1483 caltinay 3711 for (index_t k2=0; k2 < m_NE2; ++k2) {
1484     for (index_t k1=0; k1 < m_NE1; ++k1) {
1485     for (index_t k0=0; k0 < m_NE0; ++k0) {
1486     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1487     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1488     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1489     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1490     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1491     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1492     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1493     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1494     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1495     for (index_t i=0; i < numComp; ++i) {
1496 caltinay 3731 o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i] + f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1497 caltinay 3711 } /* end of component loop i */
1498     } /* end of k0 loop */
1499     } /* end of k1 loop */
1500     } /* end of k2 loop */
1501     /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1502     } else {
1503 caltinay 3731 /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1504     const double c0 = .0094373878376559314545;
1505     const double c1 = .035220810900864519624;
1506     const double c2 = .13144585576580214704;
1507     const double c3 = .49056261216234406855;
1508 caltinay 3711 #pragma omp parallel for
1509     for (index_t k2=0; k2 < m_NE2; ++k2) {
1510     for (index_t k1=0; k1 < m_NE1; ++k1) {
1511     for (index_t k0=0; k0 < m_NE0; ++k0) {
1512     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1513     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1514     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1515     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1516     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1517     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1518     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1519     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1520     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1521     for (index_t i=0; i < numComp; ++i) {
1522 caltinay 3731 o[INDEX2(i,numComp,0)] = f_000[i]*c3 + f_111[i]*c0 + c2*(f_001[i] + f_010[i] + f_100[i]) + c1*(f_011[i] + f_101[i] + f_110[i]);
1523     o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_100[i]*c3 + c2*(f_000[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_111[i]);
1524     o[INDEX2(i,numComp,2)] = f_010[i]*c3 + f_101[i]*c0 + c2*(f_000[i] + f_011[i] + f_110[i]) + c1*(f_001[i] + f_100[i] + f_111[i]);
1525     o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_110[i]*c3 + c2*(f_010[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_101[i]);
1526     o[INDEX2(i,numComp,4)] = f_001[i]*c3 + f_110[i]*c0 + c2*(f_000[i] + f_011[i] + f_101[i]) + c1*(f_010[i] + f_100[i] + f_111[i]);
1527     o[INDEX2(i,numComp,5)] = f_010[i]*c0 + f_101[i]*c3 + c2*(f_001[i] + f_100[i] + f_111[i]) + c1*(f_000[i] + f_011[i] + f_110[i]);
1528     o[INDEX2(i,numComp,6)] = f_011[i]*c3 + f_100[i]*c0 + c2*(f_001[i] + f_010[i] + f_111[i]) + c1*(f_000[i] + f_101[i] + f_110[i]);
1529     o[INDEX2(i,numComp,7)] = f_000[i]*c0 + f_111[i]*c3 + c2*(f_011[i] + f_101[i] + f_110[i]) + c1*(f_001[i] + f_010[i] + f_100[i]);
1530 caltinay 3711 } /* end of component loop i */
1531     } /* end of k0 loop */
1532     } /* end of k1 loop */
1533     } /* end of k2 loop */
1534     /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1535     }
1536 caltinay 3703 }
1537    
1538     //protected
1539 caltinay 3711 void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1540     bool reduced) const
1541 caltinay 3703 {
1542 caltinay 3704 const dim_t numComp = in.getDataPointSize();
1543 caltinay 3711 if (reduced) {
1544 caltinay 3731 const double c0 = .25;
1545     #pragma omp parallel
1546     {
1547     /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1548     if (m_faceOffset[0] > -1) {
1549     #pragma omp for nowait
1550     for (index_t k2=0; k2 < m_NE2; ++k2) {
1551     for (index_t k1=0; k1 < m_NE1; ++k1) {
1552     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1553     const register double* f_010 = in