/[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 3748 - (hide annotations)
Thu Dec 15 07:36:19 2011 UTC (7 years, 4 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 93907 byte(s)
PDE assembly in serial and 2D seems to be doing what it's supposed to when
boundary elements are not involved.

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 3748 case ReducedNodes: //FIXME: reduced
241 caltinay 3697 return &m_nodeId[0];
242     case Elements:
243 caltinay 3733 case ReducedElements:
244 caltinay 3697 return &m_elementId[0];
245 caltinay 3733 case ReducedFaceElements:
246 caltinay 3697 return &m_faceId[0];
247     default:
248     break;
249     }
250 caltinay 3691
251 caltinay 3697 stringstream msg;
252     msg << "borrowSampleReferenceIDs() not implemented for function space type "
253     << fsType;
254     throw RipleyException(msg.str());
255 caltinay 3691 }
256    
257     bool Brick::ownSample(int fsCode, index_t id) const
258     {
259     #ifdef ESYS_MPI
260     if (fsCode == Nodes) {
261 caltinay 3698 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
262     const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
263 caltinay 3691 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
264     } else
265     throw RipleyException("ownSample() only implemented for Nodes");
266     #else
267     return true;
268     #endif
269     }
270    
271 caltinay 3703 void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const
272     {
273     escript::Data& in = *const_cast<escript::Data*>(&cIn);
274     const dim_t numComp = in.getDataPointSize();
275     const double h0 = m_l0/m_gNE0;
276     const double h1 = m_l1/m_gNE1;
277     const double h2 = m_l1/m_gNE2;
278 caltinay 3731 const double C0 = .044658198738520451079;
279     const double C1 = .16666666666666666667;
280     const double C2 = .21132486540518711775;
281     const double C3 = .25;
282     const double C4 = .5;
283     const double C5 = .62200846792814621559;
284     const double C6 = .78867513459481288225;
285    
286 caltinay 3703 if (out.getFunctionSpace().getTypeCode() == Elements) {
287 caltinay 3731 /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
288 caltinay 3703 #pragma omp parallel for
289 caltinay 3707 for (index_t k2=0; k2 < m_NE2; ++k2) {
290     for (index_t k1=0; k1 < m_NE1; ++k1) {
291     for (index_t k0=0; k0 < m_NE0; ++k0) {
292 caltinay 3703 const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
293     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
294     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
295     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
296     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
297     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
298     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
299     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
300     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
301     for (index_t i=0; i < numComp; ++i) {
302 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;
303     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;
304     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;
305     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;
306     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;
307     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;
308     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;
309     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;
310     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;
311     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;
312     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;
313     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;
314     o[INDEX3(i,0,0,numComp,3)] = V0;
315     o[INDEX3(i,1,0,numComp,3)] = V4;
316     o[INDEX3(i,2,0,numComp,3)] = V8;
317     o[INDEX3(i,0,1,numComp,3)] = V0;
318     o[INDEX3(i,1,1,numComp,3)] = V5;
319     o[INDEX3(i,2,1,numComp,3)] = V9;
320     o[INDEX3(i,0,2,numComp,3)] = V1;
321     o[INDEX3(i,1,2,numComp,3)] = V4;
322     o[INDEX3(i,2,2,numComp,3)] = V10;
323     o[INDEX3(i,0,3,numComp,3)] = V1;
324     o[INDEX3(i,1,3,numComp,3)] = V5;
325     o[INDEX3(i,2,3,numComp,3)] = V11;
326     o[INDEX3(i,0,4,numComp,3)] = V2;
327     o[INDEX3(i,1,4,numComp,3)] = V6;
328     o[INDEX3(i,2,4,numComp,3)] = V8;
329     o[INDEX3(i,0,5,numComp,3)] = V2;
330     o[INDEX3(i,1,5,numComp,3)] = V7;
331     o[INDEX3(i,2,5,numComp,3)] = V9;
332     o[INDEX3(i,0,6,numComp,3)] = V3;
333     o[INDEX3(i,1,6,numComp,3)] = V6;
334     o[INDEX3(i,2,6,numComp,3)] = V10;
335     o[INDEX3(i,0,7,numComp,3)] = V3;
336     o[INDEX3(i,1,7,numComp,3)] = V7;
337     o[INDEX3(i,2,7,numComp,3)] = V11;
338 caltinay 3703 } /* end of component loop i */
339     } /* end of k0 loop */
340     } /* end of k1 loop */
341     } /* end of k2 loop */
342     /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
343 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
344 caltinay 3731 /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
345 caltinay 3711 #pragma omp parallel for
346     for (index_t k2=0; k2 < m_NE2; ++k2) {
347     for (index_t k1=0; k1 < m_NE1; ++k1) {
348     for (index_t k0=0; k0 < m_NE0; ++k0) {
349     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
350     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
351     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
352     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
353     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
354     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
355     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
356     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
357     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
358     for (index_t i=0; i < numComp; ++i) {
359 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;
360     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;
361     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;
362 caltinay 3711 } /* end of component loop i */
363     } /* end of k0 loop */
364     } /* end of k1 loop */
365     } /* end of k2 loop */
366     /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
367 caltinay 3707 } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
368 caltinay 3731 /*** GENERATOR SNIP_GRAD_FACES TOP */
369 caltinay 3722 #pragma omp parallel
370     {
371     if (m_faceOffset[0] > -1) {
372     #pragma omp for nowait
373     for (index_t k2=0; k2 < m_NE2; ++k2) {
374     for (index_t k1=0; k1 < m_NE1; ++k1) {
375     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
376     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
377     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
378     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
379     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
380     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
381     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
382     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
383     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
384     for (index_t i=0; i < numComp; ++i) {
385 caltinay 3731 const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;
386     const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / h1;
387     const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / h2;
388     const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / h2;
389     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;
390     o[INDEX3(i,1,0,numComp,3)] = V0;
391     o[INDEX3(i,2,0,numComp,3)] = V2;
392     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;
393     o[INDEX3(i,1,1,numComp,3)] = V0;
394     o[INDEX3(i,2,1,numComp,3)] = V3;
395     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;
396     o[INDEX3(i,1,2,numComp,3)] = V1;
397     o[INDEX3(i,2,2,numComp,3)] = V2;
398     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;
399     o[INDEX3(i,1,3,numComp,3)] = V1;
400     o[INDEX3(i,2,3,numComp,3)] = V3;
401 caltinay 3722 } /* end of component loop i */
402     } /* end of k1 loop */
403     } /* end of k2 loop */
404     } /* end of face 0 */
405     if (m_faceOffset[1] > -1) {
406     #pragma omp for nowait
407     for (index_t k2=0; k2 < m_NE2; ++k2) {
408     for (index_t k1=0; k1 < m_NE1; ++k1) {
409     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
410     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
411     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
412     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
413     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
414     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
415     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
416     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
417     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
418     for (index_t i=0; i < numComp; ++i) {
419 caltinay 3731 const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
420     const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
421     const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
422     const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
423     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;
424     o[INDEX3(i,1,0,numComp,3)] = V0;
425     o[INDEX3(i,2,0,numComp,3)] = V2;
426     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;
427     o[INDEX3(i,1,1,numComp,3)] = V0;
428     o[INDEX3(i,2,1,numComp,3)] = V3;
429     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;
430     o[INDEX3(i,1,2,numComp,3)] = V1;
431     o[INDEX3(i,2,2,numComp,3)] = V2;
432     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;
433     o[INDEX3(i,1,3,numComp,3)] = V1;
434     o[INDEX3(i,2,3,numComp,3)] = V3;
435 caltinay 3722 } /* end of component loop i */
436     } /* end of k1 loop */
437     } /* end of k2 loop */
438     } /* end of face 1 */
439     if (m_faceOffset[2] > -1) {
440     #pragma omp for nowait
441     for (index_t k2=0; k2 < m_NE2; ++k2) {
442     for (index_t k0=0; k0 < m_NE0; ++k0) {
443     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
444     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
445     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
446     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
447     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
448     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
449     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
450     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
451     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
452     for (index_t i=0; i < numComp; ++i) {
453 caltinay 3731 const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;
454     const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / h2;
455     const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / h2;
456     o[INDEX3(i,0,0,numComp,3)] = V0;
457     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;
458     o[INDEX3(i,2,0,numComp,3)] = V1;
459     o[INDEX3(i,0,1,numComp,3)] = V0;
460     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;
461     o[INDEX3(i,2,1,numComp,3)] = V2;
462     o[INDEX3(i,0,2,numComp,3)] = V0;
463     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;
464     o[INDEX3(i,2,2,numComp,3)] = V1;
465     o[INDEX3(i,0,3,numComp,3)] = V0;
466     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;
467     o[INDEX3(i,2,3,numComp,3)] = V2;
468 caltinay 3722 } /* end of component loop i */
469     } /* end of k0 loop */
470     } /* end of k2 loop */
471     } /* end of face 2 */
472     if (m_faceOffset[3] > -1) {
473     #pragma omp for nowait
474     for (index_t k2=0; k2 < m_NE2; ++k2) {
475     for (index_t k0=0; k0 < m_NE0; ++k0) {
476     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
477     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
478     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
479     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
480     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
481     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
482     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
483     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
484     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
485     for (index_t i=0; i < numComp; ++i) {
486 caltinay 3731 const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
487     const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
488     const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
489     const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
490     o[INDEX3(i,0,0,numComp,3)] = V0;
491     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;
492     o[INDEX3(i,2,0,numComp,3)] = V2;
493     o[INDEX3(i,0,1,numComp,3)] = V0;
494     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;
495     o[INDEX3(i,2,1,numComp,3)] = V3;
496     o[INDEX3(i,0,2,numComp,3)] = V1;
497     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;
498     o[INDEX3(i,2,2,numComp,3)] = V2;
499     o[INDEX3(i,0,3,numComp,3)] = V1;
500     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;
501     o[INDEX3(i,2,3,numComp,3)] = V3;
502 caltinay 3722 } /* end of component loop i */
503     } /* end of k0 loop */
504     } /* end of k2 loop */
505     } /* end of face 3 */
506     if (m_faceOffset[4] > -1) {
507     #pragma omp for nowait
508 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
509 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
510     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
511     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
512     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
513     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
514     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
515     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
516     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
517     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
518     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
519     for (index_t i=0; i < numComp; ++i) {
520 caltinay 3731 const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;
521     const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / h0;
522     const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / h1;
523     const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / h1;
524     o[INDEX3(i,0,0,numComp,3)] = V0;
525     o[INDEX3(i,1,0,numComp,3)] = V2;
526     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;
527     o[INDEX3(i,0,1,numComp,3)] = V0;
528     o[INDEX3(i,1,1,numComp,3)] = V3;
529     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;
530     o[INDEX3(i,0,2,numComp,3)] = V1;
531     o[INDEX3(i,1,2,numComp,3)] = V2;
532     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;
533     o[INDEX3(i,0,3,numComp,3)] = V1;
534     o[INDEX3(i,1,3,numComp,3)] = V3;
535     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;
536 caltinay 3722 } /* end of component loop i */
537     } /* end of k0 loop */
538 caltinay 3707 } /* end of k1 loop */
539 caltinay 3722 } /* end of face 4 */
540     if (m_faceOffset[5] > -1) {
541     #pragma omp for nowait
542 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
543 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
544     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
545     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
546     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
547     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
548     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
549     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
550     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
551     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
552     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
553     for (index_t i=0; i < numComp; ++i) {
554 caltinay 3731 const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
555     const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
556     const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
557     const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
558     o[INDEX3(i,0,0,numComp,3)] = V0;
559     o[INDEX3(i,1,0,numComp,3)] = V2;
560     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;
561     o[INDEX3(i,0,1,numComp,3)] = V0;
562     o[INDEX3(i,1,1,numComp,3)] = V3;
563     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;
564     o[INDEX3(i,0,2,numComp,3)] = V1;
565     o[INDEX3(i,1,2,numComp,3)] = V2;
566     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;
567     o[INDEX3(i,0,3,numComp,3)] = V1;
568     o[INDEX3(i,1,3,numComp,3)] = V3;
569     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;
570 caltinay 3722 } /* end of component loop i */
571     } /* end of k0 loop */
572 caltinay 3707 } /* end of k1 loop */
573 caltinay 3722 } /* end of face 5 */
574     } // end of parallel section
575 caltinay 3707 /* GENERATOR SNIP_GRAD_FACES BOTTOM */
576 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
577 caltinay 3731 /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
578 caltinay 3722 #pragma omp parallel
579     {
580     if (m_faceOffset[0] > -1) {
581     #pragma omp for nowait
582     for (index_t k2=0; k2 < m_NE2; ++k2) {
583     for (index_t k1=0; k1 < m_NE1; ++k1) {
584     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
585     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
586     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
587     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
588     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
589     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
590     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
591     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
592     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
593     for (index_t i=0; i < numComp; ++i) {
594 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;
595     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;
596     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;
597 caltinay 3722 } /* end of component loop i */
598     } /* end of k1 loop */
599     } /* end of k2 loop */
600     } /* end of face 0 */
601     if (m_faceOffset[1] > -1) {
602     #pragma omp for nowait
603     for (index_t k2=0; k2 < m_NE2; ++k2) {
604     for (index_t k1=0; k1 < m_NE1; ++k1) {
605     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
606     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
607     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
608     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
609     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
610     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
611     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
612     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
613     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
614     for (index_t i=0; i < numComp; ++i) {
615 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;
616     o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;
617     o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;
618 caltinay 3722 } /* end of component loop i */
619     } /* end of k1 loop */
620     } /* end of k2 loop */
621     } /* end of face 1 */
622     if (m_faceOffset[2] > -1) {
623     #pragma omp for nowait
624     for (index_t k2=0; k2 < m_NE2; ++k2) {
625     for (index_t k0=0; k0 < m_NE0; ++k0) {
626     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
627     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
628     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
629     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
630     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
631     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
632     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
633     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
634     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
635     for (index_t i=0; i < numComp; ++i) {
636 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
637     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;
638     o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;
639 caltinay 3722 } /* end of component loop i */
640     } /* end of k0 loop */
641     } /* end of k2 loop */
642     } /* end of face 2 */
643     if (m_faceOffset[3] > -1) {
644     #pragma omp for nowait
645     for (index_t k2=0; k2 < m_NE2; ++k2) {
646     for (index_t k0=0; k0 < m_NE0; ++k0) {
647     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
648     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
649     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
650     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
651     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
652     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
653     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
654     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
655     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
656     for (index_t i=0; i < numComp; ++i) {
657 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
658     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;
659     o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;
660 caltinay 3722 } /* end of component loop i */
661     } /* end of k0 loop */
662     } /* end of k2 loop */
663     } /* end of face 3 */
664     if (m_faceOffset[4] > -1) {
665     #pragma omp for nowait
666 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
667 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
668     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
669     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
670     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
671     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
672     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
673     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
674     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
675     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
676     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
677     for (index_t i=0; i < numComp; ++i) {
678 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
679     o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;
680     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;
681 caltinay 3722 } /* end of component loop i */
682     } /* end of k0 loop */
683 caltinay 3711 } /* end of k1 loop */
684 caltinay 3722 } /* end of face 4 */
685     if (m_faceOffset[5] > -1) {
686     #pragma omp for nowait
687 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
688 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
689     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
690     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
691     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
692     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
693     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
694     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
695     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
696     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
697     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
698     for (index_t i=0; i < numComp; ++i) {
699 caltinay 3731 o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
700     o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;
701     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;
702 caltinay 3722 } /* end of component loop i */
703     } /* end of k0 loop */
704 caltinay 3711 } /* end of k1 loop */
705 caltinay 3722 } /* end of face 5 */
706     } // end of parallel section
707 caltinay 3711 /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
708 caltinay 3703 } else {
709 caltinay 3707 stringstream msg;
710     msg << "setToGradient() not implemented for "
711     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
712     throw RipleyException(msg.str());
713 caltinay 3703 }
714     }
715    
716 caltinay 3713 void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
717     {
718     escript::Data& in = *const_cast<escript::Data*>(&arg);
719     const dim_t numComp = in.getDataPointSize();
720     const double h0 = m_l0/m_gNE0;
721     const double h1 = m_l1/m_gNE1;
722     const double h2 = m_l2/m_gNE2;
723     if (arg.getFunctionSpace().getTypeCode() == Elements) {
724     const double w_0 = h0*h1*h2/8.;
725     #pragma omp parallel
726     {
727     vector<double> int_local(numComp, 0);
728 caltinay 3722 #pragma omp for nowait
729 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
730     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
731     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
732     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
733     for (index_t i=0; i < numComp; ++i) {
734     const register double f_0 = f[INDEX2(i,0,numComp)];
735     const register double f_1 = f[INDEX2(i,1,numComp)];
736     const register double f_2 = f[INDEX2(i,2,numComp)];
737     const register double f_3 = f[INDEX2(i,3,numComp)];
738     const register double f_4 = f[INDEX2(i,4,numComp)];
739     const register double f_5 = f[INDEX2(i,5,numComp)];
740     const register double f_6 = f[INDEX2(i,6,numComp)];
741     const register double f_7 = f[INDEX2(i,7,numComp)];
742     int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
743     } /* end of component loop i */
744     } /* end of k0 loop */
745     } /* end of k1 loop */
746     } /* end of k2 loop */
747    
748     #pragma omp critical
749     for (index_t i=0; i<numComp; i++)
750     integrals[i]+=int_local[i];
751 caltinay 3722 } // end of parallel section
752 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
753     const double w_0 = h0*h1*h2;
754     #pragma omp parallel
755     {
756     vector<double> int_local(numComp, 0);
757 caltinay 3722 #pragma omp for nowait
758 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
759     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
760     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
761     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
762     for (index_t i=0; i < numComp; ++i) {
763     int_local[i]+=f[i]*w_0;
764     } /* end of component loop i */
765     } /* end of k0 loop */
766     } /* end of k1 loop */
767     } /* end of k2 loop */
768    
769     #pragma omp critical
770     for (index_t i=0; i<numComp; i++)
771     integrals[i]+=int_local[i];
772 caltinay 3722 } // end of parallel section
773 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
774     const double w_0 = h1*h2/4.;
775     const double w_1 = h0*h2/4.;
776     const double w_2 = h0*h1/4.;
777     #pragma omp parallel
778     {
779     vector<double> int_local(numComp, 0);
780     if (m_faceOffset[0] > -1) {
781 caltinay 3722 #pragma omp for nowait
782 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
783     for (index_t k1=0; k1 < m_NE1; ++k1) {
784     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
785     for (index_t i=0; i < numComp; ++i) {
786     const register double f_0 = f[INDEX2(i,0,numComp)];
787     const register double f_1 = f[INDEX2(i,1,numComp)];
788     const register double f_2 = f[INDEX2(i,2,numComp)];
789     const register double f_3 = f[INDEX2(i,3,numComp)];
790     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
791     } /* end of component loop i */
792     } /* end of k1 loop */
793     } /* end of k2 loop */
794     }
795    
796     if (m_faceOffset[1] > -1) {
797 caltinay 3722 #pragma omp for nowait
798 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
799     for (index_t k1=0; k1 < m_NE1; ++k1) {
800     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
801     for (index_t i=0; i < numComp; ++i) {
802     const register double f_0 = f[INDEX2(i,0,numComp)];
803     const register double f_1 = f[INDEX2(i,1,numComp)];
804     const register double f_2 = f[INDEX2(i,2,numComp)];
805     const register double f_3 = f[INDEX2(i,3,numComp)];
806     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
807     } /* end of component loop i */
808     } /* end of k1 loop */
809     } /* end of k2 loop */
810     }
811    
812     if (m_faceOffset[2] > -1) {
813 caltinay 3722 #pragma omp for nowait
814 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
815     for (index_t k0=0; k0 < m_NE0; ++k0) {
816     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
817     for (index_t i=0; i < numComp; ++i) {
818     const register double f_0 = f[INDEX2(i,0,numComp)];
819     const register double f_1 = f[INDEX2(i,1,numComp)];
820     const register double f_2 = f[INDEX2(i,2,numComp)];
821     const register double f_3 = f[INDEX2(i,3,numComp)];
822     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
823     } /* end of component loop i */
824     } /* end of k1 loop */
825     } /* end of k2 loop */
826     }
827    
828     if (m_faceOffset[3] > -1) {
829 caltinay 3722 #pragma omp for nowait
830 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
831     for (index_t k0=0; k0 < m_NE0; ++k0) {
832     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
833     for (index_t i=0; i < numComp; ++i) {
834     const register double f_0 = f[INDEX2(i,0,numComp)];
835     const register double f_1 = f[INDEX2(i,1,numComp)];
836     const register double f_2 = f[INDEX2(i,2,numComp)];
837     const register double f_3 = f[INDEX2(i,3,numComp)];
838     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
839     } /* end of component loop i */
840     } /* end of k1 loop */
841     } /* end of k2 loop */
842     }
843    
844     if (m_faceOffset[4] > -1) {
845 caltinay 3722 #pragma omp for nowait
846 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
847     for (index_t k0=0; k0 < m_NE0; ++k0) {
848     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
849     for (index_t i=0; i < numComp; ++i) {
850     const register double f_0 = f[INDEX2(i,0,numComp)];
851     const register double f_1 = f[INDEX2(i,1,numComp)];
852     const register double f_2 = f[INDEX2(i,2,numComp)];
853     const register double f_3 = f[INDEX2(i,3,numComp)];
854     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
855     } /* end of component loop i */
856     } /* end of k1 loop */
857     } /* end of k2 loop */
858     }
859    
860     if (m_faceOffset[5] > -1) {
861 caltinay 3722 #pragma omp for nowait
862 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
863     for (index_t k0=0; k0 < m_NE0; ++k0) {
864     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
865     for (index_t i=0; i < numComp; ++i) {
866     const register double f_0 = f[INDEX2(i,0,numComp)];
867     const register double f_1 = f[INDEX2(i,1,numComp)];
868     const register double f_2 = f[INDEX2(i,2,numComp)];
869     const register double f_3 = f[INDEX2(i,3,numComp)];
870     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
871     } /* end of component loop i */
872     } /* end of k1 loop */
873     } /* end of k2 loop */
874     }
875    
876     #pragma omp critical
877     for (index_t i=0; i<numComp; i++)
878     integrals[i]+=int_local[i];
879 caltinay 3722 } // end of parallel section
880 caltinay 3713
881     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
882     const double w_0 = h1*h2;
883     const double w_1 = h0*h2;
884     const double w_2 = h0*h1;
885     #pragma omp parallel
886     {
887     vector<double> int_local(numComp, 0);
888     if (m_faceOffset[0] > -1) {
889 caltinay 3722 #pragma omp for nowait
890 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
891     for (index_t k1=0; k1 < m_NE1; ++k1) {
892     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
893     for (index_t i=0; i < numComp; ++i) {
894     int_local[i]+=f[i]*w_0;
895     } /* end of component loop i */
896     } /* end of k1 loop */
897     } /* end of k2 loop */
898     }
899    
900     if (m_faceOffset[1] > -1) {
901 caltinay 3722 #pragma omp for nowait
902 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
903     for (index_t k1=0; k1 < m_NE1; ++k1) {
904     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
905     for (index_t i=0; i < numComp; ++i) {
906     int_local[i]+=f[i]*w_0;
907     } /* end of component loop i */
908     } /* end of k1 loop */
909     } /* end of k2 loop */
910     }
911    
912     if (m_faceOffset[2] > -1) {
913 caltinay 3722 #pragma omp for nowait
914 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
915     for (index_t k0=0; k0 < m_NE0; ++k0) {
916     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
917     for (index_t i=0; i < numComp; ++i) {
918     int_local[i]+=f[i]*w_1;
919     } /* end of component loop i */
920     } /* end of k1 loop */
921     } /* end of k2 loop */
922     }
923    
924     if (m_faceOffset[3] > -1) {
925 caltinay 3722 #pragma omp for nowait
926 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
927     for (index_t k0=0; k0 < m_NE0; ++k0) {
928     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
929     for (index_t i=0; i < numComp; ++i) {
930     int_local[i]+=f[i]*w_1;
931     } /* end of component loop i */
932     } /* end of k1 loop */
933     } /* end of k2 loop */
934     }
935    
936     if (m_faceOffset[4] > -1) {
937 caltinay 3722 #pragma omp for nowait
938 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
939     for (index_t k0=0; k0 < m_NE0; ++k0) {
940     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
941     for (index_t i=0; i < numComp; ++i) {
942     int_local[i]+=f[i]*w_2;
943     } /* end of component loop i */
944     } /* end of k1 loop */
945     } /* end of k2 loop */
946     }
947    
948     if (m_faceOffset[5] > -1) {
949 caltinay 3722 #pragma omp for nowait
950 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
951     for (index_t k0=0; k0 < m_NE0; ++k0) {
952     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
953     for (index_t i=0; i < numComp; ++i) {
954     int_local[i]+=f[i]*w_2;
955     } /* end of component loop i */
956     } /* end of k1 loop */
957     } /* end of k2 loop */
958     }
959    
960     #pragma omp critical
961     for (index_t i=0; i<numComp; i++)
962     integrals[i]+=int_local[i];
963 caltinay 3722 } // end of parallel section
964 caltinay 3713
965     } else {
966     stringstream msg;
967     msg << "setToIntegrals() not implemented for "
968     << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
969     throw RipleyException(msg.str());
970     }
971     }
972    
973 caltinay 3722 void Brick::setToNormal(escript::Data& out) const
974     {
975     if (out.getFunctionSpace().getTypeCode() == FaceElements) {
976     #pragma omp parallel
977     {
978     if (m_faceOffset[0] > -1) {
979     #pragma omp for nowait
980     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
981     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
982     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
983     // set vector at four quadrature points
984     *o++ = -1.; *o++ = 0.; *o++ = 0.;
985     *o++ = -1.; *o++ = 0.; *o++ = 0.;
986     *o++ = -1.; *o++ = 0.; *o++ = 0.;
987     *o++ = -1.; *o++ = 0.; *o = 0.;
988     }
989     }
990     }
991    
992     if (m_faceOffset[1] > -1) {
993     #pragma omp for nowait
994     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
995     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
996     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
997     // set vector at four quadrature points
998     *o++ = 1.; *o++ = 0.; *o++ = 0.;
999     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1000     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1001     *o++ = 1.; *o++ = 0.; *o = 0.;
1002     }
1003     }
1004     }
1005    
1006     if (m_faceOffset[2] > -1) {
1007     #pragma omp for nowait
1008     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1009     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1010     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1011     // set vector at four quadrature points
1012     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1013     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1014     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1015     *o++ = 0.; *o++ = -1.; *o = 0.;
1016     }
1017     }
1018     }
1019    
1020     if (m_faceOffset[3] > -1) {
1021     #pragma omp for nowait
1022     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1023     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1024     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1025     // set vector at four quadrature points
1026     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1027     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1028     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1029     *o++ = 0.; *o++ = 1.; *o = 0.;
1030     }
1031     }
1032     }
1033    
1034     if (m_faceOffset[4] > -1) {
1035     #pragma omp for nowait
1036     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1037     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1038     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1039     // set vector at four quadrature points
1040     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1041     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1042     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1043     *o++ = 0.; *o++ = 0.; *o = -1.;
1044     }
1045     }
1046     }
1047    
1048     if (m_faceOffset[5] > -1) {
1049     #pragma omp for nowait
1050     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1051     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1052     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1053     // set vector at four quadrature points
1054     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1055     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1056     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1057     *o++ = 0.; *o++ = 0.; *o = 1.;
1058     }
1059     }
1060     }
1061     } // end of parallel section
1062     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1063     #pragma omp parallel
1064     {
1065     if (m_faceOffset[0] > -1) {
1066     #pragma omp for nowait
1067     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1068     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1069     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1070     *o++ = -1.;
1071     *o++ = 0.;
1072     *o = 0.;
1073     }
1074     }
1075     }
1076    
1077     if (m_faceOffset[1] > -1) {
1078     #pragma omp for nowait
1079     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1080     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1081     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1082     *o++ = 1.;
1083     *o++ = 0.;
1084     *o = 0.;
1085     }
1086     }
1087     }
1088    
1089     if (m_faceOffset[2] > -1) {
1090     #pragma omp for nowait
1091     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1092     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1093     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1094     *o++ = 0.;
1095     *o++ = -1.;
1096     *o = 0.;
1097     }
1098     }
1099     }
1100    
1101     if (m_faceOffset[3] > -1) {
1102     #pragma omp for nowait
1103     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1104     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1105     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1106     *o++ = 0.;
1107     *o++ = 1.;
1108     *o = 0.;
1109     }
1110     }
1111     }
1112    
1113     if (m_faceOffset[4] > -1) {
1114     #pragma omp for nowait
1115     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1116     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1117     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1118     *o++ = 0.;
1119     *o++ = 0.;
1120     *o = -1.;
1121     }
1122     }
1123     }
1124    
1125     if (m_faceOffset[5] > -1) {
1126     #pragma omp for nowait
1127     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1128     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1129     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1130     *o++ = 0.;
1131     *o++ = 0.;
1132     *o = 1.;
1133     }
1134     }
1135     }
1136     } // end of parallel section
1137    
1138     } else {
1139     stringstream msg;
1140     msg << "setToNormal() not implemented for "
1141     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
1142     throw RipleyException(msg.str());
1143     }
1144     }
1145    
1146 caltinay 3691 Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
1147     bool reducedColOrder) const
1148     {
1149     if (reducedRowOrder || reducedColOrder)
1150     throw RipleyException("getPattern() not implemented for reduced order");
1151    
1152     throw RipleyException("getPattern() not implemented");
1153     }
1154    
1155     void Brick::Print_Mesh_Info(const bool full) const
1156     {
1157     RipleyDomain::Print_Mesh_Info(full);
1158     if (full) {
1159     cout << " Id Coordinates" << endl;
1160     cout.precision(15);
1161     cout.setf(ios::scientific, ios::floatfield);
1162 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
1163     pair<double,double> ydy = getFirstCoordAndSpacing(1);
1164     pair<double,double> zdz = getFirstCoordAndSpacing(2);
1165 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
1166     cout << " " << setw(5) << m_nodeId[i]
1167 caltinay 3698 << " " << xdx.first+(i%m_N0)*xdx.second
1168     << " " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
1169     << " " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
1170 caltinay 3691 }
1171     }
1172     }
1173    
1174 caltinay 3698 IndexVector Brick::getNumNodesPerDim() const
1175     {
1176     IndexVector ret;
1177     ret.push_back(m_N0);
1178     ret.push_back(m_N1);
1179     ret.push_back(m_N2);
1180     return ret;
1181     }
1182    
1183     IndexVector Brick::getNumElementsPerDim() const
1184     {
1185     IndexVector ret;
1186     ret.push_back(m_NE0);
1187     ret.push_back(m_NE1);
1188     ret.push_back(m_NE2);
1189     return ret;
1190     }
1191    
1192     IndexVector Brick::getNumFacesPerBoundary() const
1193     {
1194     IndexVector ret(6, 0);
1195     //left
1196     if (m_offset0==0)
1197     ret[0]=m_NE1*m_NE2;
1198     //right
1199     if (m_mpiInfo->rank%m_NX==m_NX-1)
1200     ret[1]=m_NE1*m_NE2;
1201     //bottom
1202     if (m_offset1==0)
1203     ret[2]=m_NE0*m_NE2;
1204     //top
1205     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1206     ret[3]=m_NE0*m_NE2;
1207     //front
1208     if (m_offset2==0)
1209     ret[4]=m_NE0*m_NE1;
1210     //back
1211     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1212     ret[5]=m_NE0*m_NE1;
1213     return ret;
1214     }
1215    
1216     pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
1217     {
1218     if (dim==0)
1219     return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
1220     else if (dim==1)
1221     return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
1222     else if (dim==2)
1223     return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
1224    
1225     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
1226     }
1227    
1228 caltinay 3744 //protected
1229 caltinay 3691 dim_t Brick::getNumFaceElements() const
1230     {
1231     dim_t n=0;
1232     //left
1233     if (m_offset0==0)
1234     n+=m_NE1*m_NE2;
1235     //right
1236     if (m_mpiInfo->rank%m_NX==m_NX-1)
1237     n+=m_NE1*m_NE2;
1238     //bottom
1239     if (m_offset1==0)
1240     n+=m_NE0*m_NE2;
1241     //top
1242     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1243     n+=m_NE0*m_NE2;
1244     //front
1245     if (m_offset2==0)
1246     n+=m_NE0*m_NE1;
1247     //back
1248     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1249     n+=m_NE0*m_NE1;
1250    
1251     return n;
1252     }
1253    
1254     //protected
1255     void Brick::assembleCoordinates(escript::Data& arg) const
1256     {
1257     escriptDataC x = arg.getDataC();
1258     int numDim = m_numDim;
1259     if (!isDataPointShapeEqual(&x, 1, &numDim))
1260     throw RipleyException("setToX: Invalid Data object shape");
1261     if (!numSamplesEqual(&x, 1, getNumNodes()))
1262     throw RipleyException("setToX: Illegal number of samples in Data object");
1263    
1264 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
1265     pair<double,double> ydy = getFirstCoordAndSpacing(1);
1266     pair<double,double> zdz = getFirstCoordAndSpacing(2);
1267 caltinay 3691 arg.requireWrite();
1268     #pragma omp parallel for
1269     for (dim_t i2 = 0; i2 < m_N2; i2++) {
1270     for (dim_t i1 = 0; i1 < m_N1; i1++) {
1271     for (dim_t i0 = 0; i0 < m_N0; i0++) {
1272     double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
1273 caltinay 3698 point[0] = xdx.first+i0*xdx.second;
1274     point[1] = ydy.first+i1*ydy.second;
1275     point[2] = zdz.first+i2*zdz.second;
1276 caltinay 3691 }
1277     }
1278     }
1279     }
1280    
1281     //private
1282     void Brick::populateSampleIds()
1283     {
1284 caltinay 3697 // identifiers are ordered from left to right, bottom to top, front to back
1285     // on each rank, except for the shared nodes which are owned by the rank
1286     // below / to the left / to the front of the current rank
1287 caltinay 3698
1288     // build node distribution vector first.
1289     // m_nodeDistribution[i] is the first node id on rank i, that is
1290     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
1291     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1292     m_nodeDistribution[1]=getNumNodes();
1293     for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
1294     const index_t x = k%m_NX;
1295     const index_t y = k%(m_NX*m_NY)/m_NX;
1296     const index_t z = k/(m_NX*m_NY);
1297     index_t numNodes=getNumNodes();
1298     if (x>0)
1299     numNodes-=m_N1*m_N2;
1300     if (y>0)
1301     numNodes-=m_N0*m_N2;
1302     if (z>0)
1303     numNodes-=m_N0*m_N1;
1304     // if an edge was subtracted twice add it back
1305     if (x>0 && y>0)
1306     numNodes+=m_N2;
1307     if (x>0 && z>0)
1308     numNodes+=m_N1;
1309     if (y>0 && z>0)
1310     numNodes+=m_N0;
1311     // the corner node was removed 3x and added back 3x, so subtract it
1312     if (x>0 && y>0 && z>0)
1313     numNodes--;
1314     m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
1315     }
1316     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
1317    
1318 caltinay 3691 m_nodeId.resize(getNumNodes());
1319 caltinay 3698
1320     // the bottom, left and front planes are not owned by this rank so the
1321     // identifiers need to be computed accordingly
1322     const index_t left = (m_offset0==0 ? 0 : 1);
1323     const index_t bottom = (m_offset1==0 ? 0 : 1);
1324     const index_t front = (m_offset2==0 ? 0 : 1);
1325    
1326     // case 1: all nodes on left plane are owned by rank on the left
1327     if (left>0) {
1328     const int neighbour=m_mpiInfo->rank-1;
1329     const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1330     const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1331 caltinay 3691 #pragma omp parallel for
1332 caltinay 3698 for (dim_t i2=front; i2<m_N2; i2++) {
1333     for (dim_t i1=bottom; i1<m_N1; i1++) {
1334     m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1335     + (i1-bottom+1)*leftN0
1336     + (i2-front)*leftN0*leftN1 - 1;
1337     }
1338     }
1339 caltinay 3691 }
1340 caltinay 3698 // case 2: all nodes on bottom plane are owned by rank below
1341     if (bottom>0) {
1342     const int neighbour=m_mpiInfo->rank-m_NX;
1343     const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1344     const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1345     #pragma omp parallel for
1346     for (dim_t i2=front; i2<m_N2; i2++) {
1347     for (dim_t i0=left; i0<m_N0; i0++) {
1348     m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1349     + bottomN0*(bottomN1-1)
1350     + (i2-front)*bottomN0*bottomN1 + i0-left;
1351     }
1352     }
1353     }
1354     // case 3: all nodes on front plane are owned by rank in front
1355     if (front>0) {
1356     const int neighbour=m_mpiInfo->rank-m_NX*m_NY;
1357     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1358     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1359     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1360     #pragma omp parallel for
1361     for (dim_t i1=bottom; i1<m_N1; i1++) {
1362     for (dim_t i0=left; i0<m_N0; i0++) {
1363     m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]
1364     + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;
1365     }
1366     }
1367     }
1368     // case 4: nodes on front bottom edge are owned by the corresponding rank
1369     if (front>0 && bottom>0) {
1370     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);
1371     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1372     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1373     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1374     #pragma omp parallel for
1375     for (dim_t i0=left; i0<m_N0; i0++) {
1376     m_nodeId[i0]=m_nodeDistribution[neighbour]
1377     + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;
1378     }
1379     }
1380     // case 5: nodes on left bottom edge are owned by the corresponding rank
1381     if (left>0 && bottom>0) {
1382     const int neighbour=m_mpiInfo->rank-m_NX-1;
1383     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1384     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1385     #pragma omp parallel for
1386     for (dim_t i2=front; i2<m_N2; i2++) {
1387     m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
1388     + (1+i2-front)*N0*N1-1;
1389     }
1390     }
1391     // case 6: nodes on left front edge are owned by the corresponding rank
1392     if (left>0 && front>0) {
1393     const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;
1394     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1395     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1396     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
1397     #pragma omp parallel for
1398     for (dim_t i1=bottom; i1<m_N1; i1++) {
1399     m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
1400     + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;
1401     }
1402     }
1403     // case 7: bottom-left-front corner node owned by corresponding rank
1404     if (left>0 && bottom>0 && front>0) {
1405     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;
1406     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
1407     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
1408     const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);
1409     m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;
1410     }
1411 caltinay 3697
1412 caltinay 3698 // the rest of the id's are contiguous
1413     const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
1414     #pragma omp parallel for
1415     for (dim_t i2=front; i2<m_N2; i2++) {
1416     for (dim_t i1=bottom; i1<m_N1; i1++) {
1417     for (dim_t i0=left; i0<m_N0; i0++) {
1418     m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
1419     +(i1-bottom)*(m_N0-left)
1420     +(i2-front)*(m_N0-left)*(m_N1-bottom);
1421     }
1422     }
1423     }
1424 caltinay 3735 m_nodeTags.assign(getNumNodes(), 0);
1425     updateTagsInUse(Nodes);
1426 caltinay 3698
1427 caltinay 3697 // elements
1428     m_elementId.resize(getNumElements());
1429     #pragma omp parallel for
1430     for (dim_t k=0; k<getNumElements(); k++) {
1431     m_elementId[k]=k;
1432     }
1433 caltinay 3735 m_elementTags.assign(getNumElements(), 0);
1434     updateTagsInUse(Elements);
1435 caltinay 3697
1436     // face elements
1437     m_faceId.resize(getNumFaceElements());
1438     #pragma omp parallel for
1439     for (dim_t k=0; k<getNumFaceElements(); k++) {
1440     m_faceId[k]=k;
1441     }
1442 caltinay 3704
1443 caltinay 3722 // generate face offset vector and set face tags
1444 caltinay 3704 const IndexVector facesPerEdge = getNumFacesPerBoundary();
1445 caltinay 3722 const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
1446     const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
1447 caltinay 3704 m_faceOffset.assign(facesPerEdge.size(), -1);
1448 caltinay 3722 m_faceTags.clear();
1449 caltinay 3704 index_t offset=0;
1450     for (size_t i=0; i<facesPerEdge.size(); i++) {
1451     if (facesPerEdge[i]>0) {
1452     m_faceOffset[i]=offset;
1453     offset+=facesPerEdge[i];
1454 caltinay 3722 m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1455 caltinay 3704 }
1456     }
1457 caltinay 3722 setTagMap("left", LEFT);
1458     setTagMap("right", RIGHT);
1459     setTagMap("bottom", BOTTOM);
1460     setTagMap("top", TOP);
1461     setTagMap("front", FRONT);
1462     setTagMap("back", BACK);
1463     updateTagsInUse(FaceElements);
1464 caltinay 3691 }
1465    
1466 caltinay 3703 //protected
1467 caltinay 3711 void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,
1468     bool reduced) const
1469 caltinay 3703 {
1470     const dim_t numComp = in.getDataPointSize();
1471 caltinay 3711 if (reduced) {
1472 caltinay 3731 /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1473     const double c0 = .125;
1474 caltinay 3703 #pragma omp parallel for
1475 caltinay 3711 for (index_t k2=0; k2 < m_NE2; ++k2) {
1476     for (index_t k1=0; k1 < m_NE1; ++k1) {
1477     for (index_t k0=0; k0 < m_NE0; ++k0) {
1478     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1479     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1480     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1481     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1482     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1483     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1484     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1485     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1486     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1487     for (index_t i=0; i < numComp; ++i) {
1488 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]);
1489 caltinay 3711 } /* end of component loop i */
1490     } /* end of k0 loop */
1491     } /* end of k1 loop */
1492     } /* end of k2 loop */
1493     /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1494     } else {
1495 caltinay 3731 /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1496     const double c0 = .0094373878376559314545;
1497     const double c1 = .035220810900864519624;
1498     const double c2 = .13144585576580214704;
1499     const double c3 = .49056261216234406855;
1500 caltinay 3711 #pragma omp parallel for
1501     for (index_t k2=0; k2 < m_NE2; ++k2) {
1502     for (index_t k1=0; k1 < m_NE1; ++k1) {
1503     for (index_t k0=0; k0 < m_NE0; ++k0) {
1504     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1505     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1506     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1507     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1508     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1509     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1510     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1511     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1512     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1513     for (index_t i=0; i < numComp; ++i) {
1514 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]);
1515     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]);
1516     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]);
1517     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]);
1518     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]);
1519     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]);
1520     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]);
1521     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]);
1522 caltinay 3711 } /* end of component loop i */
1523     } /* end of k0 loop */
1524     } /* end of k1 loop */
1525     } /* end of k2 loop */
1526     /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1527     }
1528 caltinay 3703 }
1529    
1530     //protected
1531 caltinay 3711 void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1532     bool reduced) const
1533 caltinay 3703 {
1534 caltinay 3704 const dim_t numComp = in.getDataPointSize();
1535 caltinay 3711 if (reduced) {
1536 caltinay 3731 const double c0 = .25;
1537     #pragma omp parallel
1538     {
1539     /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1540     if (m_faceOffset[0] > -1) {
1541     #pragma omp for nowait
1542     for (index_t k2=0; k2 < m_NE2; ++k2) {
1543     for (index_t k1=0; k1 < m_NE1; ++k1) {
1544     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1545     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1546     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1547     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1548     double* o = out