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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3722 - (hide annotations)
Wed Dec 7 05:53:22 2011 UTC (7 years, 4 months ago) by caltinay
File size: 117887 byte(s)
All "util on ripley" tests pass now :-)
-added support for node/element/face tags
-implemented setToNormal()
-added a bunch of omp nowait

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     if (dynamic_cast<const Brick*>(&other))
81     return this==&other;
82    
83     return false;
84     }
85    
86     void Brick::dump(const string& fileName) const
87     {
88     #if USE_SILO
89     string fn(fileName);
90     if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
91     fn+=".silo";
92     }
93    
94     const int NUM_SILO_FILES = 1;
95     const char* blockDirFmt = "/block%04d";
96     int driver=DB_HDF5;
97     string siloPath;
98     DBfile* dbfile = NULL;
99    
100     #ifdef ESYS_MPI
101     PMPIO_baton_t* baton = NULL;
102     #endif
103    
104     if (m_mpiInfo->size > 1) {
105     #ifdef ESYS_MPI
106     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
107     0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
108     PMPIO_DefaultClose, (void*)&driver);
109     // try the fallback driver in case of error
110     if (!baton && driver != DB_PDB) {
111     driver = DB_PDB;
112     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
113     0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
114     PMPIO_DefaultClose, (void*)&driver);
115     }
116     if (baton) {
117     char str[64];
118     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
119     siloPath = str;
120     dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
121     }
122     #endif
123     } else {
124     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
125     getDescription().c_str(), driver);
126     // try the fallback driver in case of error
127     if (!dbfile && driver != DB_PDB) {
128     driver = DB_PDB;
129     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
130     getDescription().c_str(), driver);
131     }
132     }
133    
134     if (!dbfile)
135     throw RipleyException("dump: Could not create Silo file");
136    
137     /*
138     if (driver==DB_HDF5) {
139     // gzip level 1 already provides good compression with minimal
140     // performance penalty. Some tests showed that gzip levels >3 performed
141     // rather badly on escript data both in terms of time and space
142     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
143     }
144     */
145    
146     boost::scoped_ptr<double> x(new double[m_N0]);
147     boost::scoped_ptr<double> y(new double[m_N1]);
148     boost::scoped_ptr<double> z(new double[m_N2]);
149     double* coords[3] = { x.get(), y.get(), z.get() };
150 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
151     pair<double,double> ydy = getFirstCoordAndSpacing(1);
152     pair<double,double> zdz = getFirstCoordAndSpacing(2);
153 caltinay 3691 #pragma omp parallel
154     {
155     #pragma omp for
156     for (dim_t i0 = 0; i0 < m_N0; i0++) {
157 caltinay 3698 coords[0][i0]=xdx.first+i0*xdx.second;
158 caltinay 3691 }
159     #pragma omp for
160     for (dim_t i1 = 0; i1 < m_N1; i1++) {
161 caltinay 3698 coords[1][i1]=ydy.first+i1*ydy.second;
162 caltinay 3691 }
163     #pragma omp for
164     for (dim_t i2 = 0; i2 < m_N2; i2++) {
165 caltinay 3698 coords[2][i2]=zdz.first+i2*zdz.second;
166 caltinay 3691 }
167     }
168 caltinay 3698 IndexVector dims = getNumNodesPerDim();
169     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,
170 caltinay 3691 DB_COLLINEAR, NULL);
171    
172 caltinay 3698 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
173     NULL, 0, DB_INT, DB_NODECENT, NULL);
174 caltinay 3691
175 caltinay 3698 // write element ids
176     dims = getNumElementsPerDim();
177     DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
178     &dims[0], 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
179    
180     // rank 0 writes multimesh and multivar
181 caltinay 3691 if (m_mpiInfo->rank == 0) {
182     vector<string> tempstrings;
183     vector<char*> names;
184     for (dim_t i=0; i<m_mpiInfo->size; i++) {
185     stringstream path;
186     path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
187     tempstrings.push_back(path.str());
188     names.push_back((char*)tempstrings.back().c_str());
189     }
190     vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
191     DBSetDir(dbfile, "/");
192     DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
193     &types[0], NULL);
194     tempstrings.clear();
195     names.clear();
196     for (dim_t i=0; i<m_mpiInfo->size; i++) {
197     stringstream path;
198     path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
199     tempstrings.push_back(path.str());
200     names.push_back((char*)tempstrings.back().c_str());
201     }
202     types.assign(m_mpiInfo->size, DB_QUADVAR);
203     DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
204     &types[0], NULL);
205 caltinay 3698 tempstrings.clear();
206     names.clear();
207     for (dim_t i=0; i<m_mpiInfo->size; i++) {
208     stringstream path;
209     path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
210     tempstrings.push_back(path.str());
211     names.push_back((char*)tempstrings.back().c_str());
212     }
213     DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
214     &types[0], NULL);
215 caltinay 3691 }
216    
217     if (m_mpiInfo->size > 1) {
218     #ifdef ESYS_MPI
219     PMPIO_HandOffBaton(baton, dbfile);
220     PMPIO_Finish(baton);
221     #endif
222     } else {
223     DBClose(dbfile);
224     }
225    
226     #else // USE_SILO
227     throw RipleyException("dump(): no Silo support");
228     #endif
229     }
230    
231     const int* Brick::borrowSampleReferenceIDs(int fsType) const
232     {
233 caltinay 3697 switch (fsType) {
234     case Nodes:
235     return &m_nodeId[0];
236     case Elements:
237     return &m_elementId[0];
238     case FaceElements:
239     return &m_faceId[0];
240     default:
241     break;
242     }
243 caltinay 3691
244 caltinay 3697 stringstream msg;
245     msg << "borrowSampleReferenceIDs() not implemented for function space type "
246     << fsType;
247     throw RipleyException(msg.str());
248 caltinay 3691 }
249    
250     bool Brick::ownSample(int fsCode, index_t id) const
251     {
252     #ifdef ESYS_MPI
253     if (fsCode == Nodes) {
254 caltinay 3698 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
255     const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
256 caltinay 3691 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
257     } else
258     throw RipleyException("ownSample() only implemented for Nodes");
259     #else
260     return true;
261     #endif
262     }
263    
264 caltinay 3703 void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const
265     {
266     escript::Data& in = *const_cast<escript::Data*>(&cIn);
267     const dim_t numComp = in.getDataPointSize();
268     const double h0 = m_l0/m_gNE0;
269     const double h1 = m_l1/m_gNE1;
270     const double h2 = m_l1/m_gNE2;
271     if (out.getFunctionSpace().getTypeCode() == Elements) {
272     /* GENERATOR SNIP_GRAD_ELEMENTS TOP */
273     const double tmp0_22 = -0.044658198738520451079/h1;
274     const double tmp0_16 = 0.16666666666666666667/h0;
275     const double tmp0_33 = -0.62200846792814621559/h1;
276     const double tmp0_0 = -0.62200846792814621559/h0;
277     const double tmp0_21 = -0.16666666666666666667/h1;
278     const double tmp0_17 = 0.62200846792814621559/h0;
279     const double tmp0_52 = -0.044658198738520451079/h2;
280     const double tmp0_1 = -0.16666666666666666667/h0;
281     const double tmp0_20 = -0.62200846792814621559/h1;
282     const double tmp0_14 = -0.044658198738520451079/h0;
283     const double tmp0_53 = -0.62200846792814621559/h2;
284     const double tmp0_49 = 0.16666666666666666667/h2;
285     const double tmp0_2 = 0.16666666666666666667/h0;
286     const double tmp0_27 = -0.044658198738520451079/h1;
287     const double tmp0_15 = -0.16666666666666666667/h0;
288     const double tmp0_50 = -0.16666666666666666667/h2;
289     const double tmp0_48 = 0.62200846792814621559/h2;
290     const double tmp0_3 = 0.044658198738520451079/h0;
291     const double tmp0_26 = -0.16666666666666666667/h1;
292     const double tmp0_12 = -0.62200846792814621559/h0;
293     const double tmp0_51 = 0.044658198738520451079/h2;
294     const double tmp0_25 = 0.62200846792814621559/h1;
295     const double tmp0_13 = 0.16666666666666666667/h0;
296     const double tmp0_56 = 0.16666666666666666667/h2;
297     const double tmp0_24 = 0.16666666666666666667/h1;
298     const double tmp0_10 = 0.62200846792814621559/h0;
299     const double tmp0_57 = 0.62200846792814621559/h2;
300     const double tmp0_11 = -0.16666666666666666667/h0;
301     const double tmp0_54 = -0.044658198738520451079/h2;
302     const double tmp0_38 = 0.16666666666666666667/h1;
303     const double tmp0_34 = -0.044658198738520451079/h1;
304     const double tmp0_42 = 0.16666666666666666667/h2;
305     const double tmp0_35 = -0.16666666666666666667/h1;
306     const double tmp0_36 = -0.62200846792814621559/h1;
307     const double tmp0_41 = 0.62200846792814621559/h2;
308     const double tmp0_8 = 0.044658198738520451079/h0;
309     const double tmp0_37 = 0.62200846792814621559/h1;
310     const double tmp0_29 = 0.16666666666666666667/h1;
311     const double tmp0_40 = -0.62200846792814621559/h2;
312     const double tmp0_9 = 0.16666666666666666667/h0;
313     const double tmp0_30 = 0.62200846792814621559/h1;
314     const double tmp0_28 = -0.16666666666666666667/h1;
315     const double tmp0_43 = 0.044658198738520451079/h2;
316     const double tmp0_32 = 0.16666666666666666667/h1;
317     const double tmp0_31 = 0.044658198738520451079/h1;
318     const double tmp0_39 = 0.044658198738520451079/h1;
319     const double tmp0_58 = -0.62200846792814621559/h2;
320     const double tmp0_55 = 0.044658198738520451079/h2;
321     const double tmp0_18 = -0.62200846792814621559/h0;
322     const double tmp0_45 = -0.16666666666666666667/h2;
323     const double tmp0_59 = -0.16666666666666666667/h2;
324     const double tmp0_4 = -0.044658198738520451079/h0;
325     const double tmp0_19 = 0.044658198738520451079/h0;
326     const double tmp0_44 = -0.044658198738520451079/h2;
327     const double tmp0_5 = 0.62200846792814621559/h0;
328     const double tmp0_47 = 0.16666666666666666667/h2;
329     const double tmp0_6 = -0.16666666666666666667/h0;
330     const double tmp0_23 = 0.044658198738520451079/h1;
331     const double tmp0_46 = -0.16666666666666666667/h2;
332     const double tmp0_7 = -0.044658198738520451079/h0;
333     #pragma omp parallel for
334 caltinay 3707 for (index_t k2=0; k2 < m_NE2; ++k2) {
335     for (index_t k1=0; k1 < m_NE1; ++k1) {
336     for (index_t k0=0; k0 < m_NE0; ++k0) {
337 caltinay 3703 const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
338     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
339     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
340     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
341     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
342     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
343     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
344     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
345     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
346     for (index_t i=0; i < numComp; ++i) {
347     o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
348     o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);
349     o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);
350     o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
351     o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
352     o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
353     o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
354     o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);
355     o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
356     o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
357     o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
358     o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);
359     o[INDEX3(i,0,4,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
360     o[INDEX3(i,1,4,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;
361     o[INDEX3(i,2,4,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);
362     o[INDEX3(i,0,5,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
363     o[INDEX3(i,1,5,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);
364     o[INDEX3(i,2,5,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
365     o[INDEX3(i,0,6,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
366     o[INDEX3(i,1,6,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;
367     o[INDEX3(i,2,6,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
368     o[INDEX3(i,0,7,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
369     o[INDEX3(i,1,7,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);
370     o[INDEX3(i,2,7,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);
371     } /* end of component loop i */
372     } /* end of k0 loop */
373     } /* end of k1 loop */
374     } /* end of k2 loop */
375     /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
376 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
377     /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
378     const double tmp0_0 = -0.25/h0;
379     const double tmp0_1 = 0.25/h0;
380     const double tmp0_2 = -0.25/h1;
381     const double tmp0_3 = 0.25/h1;
382 caltinay 3722 const double tmp0_4 = -0.25/h2;
383     const double tmp0_5 = 0.25/h2;
384 caltinay 3711 #pragma omp parallel for
385     for (index_t k2=0; k2 < m_NE2; ++k2) {
386     for (index_t k1=0; k1 < m_NE1; ++k1) {
387     for (index_t k0=0; k0 < m_NE0; ++k0) {
388     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
389     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
390     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
391     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
392     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
393     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
394     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
395     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
396     double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
397     for (index_t i=0; i < numComp; ++i) {
398     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
399     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i] + f_100[i] + f_101[i]) + tmp0_3*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
400     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i] + f_100[i] + f_110[i]) + tmp0_5*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
401     } /* end of component loop i */
402     } /* end of k0 loop */
403     } /* end of k1 loop */
404     } /* end of k2 loop */
405     /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
406 caltinay 3707 } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
407     /* GENERATOR SNIP_GRAD_FACES TOP */
408 caltinay 3722 #pragma omp parallel
409     {
410     if (m_faceOffset[0] > -1) {
411     const double tmp0_22 = 0.21132486540518711775/h1;
412     const double tmp0_16 = 0.16666666666666666667/h0;
413     const double tmp0_33 = 0.21132486540518711775/h2;
414     const double tmp0_0 = -0.62200846792814621559/h0;
415     const double tmp0_21 = -0.21132486540518711775/h1;
416     const double tmp0_17 = 0.62200846792814621559/h0;
417     const double tmp0_1 = -0.16666666666666666667/h0;
418     const double tmp0_20 = -0.78867513459481288225/h1;
419     const double tmp0_14 = -0.044658198738520451079/h0;
420     const double tmp0_2 = 0.16666666666666666667/h0;
421     const double tmp0_27 = 0.21132486540518711775/h1;
422     const double tmp0_15 = -0.16666666666666666667/h0;
423     const double tmp0_3 = 0.044658198738520451079/h0;
424     const double tmp0_26 = 0.78867513459481288225/h1;
425     const double tmp0_12 = -0.62200846792814621559/h0;
426     const double tmp0_25 = -0.78867513459481288225/h1;
427     const double tmp0_13 = 0.16666666666666666667/h0;
428     const double tmp0_24 = -0.21132486540518711775/h1;
429     const double tmp0_10 = 0.62200846792814621559/h0;
430     const double tmp0_11 = -0.16666666666666666667/h0;
431     const double tmp0_34 = 0.78867513459481288225/h2;
432     const double tmp0_35 = -0.78867513459481288225/h2;
433     const double tmp0_8 = 0.044658198738520451079/h0;
434     const double tmp0_29 = 0.78867513459481288225/h2;
435     const double tmp0_9 = 0.16666666666666666667/h0;
436     const double tmp0_30 = 0.21132486540518711775/h2;
437     const double tmp0_28 = -0.78867513459481288225/h2;
438     const double tmp0_32 = -0.21132486540518711775/h2;
439     const double tmp0_31 = -0.21132486540518711775/h2;
440     const double tmp0_18 = -0.62200846792814621559/h0;
441     const double tmp0_4 = -0.044658198738520451079/h0;
442     const double tmp0_19 = 0.044658198738520451079/h0;
443     const double tmp0_5 = 0.62200846792814621559/h0;
444     const double tmp0_6 = -0.16666666666666666667/h0;
445     const double tmp0_23 = 0.78867513459481288225/h1;
446     const double tmp0_7 = -0.044658198738520451079/h0;
447     #pragma omp for nowait
448     for (index_t k2=0; k2 < m_NE2; ++k2) {
449     for (index_t k1=0; k1 < m_NE1; ++k1) {
450     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
451     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
452     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
453     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
454     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
455     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
456     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
457     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
458     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
459     for (index_t i=0; i < numComp; ++i) {
460     o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
461     o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_20 + f_001[i]*tmp0_21 + f_010[i]*tmp0_23 + f_011[i]*tmp0_22;
462     o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_010[i]*tmp0_31 + f_011[i]*tmp0_30;
463     o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
464     o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_20 + f_001[i]*tmp0_21 + f_010[i]*tmp0_23 + f_011[i]*tmp0_22;
465     o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_010[i]*tmp0_35 + f_011[i]*tmp0_34;
466     o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
467     o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_24 + f_001[i]*tmp0_25 + f_010[i]*tmp0_27 + f_011[i]*tmp0_26;
468     o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_010[i]*tmp0_31 + f_011[i]*tmp0_30;
469     o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
470     o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_24 + f_001[i]*tmp0_25 + f_010[i]*tmp0_27 + f_011[i]*tmp0_26;
471     o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_010[i]*tmp0_35 + f_011[i]*tmp0_34;
472     } /* end of component loop i */
473     } /* end of k1 loop */
474     } /* end of k2 loop */
475     } /* end of face 0 */
476     if (m_faceOffset[1] > -1) {
477     const double tmp0_22 = 0.78867513459481288225/h1;
478     const double tmp0_16 = 0.16666666666666666667/h0;
479     const double tmp0_33 = 0.78867513459481288225/h2;
480     const double tmp0_0 = -0.62200846792814621559/h0;
481     const double tmp0_21 = 0.21132486540518711775/h1;
482     const double tmp0_17 = 0.62200846792814621559/h0;
483     const double tmp0_1 = -0.16666666666666666667/h0;
484     const double tmp0_20 = -0.21132486540518711775/h1;
485     const double tmp0_14 = -0.044658198738520451079/h0;
486     const double tmp0_2 = 0.16666666666666666667/h0;
487     const double tmp0_27 = -0.21132486540518711775/h1;
488     const double tmp0_15 = -0.16666666666666666667/h0;
489     const double tmp0_3 = 0.044658198738520451079/h0;
490     const double tmp0_26 = 0.21132486540518711775/h1;
491     const double tmp0_12 = -0.62200846792814621559/h0;
492     const double tmp0_25 = 0.78867513459481288225/h1;
493     const double tmp0_13 = 0.16666666666666666667/h0;
494     const double tmp0_24 = -0.78867513459481288225/h1;
495     const double tmp0_10 = 0.62200846792814621559/h0;
496     const double tmp0_11 = -0.16666666666666666667/h0;
497     const double tmp0_34 = -0.78867513459481288225/h2;
498     const double tmp0_35 = -0.21132486540518711775/h2;
499     const double tmp0_8 = 0.044658198738520451079/h0;
500     const double tmp0_29 = 0.21132486540518711775/h2;
501     const double tmp0_9 = 0.16666666666666666667/h0;
502     const double tmp0_30 = -0.21132486540518711775/h2;
503     const double tmp0_28 = 0.78867513459481288225/h2;
504     const double tmp0_32 = 0.21132486540518711775/h2;
505     const double tmp0_31 = -0.78867513459481288225/h2;
506     const double tmp0_18 = -0.62200846792814621559/h0;
507     const double tmp0_4 = -0.044658198738520451079/h0;
508     const double tmp0_19 = 0.044658198738520451079/h0;
509     const double tmp0_5 = 0.62200846792814621559/h0;
510     const double tmp0_6 = -0.16666666666666666667/h0;
511     const double tmp0_23 = -0.78867513459481288225/h1;
512     const double tmp0_7 = -0.044658198738520451079/h0;
513     #pragma omp for nowait
514     for (index_t k2=0; k2 < m_NE2; ++k2) {
515     for (index_t k1=0; k1 < m_NE1; ++k1) {
516     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
517     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
518     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
519     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
520     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
521     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
522     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
523     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
524     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
525     for (index_t i=0; i < numComp; ++i) {
526     o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
527     o[INDEX3(i,1,0,numComp,3)] = f_100[i]*tmp0_23 + f_101[i]*tmp0_20 + f_110[i]*tmp0_22 + f_111[i]*tmp0_21;
528     o[INDEX3(i,2,0,numComp,3)] = f_100[i]*tmp0_31 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
529     o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
530     o[INDEX3(i,1,1,numComp,3)] = f_100[i]*tmp0_23 + f_101[i]*tmp0_20 + f_110[i]*tmp0_22 + f_111[i]*tmp0_21;
531     o[INDEX3(i,2,1,numComp,3)] = f_100[i]*tmp0_35 + f_101[i]*tmp0_32 + f_110[i]*tmp0_34 + f_111[i]*tmp0_33;
532     o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
533     o[INDEX3(i,1,2,numComp,3)] = f_100[i]*tmp0_27 + f_101[i]*tmp0_24 + f_110[i]*tmp0_26 + f_111[i]*tmp0_25;
534     o[INDEX3(i,2,2,numComp,3)] = f_100[i]*tmp0_31 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
535     o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
536     o[INDEX3(i,1,3,numComp,3)] = f_100[i]*tmp0_27 + f_101[i]*tmp0_24 + f_110[i]*tmp0_26 + f_111[i]*tmp0_25;
537     o[INDEX3(i,2,3,numComp,3)] = f_100[i]*tmp0_35 + f_101[i]*tmp0_32 + f_110[i]*tmp0_34 + f_111[i]*tmp0_33;
538     } /* end of component loop i */
539     } /* end of k1 loop */
540     } /* end of k2 loop */
541     } /* end of face 1 */
542     if (m_faceOffset[2] > -1) {
543     const double tmp0_22 = -0.044658198738520451079/h1;
544     const double tmp0_16 = -0.16666666666666666667/h1;
545     const double tmp0_33 = 0.21132486540518711775/h2;
546     const double tmp0_0 = -0.78867513459481288225/h0;
547     const double tmp0_21 = 0.16666666666666666667/h1;
548     const double tmp0_17 = -0.62200846792814621559/h1;
549     const double tmp0_1 = -0.21132486540518711775/h0;
550     const double tmp0_20 = 0.044658198738520451079/h1;
551     const double tmp0_14 = -0.16666666666666666667/h1;
552     const double tmp0_2 = 0.21132486540518711775/h0;
553     const double tmp0_27 = 0.044658198738520451079/h1;
554     const double tmp0_15 = -0.044658198738520451079/h1;
555     const double tmp0_3 = 0.78867513459481288225/h0;
556     const double tmp0_26 = 0.16666666666666666667/h1;
557     const double tmp0_12 = 0.16666666666666666667/h1;
558     const double tmp0_25 = 0.62200846792814621559/h1;
559     const double tmp0_13 = 0.62200846792814621559/h1;
560     const double tmp0_24 = -0.62200846792814621559/h1;
561     const double tmp0_10 = -0.044658198738520451079/h1;
562     const double tmp0_11 = 0.044658198738520451079/h1;
563     const double tmp0_34 = 0.78867513459481288225/h2;
564     const double tmp0_35 = -0.78867513459481288225/h2;
565     const double tmp0_8 = -0.62200846792814621559/h1;
566     const double tmp0_29 = 0.78867513459481288225/h2;
567     const double tmp0_9 = -0.16666666666666666667/h1;
568     const double tmp0_30 = 0.21132486540518711775/h2;
569     const double tmp0_28 = -0.78867513459481288225/h2;
570     const double tmp0_32 = -0.21132486540518711775/h2;
571     const double tmp0_31 = -0.21132486540518711775/h2;
572     const double tmp0_18 = 0.16666666666666666667/h1;
573     const double tmp0_4 = -0.21132486540518711775/h0;
574     const double tmp0_19 = 0.62200846792814621559/h1;
575     const double tmp0_5 = -0.78867513459481288225/h0;
576     const double tmp0_6 = 0.78867513459481288225/h0;
577     const double tmp0_23 = -0.16666666666666666667/h1;
578     const double tmp0_7 = 0.21132486540518711775/h0;
579     #pragma omp for nowait
580     for (index_t k2=0; k2 < m_NE2; ++k2) {
581     for (index_t k0=0; k0 < m_NE0; ++k0) {
582     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
583     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
584     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
585     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
586     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
587     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
588     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
589     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
590     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
591     for (index_t i=0; i < numComp; ++i) {
592     o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_001[i]*tmp0_1 + f_100[i]*tmp0_3 + f_101[i]*tmp0_2;
593     o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_13 + f_101[i]*tmp0_10 + f_111[i]*tmp0_11 + tmp0_12*(f_011[i] + f_110[i]) + tmp0_9*(f_001[i] + f_100[i]);
594     o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_100[i]*tmp0_31 + f_101[i]*tmp0_30;
595     o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_001[i]*tmp0_1 + f_100[i]*tmp0_3 + f_101[i]*tmp0_2;
596     o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_14 + f_001[i]*tmp0_15 + f_010[i]*tmp0_21 + f_011[i]*tmp0_20 + f_100[i]*tmp0_17 + f_101[i]*tmp0_16 + f_110[i]*tmp0_19 + f_111[i]*tmp0_18;
597     o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_100[i]*tmp0_35 + f_101[i]*tmp0_34;
598     o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_4 + f_001[i]*tmp0_5 + f_100[i]*tmp0_7 + f_101[i]*tmp0_6;
599     o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_14 + f_001[i]*tmp0_17 + f_010[i]*tmp0_21 + f_011[i]*tmp0_19 + f_100[i]*tmp0_15 + f_101[i]*tmp0_16 + f_110[i]*tmp0_20 + f_111[i]*tmp0_18;
600     o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_28 + f_001[i]*tmp0_29 + f_100[i]*tmp0_31 + f_101[i]*tmp0_30;
601     o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_4 + f_001[i]*tmp0_5 + f_100[i]*tmp0_7 + f_101[i]*tmp0_6;
602     o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_22 + f_010[i]*tmp0_27 + f_101[i]*tmp0_24 + f_111[i]*tmp0_25 + tmp0_23*(f_001[i] + f_100[i]) + tmp0_26*(f_011[i] + f_110[i]);
603     o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_32 + f_001[i]*tmp0_33 + f_100[i]*tmp0_35 + f_101[i]*tmp0_34;
604     } /* end of component loop i */
605     } /* end of k0 loop */
606     } /* end of k2 loop */
607     } /* end of face 2 */
608     if (m_faceOffset[3] > -1) {
609     const double tmp0_22 = 0.16666666666666666667/h1;
610     const double tmp0_16 = 0.16666666666666666667/h1;
611     const double tmp0_33 = -0.78867513459481288225/h2;
612     const double tmp0_0 = -0.21132486540518711775/h0;
613     const double tmp0_21 = -0.62200846792814621559/h1;
614     const double tmp0_17 = 0.16666666666666666667/h1;
615     const double tmp0_1 = 0.78867513459481288225/h0;
616     const double tmp0_20 = -0.16666666666666666667/h1;
617     const double tmp0_14 = 0.044658198738520451079/h1;
618     const double tmp0_2 = -0.78867513459481288225/h0;
619     const double tmp0_27 = -0.62200846792814621559/h1;
620     const double tmp0_15 = 0.62200846792814621559/h1;
621     const double tmp0_3 = 0.21132486540518711775/h0;
622     const double tmp0_26 = -0.16666666666666666667/h1;
623     const double tmp0_12 = -0.16666666666666666667/h1;
624     const double tmp0_25 = -0.044658198738520451079/h1;
625     const double tmp0_13 = -0.044658198738520451079/h1;
626     const double tmp0_24 = 0.62200846792814621559/h1;
627     const double tmp0_10 = 0.044658198738520451079/h1;
628     const double tmp0_11 = -0.62200846792814621559/h1;
629     const double tmp0_34 = -0.21132486540518711775/h2;
630     const double tmp0_35 = 0.78867513459481288225/h2;
631     const double tmp0_8 = 0.16666666666666666667/h1;
632     const double tmp0_29 = -0.21132486540518711775/h2;
633     const double tmp0_9 = 0.62200846792814621559/h1;
634     const double tmp0_30 = -0.78867513459481288225/h2;
635     const double tmp0_28 = 0.78867513459481288225/h2;
636     const double tmp0_32 = 0.21132486540518711775/h2;
637     const double tmp0_31 = 0.21132486540518711775/h2;
638     const double tmp0_18 = -0.16666666666666666667/h1;
639     const double tmp0_4 = -0.78867513459481288225/h0;
640     const double tmp0_19 = -0.044658198738520451079/h1;
641     const double tmp0_5 = 0.21132486540518711775/h0;
642     const double tmp0_6 = -0.21132486540518711775/h0;
643     const double tmp0_23 = 0.044658198738520451079/h1;
644     const double tmp0_7 = 0.78867513459481288225/h0;
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_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
654     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,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     o[INDEX3(i,0,0,numComp,3)] = f_010[i]*tmp0_2 + f_011[i]*tmp0_0 + f_110[i]*tmp0_1 + f_111[i]*tmp0_3;
659     o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_11 + f_010[i]*tmp0_9 + f_101[i]*tmp0_13 + f_111[i]*tmp0_10 + tmp0_12*(f_001[i] + f_100[i]) + tmp0_8*(f_011[i] + f_110[i]);
660     o[INDEX3(i,2,0,numComp,3)] = f_010[i]*tmp0_30 + f_011[i]*tmp0_28 + f_110[i]*tmp0_29 + f_111[i]*tmp0_31;
661     o[INDEX3(i,0,1,numComp,3)] = f_010[i]*tmp0_2 + f_011[i]*tmp0_0 + f_110[i]*tmp0_1 + f_111[i]*tmp0_3;
662     o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_18 + f_001[i]*tmp0_19 + f_010[i]*tmp0_16 + f_011[i]*tmp0_14 + f_100[i]*tmp0_21 + f_101[i]*tmp0_20 + f_110[i]*tmp0_15 + f_111[i]*tmp0_17;
663     o[INDEX3(i,2,1,numComp,3)] = f_010[i]*tmp0_34 + f_011[i]*tmp0_32 + f_110[i]*tmp0_33 + f_111[i]*tmp0_35;
664     o[INDEX3(i,0,2,numComp,3)] = f_010[i]*tmp0_6 + f_011[i]*tmp0_4 + f_110[i]*tmp0_5 + f_111[i]*tmp0_7;
665     o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_18 + f_001[i]*tmp0_21 + f_010[i]*tmp0_16 + f_011[i]*tmp0_15 + f_100[i]*tmp0_19 + f_101[i]*tmp0_20 + f_110[i]*tmp0_14 + f_111[i]*tmp0_17;
666     o[INDEX3(i,2,2,numComp,3)] = f_010[i]*tmp0_30 + f_011[i]*tmp0_28 + f_110[i]*tmp0_29 + f_111[i]*tmp0_31;
667     o[INDEX3(i,0,3,numComp,3)] = f_010[i]*tmp0_6 + f_011[i]*tmp0_4 + f_110[i]*tmp0_5 + f_111[i]*tmp0_7;
668     o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_25 + f_010[i]*tmp0_23 + f_101[i]*tmp0_27 + f_111[i]*tmp0_24 + tmp0_22*(f_011[i] + f_110[i]) + tmp0_26*(f_001[i] + f_100[i]);
669     o[INDEX3(i,2,3,numComp,3)] = f_010[i]*tmp0_34 + f_011[i]*tmp0_32 + f_110[i]*tmp0_33 + f_111[i]*tmp0_35;
670     } /* end of component loop i */
671     } /* end of k0 loop */
672     } /* end of k2 loop */
673     } /* end of face 3 */
674     if (m_faceOffset[4] > -1) {
675     const double tmp0_22 = -0.16666666666666666667/h2;
676     const double tmp0_16 = -0.62200846792814621559/h2;
677     const double tmp0_33 = 0.044658198738520451079/h2;
678     const double tmp0_0 = -0.78867513459481288225/h0;
679     const double tmp0_21 = 0.044658198738520451079/h2;
680     const double tmp0_17 = -0.16666666666666666667/h2;
681     const double tmp0_1 = 0.78867513459481288225/h0;
682     const double tmp0_20 = 0.16666666666666666667/h2;
683     const double tmp0_14 = 0.78867513459481288225/h1;
684     const double tmp0_2 = 0.21132486540518711775/h0;
685     const double tmp0_27 = 0.62200846792814621559/h2;
686     const double tmp0_15 = 0.21132486540518711775/h1;
687     const double tmp0_3 = -0.21132486540518711775/h0;
688     const double tmp0_26 = 0.16666666666666666667/h2;
689     const double tmp0_12 = -0.21132486540518711775/h1;
690     const double tmp0_25 = -0.044658198738520451079/h2;
691     const double tmp0_13 = -0.78867513459481288225/h1;
692     const double tmp0_24 = -0.16666666666666666667/h2;
693     const double tmp0_10 = 0.21132486540518711775/h1;
694     const double tmp0_11 = 0.78867513459481288225/h1;
695     const double tmp0_34 = 0.16666666666666666667/h2;
696     const double tmp0_35 = 0.62200846792814621559/h2;
697     const double tmp0_8 = -0.78867513459481288225/h1;
698     const double tmp0_29 = 0.16666666666666666667/h2;
699     const double tmp0_9 = -0.21132486540518711775/h1;
700     const double tmp0_30 = -0.044658198738520451079/h2;
701     const double tmp0_28 = 0.044658198738520451079/h2;
702     const double tmp0_32 = -0.62200846792814621559/h2;
703     const double tmp0_31 = -0.16666666666666666667/h2;
704     const double tmp0_18 = -0.044658198738520451079/h2;
705     const double tmp0_4 = -0.21132486540518711775/h0;
706     const double tmp0_19 = 0.62200846792814621559/h2;
707     const double tmp0_5 = 0.21132486540518711775/h0;
708     const double tmp0_6 = 0.78867513459481288225/h0;
709     const double tmp0_23 = -0.62200846792814621559/h2;
710     const double tmp0_7 = -0.78867513459481288225/h0;
711     #pragma omp for nowait
712 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
713 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
714     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
715     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
716     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
717     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
718     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
719     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
720     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
721     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
722     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
723     for (index_t i=0; i < numComp; ++i) {
724     o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_010[i]*tmp0_3 + f_100[i]*tmp0_1 + f_110[i]*tmp0_2;
725     o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_11 + f_100[i]*tmp0_9 + f_110[i]*tmp0_10;
726     o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_16 + f_001[i]*tmp0_19 + f_110[i]*tmp0_18 + f_111[i]*tmp0_21 + tmp0_17*(f_010[i] + f_100[i]) + tmp0_20*(f_011[i] + f_101[i]);
727     o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_010[i]*tmp0_3 + f_100[i]*tmp0_1 + f_110[i]*tmp0_2;
728     o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_12 + f_010[i]*tmp0_15 + f_100[i]*tmp0_13 + f_110[i]*tmp0_14;
729     o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_22 + f_001[i]*tmp0_26 + f_010[i]*tmp0_25 + f_011[i]*tmp0_28 + f_100[i]*tmp0_23 + f_101[i]*tmp0_27 + f_110[i]*tmp0_24 + f_111[i]*tmp0_29;
730     o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_4 + f_010[i]*tmp0_7 + f_100[i]*tmp0_5 + f_110[i]*tmp0_6;
731     o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_8 + f_010[i]*tmp0_11 + f_100[i]*tmp0_9 + f_110[i]*tmp0_10;
732     o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_22 + f_001[i]*tmp0_26 + f_010[i]*tmp0_23 + f_011[i]*tmp0_27 + f_100[i]*tmp0_25 + f_101[i]*tmp0_28 + f_110[i]*tmp0_24 + f_111[i]*tmp0_29;
733     o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_4 + f_010[i]*tmp0_7 + f_100[i]*tmp0_5 + f_110[i]*tmp0_6;
734     o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_12 + f_010[i]*tmp0_15 + f_100[i]*tmp0_13 + f_110[i]*tmp0_14;
735     o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_30 + f_001[i]*tmp0_33 + f_110[i]*tmp0_32 + f_111[i]*tmp0_35 + tmp0_31*(f_010[i] + f_100[i]) + tmp0_34*(f_011[i] + f_101[i]);
736     } /* end of component loop i */
737     } /* end of k0 loop */
738 caltinay 3707 } /* end of k1 loop */
739 caltinay 3722 } /* end of face 4 */
740     if (m_faceOffset[5] > -1) {
741     const double tmp0_22 = 0.16666666666666666667/h2;
742     const double tmp0_16 = 0.62200846792814621559/h2;
743     const double tmp0_33 = -0.044658198738520451079/h2;
744     const double tmp0_0 = -0.78867513459481288225/h0;
745     const double tmp0_21 = -0.16666666666666666667/h2;
746     const double tmp0_17 = 0.16666666666666666667/h2;
747     const double tmp0_1 = 0.78867513459481288225/h0;
748     const double tmp0_20 = -0.044658198738520451079/h2;
749     const double tmp0_14 = 0.21132486540518711775/h1;
750     const double tmp0_2 = -0.21132486540518711775/h0;
751     const double tmp0_27 = -0.16666666666666666667/h2;
752     const double tmp0_15 = 0.78867513459481288225/h1;
753     const double tmp0_3 = 0.21132486540518711775/h0;
754     const double tmp0_26 = -0.16666666666666666667/h2;
755     const double tmp0_12 = -0.21132486540518711775/h1;
756     const double tmp0_25 = 0.16666666666666666667/h2;
757     const double tmp0_13 = -0.78867513459481288225/h1;
758     const double tmp0_24 = 0.044658198738520451079/h2;
759     const double tmp0_10 = 0.78867513459481288225/h1;
760     const double tmp0_11 = 0.21132486540518711775/h1;
761     const double tmp0_34 = -0.62200846792814621559/h2;
762     const double tmp0_35 = -0.16666666666666666667/h2;
763     const double tmp0_8 = -0.78867513459481288225/h1;
764     const double tmp0_29 = -0.62200846792814621559/h2;
765     const double tmp0_9 = -0.21132486540518711775/h1;
766     const double tmp0_30 = 0.044658198738520451079/h2;
767     const double tmp0_28 = -0.044658198738520451079/h2;
768     const double tmp0_32 = 0.62200846792814621559/h2;
769     const double tmp0_31 = 0.16666666666666666667/h2;
770     const double tmp0_18 = 0.044658198738520451079/h2;
771     const double tmp0_4 = -0.21132486540518711775/h0;
772     const double tmp0_19 = -0.62200846792814621559/h2;
773     const double tmp0_5 = 0.21132486540518711775/h0;
774     const double tmp0_6 = -0.78867513459481288225/h0;
775     const double tmp0_23 = 0.62200846792814621559/h2;
776     const double tmp0_7 = 0.78867513459481288225/h0;
777     #pragma omp for nowait
778 caltinay 3707 for (index_t k1=0; k1 < m_NE1; ++k1) {
779 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
780     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
781     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
782     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
783     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
784     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
785     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
786     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
787     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
788     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
789     for (index_t i=0; i < numComp; ++i) {
790     o[INDEX3(i,0,0,numComp,3)] = f_001[i]*tmp0_0 + f_011[i]*tmp0_2 + f_101[i]*tmp0_1 + f_111[i]*tmp0_3;
791     o[INDEX3(i,1,0,numComp,3)] = f_001[i]*tmp0_8 + f_011[i]*tmp0_10 + f_101[i]*tmp0_9 + f_111[i]*tmp0_11;
792     o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_19 + f_001[i]*tmp0_16 + f_110[i]*tmp0_20 + f_111[i]*tmp0_18 + tmp0_17*(f_011[i] + f_101[i]) + tmp0_21*(f_010[i] + f_100[i]);
793     o[INDEX3(i,0,1,numComp,3)] = f_001[i]*tmp0_0 + f_011[i]*tmp0_2 + f_101[i]*tmp0_1 + f_111[i]*tmp0_3;
794     o[INDEX3(i,1,1,numComp,3)] = f_001[i]*tmp0_12 + f_011[i]*tmp0_14 + f_101[i]*tmp0_13 + f_111[i]*tmp0_15;
795     o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_22 + f_010[i]*tmp0_28 + f_011[i]*tmp0_24 + f_100[i]*tmp0_29 + f_101[i]*tmp0_23 + f_110[i]*tmp0_27 + f_111[i]*tmp0_25;
796     o[INDEX3(i,0,2,numComp,3)] = f_001[i]*tmp0_4 + f_011[i]*tmp0_6 + f_101[i]*tmp0_5 + f_111[i]*tmp0_7;
797     o[INDEX3(i,1,2,numComp,3)] = f_001[i]*tmp0_8 + f_011[i]*tmp0_10 + f_101[i]*tmp0_9 + f_111[i]*tmp0_11;
798     o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_22 + f_010[i]*tmp0_29 + f_011[i]*tmp0_23 + f_100[i]*tmp0_28 + f_101[i]*tmp0_24 + f_110[i]*tmp0_27 + f_111[i]*tmp0_25;
799     o[INDEX3(i,0,3,numComp,3)] = f_001[i]*tmp0_4 + f_011[i]*tmp0_6 + f_101[i]*tmp0_5 + f_111[i]*tmp0_7;
800     o[INDEX3(i,1,3,numComp,3)] = f_001[i]*tmp0_12 + f_011[i]*tmp0_14 + f_101[i]*tmp0_13 + f_111[i]*tmp0_15;
801     o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_33 + f_001[i]*tmp0_30 + f_110[i]*tmp0_34 + f_111[i]*tmp0_32 + tmp0_31*(f_011[i] + f_101[i]) + tmp0_35*(f_010[i] + f_100[i]);
802     } /* end of component loop i */
803     } /* end of k0 loop */
804 caltinay 3707 } /* end of k1 loop */
805 caltinay 3722 } /* end of face 5 */
806     } // end of parallel section
807 caltinay 3707 /* GENERATOR SNIP_GRAD_FACES BOTTOM */
808 caltinay 3711 } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
809     /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
810 caltinay 3722 #pragma omp parallel
811     {
812     if (m_faceOffset[0] > -1) {
813     const double tmp0_0 = -0.25/h0;
814     const double tmp0_4 = -0.5/h2;
815     const double tmp0_1 = 0.25/h0;
816     const double tmp0_5 = 0.5/h2;
817     const double tmp0_2 = -0.5/h1;
818     const double tmp0_3 = 0.5/h1;
819     #pragma omp for nowait
820     for (index_t k2=0; k2 < m_NE2; ++k2) {
821     for (index_t k1=0; k1 < m_NE1; ++k1) {
822     const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
823     const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
824     const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
825     const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
826     const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
827     const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
828     const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
829     const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
830     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
831     for (index_t i=0; i < numComp; ++i) {
832     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
833     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i]) + tmp0_3*(f_010[i] + f_011[i]);
834     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i]) + tmp0_5*(f_001[i] + f_011[i]);
835     } /* end of component loop i */
836     } /* end of k1 loop */
837     } /* end of k2 loop */
838     } /* end of face 0 */
839     if (m_faceOffset[1] > -1) {
840     const double tmp0_0 = -0.25/h0;
841     const double tmp0_4 = 0.5/h2;
842     const double tmp0_1 = 0.25/h0;
843     const double tmp0_5 = -0.5/h2;
844     const double tmp0_2 = -0.5/h1;
845     const double tmp0_3 = 0.5/h1;
846     #pragma omp for nowait
847     for (index_t k2=0; k2 < m_NE2; ++k2) {
848     for (index_t k1=0; k1 < m_NE1; ++k1) {
849     const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
850     const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
851     const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
852     const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
853     const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
854     const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
855     const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
856     const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
857     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
858     for (index_t i=0; i < numComp; ++i) {
859     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]) + tmp0_1*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
860     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_100[i] + f_101[i]) + tmp0_3*(f_110[i] + f_111[i]);
861     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_101[i] + f_111[i]) + tmp0_5*(f_100[i] + f_110[i]);
862     } /* end of component loop i */
863     } /* end of k1 loop */
864     } /* end of k2 loop */
865     } /* end of face 1 */
866     if (m_faceOffset[2] > -1) {
867     const double tmp0_0 = -0.5/h0;
868     const double tmp0_4 = -0.5/h2;
869     const double tmp0_1 = 0.5/h0;
870     const double tmp0_5 = 0.5/h2;
871     const double tmp0_2 = -0.25/h1;
872     const double tmp0_3 = 0.25/h1;
873     #pragma omp for nowait
874     for (index_t k2=0; k2 < m_NE2; ++k2) {
875     for (index_t k0=0; k0 < m_NE0; ++k0) {
876     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
877     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
878     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
879     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
880     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
881     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
882     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
883     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
884     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
885     for (index_t i=0; i < numComp; ++i) {
886     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_001[i]) + tmp0_1*(f_100[i] + f_101[i]);
887     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_001[i] + f_100[i] + f_101[i]) + tmp0_3*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
888     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_100[i]) + tmp0_5*(f_001[i] + f_101[i]);
889     } /* end of component loop i */
890     } /* end of k0 loop */
891     } /* end of k2 loop */
892     } /* end of face 2 */
893     if (m_faceOffset[3] > -1) {
894     const double tmp0_0 = -0.5/h0;
895     const double tmp0_4 = 0.5/h2;
896     const double tmp0_1 = 0.5/h0;
897     const double tmp0_5 = -0.5/h2;
898     const double tmp0_2 = 0.25/h1;
899     const double tmp0_3 = -0.25/h1;
900     #pragma omp for nowait
901     for (index_t k2=0; k2 < m_NE2; ++k2) {
902     for (index_t k0=0; k0 < m_NE0; ++k0) {
903     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
904     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
905     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
906     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
907     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
908     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
909     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
910     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
911     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
912     for (index_t i=0; i < numComp; ++i) {
913     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_010[i] + f_011[i]) + tmp0_1*(f_110[i] + f_111[i]);
914     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_010[i] + f_011[i] + f_110[i] + f_111[i]) + tmp0_3*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
915     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_011[i] + f_111[i]) + tmp0_5*(f_010[i] + f_110[i]);
916     } /* end of component loop i */
917     } /* end of k0 loop */
918     } /* end of k2 loop */
919     } /* end of face 3 */
920     if (m_faceOffset[4] > -1) {
921     const double tmp0_0 = -0.5/h0;
922     const double tmp0_4 = -0.25/h2;
923     const double tmp0_1 = 0.5/h0;
924     const double tmp0_5 = 0.25/h2;
925     const double tmp0_2 = -0.5/h1;
926     const double tmp0_3 = 0.5/h1;
927     #pragma omp for nowait
928 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
929 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
930     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
931     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
932     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
933     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
934     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
935     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
936     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
937     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
938     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
939     for (index_t i=0; i < numComp; ++i) {
940     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_000[i] + f_010[i]) + tmp0_1*(f_100[i] + f_110[i]);
941     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_000[i] + f_100[i]) + tmp0_3*(f_010[i] + f_110[i]);
942     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_000[i] + f_010[i] + f_100[i] + f_110[i]) + tmp0_5*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
943     } /* end of component loop i */
944     } /* end of k0 loop */
945 caltinay 3711 } /* end of k1 loop */
946 caltinay 3722 } /* end of face 4 */
947     if (m_faceOffset[5] > -1) {
948     const double tmp0_0 = -0.5/h0;
949     const double tmp0_4 = 0.25/h2;
950     const double tmp0_1 = 0.5/h0;
951     const double tmp0_5 = -0.25/h2;
952     const double tmp0_2 = -0.5/h1;
953     const double tmp0_3 = 0.5/h1;
954     #pragma omp for nowait
955 caltinay 3711 for (index_t k1=0; k1 < m_NE1; ++k1) {
956 caltinay 3722 for (index_t k0=0; k0 < m_NE0; ++k0) {
957     const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
958     const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
959     const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
960     const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
961     const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
962     const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
963     const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
964     const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
965     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
966     for (index_t i=0; i < numComp; ++i) {
967     o[INDEX3(i,0,0,numComp,3)] = tmp0_0*(f_001[i] + f_011[i]) + tmp0_1*(f_101[i] + f_111[i]);
968     o[INDEX3(i,1,0,numComp,3)] = tmp0_2*(f_001[i] + f_101[i]) + tmp0_3*(f_011[i] + f_111[i]);
969     o[INDEX3(i,2,0,numComp,3)] = tmp0_4*(f_001[i] + f_011[i] + f_101[i] + f_111[i]) + tmp0_5*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
970     } /* end of component loop i */
971     } /* end of k0 loop */
972 caltinay 3711 } /* end of k1 loop */
973 caltinay 3722 } /* end of face 5 */
974     } // end of parallel section
975 caltinay 3711 /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
976 caltinay 3703 } else {
977 caltinay 3707 stringstream msg;
978     msg << "setToGradient() not implemented for "
979     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
980     throw RipleyException(msg.str());
981 caltinay 3703 }
982     }
983    
984 caltinay 3713 void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
985     {
986     escript::Data& in = *const_cast<escript::Data*>(&arg);
987     const dim_t numComp = in.getDataPointSize();
988     const double h0 = m_l0/m_gNE0;
989     const double h1 = m_l1/m_gNE1;
990     const double h2 = m_l2/m_gNE2;
991     if (arg.getFunctionSpace().getTypeCode() == Elements) {
992     const double w_0 = h0*h1*h2/8.;
993     #pragma omp parallel
994     {
995     vector<double> int_local(numComp, 0);
996 caltinay 3722 #pragma omp for nowait
997 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
998     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
999     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1000     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
1001     for (index_t i=0; i < numComp; ++i) {
1002     const register double f_0 = f[INDEX2(i,0,numComp)];
1003     const register double f_1 = f[INDEX2(i,1,numComp)];
1004     const register double f_2 = f[INDEX2(i,2,numComp)];
1005     const register double f_3 = f[INDEX2(i,3,numComp)];
1006     const register double f_4 = f[INDEX2(i,4,numComp)];
1007     const register double f_5 = f[INDEX2(i,5,numComp)];
1008     const register double f_6 = f[INDEX2(i,6,numComp)];
1009     const register double f_7 = f[INDEX2(i,7,numComp)];
1010     int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
1011     } /* end of component loop i */
1012     } /* end of k0 loop */
1013     } /* end of k1 loop */
1014     } /* end of k2 loop */
1015    
1016     #pragma omp critical
1017     for (index_t i=0; i<numComp; i++)
1018     integrals[i]+=int_local[i];
1019 caltinay 3722 } // end of parallel section
1020 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
1021     const double w_0 = h0*h1*h2;
1022     #pragma omp parallel
1023     {
1024     vector<double> int_local(numComp, 0);
1025 caltinay 3722 #pragma omp for nowait
1026 caltinay 3713 for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1027     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1028     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1029     const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
1030     for (index_t i=0; i < numComp; ++i) {
1031     int_local[i]+=f[i]*w_0;
1032     } /* end of component loop i */
1033     } /* end of k0 loop */
1034     } /* end of k1 loop */
1035     } /* end of k2 loop */
1036    
1037     #pragma omp critical
1038     for (index_t i=0; i<numComp; i++)
1039     integrals[i]+=int_local[i];
1040 caltinay 3722 } // end of parallel section
1041 caltinay 3713 } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
1042     const double w_0 = h1*h2/4.;
1043     const double w_1 = h0*h2/4.;
1044     const double w_2 = h0*h1/4.;
1045     #pragma omp parallel
1046     {
1047     vector<double> int_local(numComp, 0);
1048     if (m_faceOffset[0] > -1) {
1049 caltinay 3722 #pragma omp for nowait
1050 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1051     for (index_t k1=0; k1 < m_NE1; ++k1) {
1052     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1053     for (index_t i=0; i < numComp; ++i) {
1054     const register double f_0 = f[INDEX2(i,0,numComp)];
1055     const register double f_1 = f[INDEX2(i,1,numComp)];
1056     const register double f_2 = f[INDEX2(i,2,numComp)];
1057     const register double f_3 = f[INDEX2(i,3,numComp)];
1058     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1059     } /* end of component loop i */
1060     } /* end of k1 loop */
1061     } /* end of k2 loop */
1062     }
1063    
1064     if (m_faceOffset[1] > -1) {
1065 caltinay 3722 #pragma omp for nowait
1066 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1067     for (index_t k1=0; k1 < m_NE1; ++k1) {
1068     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1069     for (index_t i=0; i < numComp; ++i) {
1070     const register double f_0 = f[INDEX2(i,0,numComp)];
1071     const register double f_1 = f[INDEX2(i,1,numComp)];
1072     const register double f_2 = f[INDEX2(i,2,numComp)];
1073     const register double f_3 = f[INDEX2(i,3,numComp)];
1074     int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1075     } /* end of component loop i */
1076     } /* end of k1 loop */
1077     } /* end of k2 loop */
1078     }
1079    
1080     if (m_faceOffset[2] > -1) {
1081 caltinay 3722 #pragma omp for nowait
1082 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1083     for (index_t k0=0; k0 < m_NE0; ++k0) {
1084     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1085     for (index_t i=0; i < numComp; ++i) {
1086     const register double f_0 = f[INDEX2(i,0,numComp)];
1087     const register double f_1 = f[INDEX2(i,1,numComp)];
1088     const register double f_2 = f[INDEX2(i,2,numComp)];
1089     const register double f_3 = f[INDEX2(i,3,numComp)];
1090     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1091     } /* end of component loop i */
1092     } /* end of k1 loop */
1093     } /* end of k2 loop */
1094     }
1095    
1096     if (m_faceOffset[3] > -1) {
1097 caltinay 3722 #pragma omp for nowait
1098 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1099     for (index_t k0=0; k0 < m_NE0; ++k0) {
1100     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1101     for (index_t i=0; i < numComp; ++i) {
1102     const register double f_0 = f[INDEX2(i,0,numComp)];
1103     const register double f_1 = f[INDEX2(i,1,numComp)];
1104     const register double f_2 = f[INDEX2(i,2,numComp)];
1105     const register double f_3 = f[INDEX2(i,3,numComp)];
1106     int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1107     } /* end of component loop i */
1108     } /* end of k1 loop */
1109     } /* end of k2 loop */
1110     }
1111    
1112     if (m_faceOffset[4] > -1) {
1113 caltinay 3722 #pragma omp for nowait
1114 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
1115     for (index_t k0=0; k0 < m_NE0; ++k0) {
1116     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1117     for (index_t i=0; i < numComp; ++i) {
1118     const register double f_0 = f[INDEX2(i,0,numComp)];
1119     const register double f_1 = f[INDEX2(i,1,numComp)];
1120     const register double f_2 = f[INDEX2(i,2,numComp)];
1121     const register double f_3 = f[INDEX2(i,3,numComp)];
1122     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1123     } /* end of component loop i */
1124     } /* end of k1 loop */
1125     } /* end of k2 loop */
1126     }
1127    
1128     if (m_faceOffset[5] > -1) {
1129 caltinay 3722 #pragma omp for nowait
1130 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
1131     for (index_t k0=0; k0 < m_NE0; ++k0) {
1132     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1133     for (index_t i=0; i < numComp; ++i) {
1134     const register double f_0 = f[INDEX2(i,0,numComp)];
1135     const register double f_1 = f[INDEX2(i,1,numComp)];
1136     const register double f_2 = f[INDEX2(i,2,numComp)];
1137     const register double f_3 = f[INDEX2(i,3,numComp)];
1138     int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1139     } /* end of component loop i */
1140     } /* end of k1 loop */
1141     } /* end of k2 loop */
1142     }
1143    
1144     #pragma omp critical
1145     for (index_t i=0; i<numComp; i++)
1146     integrals[i]+=int_local[i];
1147 caltinay 3722 } // end of parallel section
1148 caltinay 3713
1149     } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1150     const double w_0 = h1*h2;
1151     const double w_1 = h0*h2;
1152     const double w_2 = h0*h1;
1153     #pragma omp parallel
1154     {
1155     vector<double> int_local(numComp, 0);
1156     if (m_faceOffset[0] > -1) {
1157 caltinay 3722 #pragma omp for nowait
1158 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1159     for (index_t k1=0; k1 < m_NE1; ++k1) {
1160     const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1161     for (index_t i=0; i < numComp; ++i) {
1162     int_local[i]+=f[i]*w_0;
1163     } /* end of component loop i */
1164     } /* end of k1 loop */
1165     } /* end of k2 loop */
1166     }
1167    
1168     if (m_faceOffset[1] > -1) {
1169 caltinay 3722 #pragma omp for nowait
1170 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1171     for (index_t k1=0; k1 < m_NE1; ++k1) {
1172     const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1173     for (index_t i=0; i < numComp; ++i) {
1174     int_local[i]+=f[i]*w_0;
1175     } /* end of component loop i */
1176     } /* end of k1 loop */
1177     } /* end of k2 loop */
1178     }
1179    
1180     if (m_faceOffset[2] > -1) {
1181 caltinay 3722 #pragma omp for nowait
1182 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1183     for (index_t k0=0; k0 < m_NE0; ++k0) {
1184     const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1185     for (index_t i=0; i < numComp; ++i) {
1186     int_local[i]+=f[i]*w_1;
1187     } /* end of component loop i */
1188     } /* end of k1 loop */
1189     } /* end of k2 loop */
1190     }
1191    
1192     if (m_faceOffset[3] > -1) {
1193 caltinay 3722 #pragma omp for nowait
1194 caltinay 3713 for (index_t k2=0; k2 < m_NE2; ++k2) {
1195     for (index_t k0=0; k0 < m_NE0; ++k0) {
1196     const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1197     for (index_t i=0; i < numComp; ++i) {
1198     int_local[i]+=f[i]*w_1;
1199     } /* end of component loop i */
1200     } /* end of k1 loop */
1201     } /* end of k2 loop */
1202     }
1203    
1204     if (m_faceOffset[4] > -1) {
1205 caltinay 3722 #pragma omp for nowait
1206 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
1207     for (index_t k0=0; k0 < m_NE0; ++k0) {
1208     const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1209     for (index_t i=0; i < numComp; ++i) {
1210     int_local[i]+=f[i]*w_2;
1211     } /* end of component loop i */
1212     } /* end of k1 loop */
1213     } /* end of k2 loop */
1214     }
1215    
1216     if (m_faceOffset[5] > -1) {
1217 caltinay 3722 #pragma omp for nowait
1218 caltinay 3713 for (index_t k1=0; k1 < m_NE1; ++k1) {
1219     for (index_t k0=0; k0 < m_NE0; ++k0) {
1220     const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1221     for (index_t i=0; i < numComp; ++i) {
1222     int_local[i]+=f[i]*w_2;
1223     } /* end of component loop i */
1224     } /* end of k1 loop */
1225     } /* end of k2 loop */
1226     }
1227    
1228     #pragma omp critical
1229     for (index_t i=0; i<numComp; i++)
1230     integrals[i]+=int_local[i];
1231 caltinay 3722 } // end of parallel section
1232 caltinay 3713
1233     } else {
1234     stringstream msg;
1235     msg << "setToIntegrals() not implemented for "
1236     << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
1237     throw RipleyException(msg.str());
1238     }
1239     }
1240    
1241 caltinay 3722 void Brick::setToNormal(escript::Data& out) const
1242     {
1243     if (out.getFunctionSpace().getTypeCode() == FaceElements) {
1244     #pragma omp parallel
1245     {
1246     if (m_faceOffset[0] > -1) {
1247     #pragma omp for nowait
1248     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1249     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1250     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1251     // set vector at four quadrature points
1252     *o++ = -1.; *o++ = 0.; *o++ = 0.;
1253     *o++ = -1.; *o++ = 0.; *o++ = 0.;
1254     *o++ = -1.; *o++ = 0.; *o++ = 0.;
1255     *o++ = -1.; *o++ = 0.; *o = 0.;
1256     }
1257     }
1258     }
1259    
1260     if (m_faceOffset[1] > -1) {
1261     #pragma omp for nowait
1262     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1263     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1264     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1265     // set vector at four quadrature points
1266     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1267     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1268     *o++ = 1.; *o++ = 0.; *o++ = 0.;
1269     *o++ = 1.; *o++ = 0.; *o = 0.;
1270     }
1271     }
1272     }
1273    
1274     if (m_faceOffset[2] > -1) {
1275     #pragma omp for nowait
1276     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1277     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1278     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1279     // set vector at four quadrature points
1280     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1281     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1282     *o++ = 0.; *o++ = -1.; *o++ = 0.;
1283     *o++ = 0.; *o++ = -1.; *o = 0.;
1284     }
1285     }
1286     }
1287    
1288     if (m_faceOffset[3] > -1) {
1289     #pragma omp for nowait
1290     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1291     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1292     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1293     // set vector at four quadrature points
1294     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1295     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1296     *o++ = 0.; *o++ = 1.; *o++ = 0.;
1297     *o++ = 0.; *o++ = 1.; *o = 0.;
1298     }
1299     }
1300     }
1301    
1302     if (m_faceOffset[4] > -1) {
1303     #pragma omp for nowait
1304     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1305     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1306     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1307     // set vector at four quadrature points
1308     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1309     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1310     *o++ = 0.; *o++ = 0.; *o++ = -1.;
1311     *o++ = 0.; *o++ = 0.; *o = -1.;
1312     }
1313     }
1314     }
1315    
1316     if (m_faceOffset[5] > -1) {
1317     #pragma omp for nowait
1318     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1319     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1320     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1321     // set vector at four quadrature points
1322     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1323     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1324     *o++ = 0.; *o++ = 0.; *o++ = 1.;
1325     *o++ = 0.; *o++ = 0.; *o = 1.;
1326     }
1327     }
1328     }
1329     } // end of parallel section
1330     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1331     #pragma omp parallel
1332     {
1333     if (m_faceOffset[0] > -1) {
1334     #pragma omp for nowait
1335     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1336     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1337     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1338     *o++ = -1.;
1339     *o++ = 0.;
1340     *o = 0.;
1341     }
1342     }
1343     }
1344    
1345     if (m_faceOffset[1] > -1) {
1346     #pragma omp for nowait
1347     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1348     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1349     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1350     *o++ = 1.;
1351     *o++ = 0.;
1352     *o = 0.;
1353     }
1354     }
1355     }
1356    
1357     if (m_faceOffset[2] > -1) {
1358     #pragma omp for nowait
1359     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1360     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1361     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1362     *o++ = 0.;
1363     *o++ = -1.;
1364     *o = 0.;
1365     }
1366     }
1367     }
1368    
1369     if (m_faceOffset[3] > -1) {
1370     #pragma omp for nowait
1371     for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1372     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1373     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1374     *o++ = 0.;
1375     *o++ = 1.;
1376     *o = 0.;
1377     }
1378     }
1379     }
1380    
1381     if (m_faceOffset[4] > -1) {
1382     #pragma omp for nowait
1383     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1384     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1385     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1386     *o++ = 0.;
1387     *o++ = 0.;
1388     *o = -1.;
1389     }
1390     }
1391     }
1392    
1393     if (m_faceOffset[5] > -1) {
1394     #pragma omp for nowait
1395     for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1396     for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1397     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1398     *o++ = 0.;
1399     *o++ = 0.;
1400     *o = 1.;
1401     }
1402     }
1403     }
1404     } // end of parallel section
1405    
1406     } else {
1407     stringstream msg;
1408     msg << "setToNormal() not implemented for "
1409     << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
1410     throw RipleyException(msg.str());
1411     }
1412     }
1413    
1414 caltinay 3691 Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
1415     bool reducedColOrder) const
1416     {
1417     if (reducedRowOrder || reducedColOrder)
1418     throw RipleyException("getPattern() not implemented for reduced order");
1419    
1420     throw RipleyException("getPattern() not implemented");
1421     }
1422    
1423     void Brick::Print_Mesh_Info(const bool full) const
1424     {
1425     RipleyDomain::Print_Mesh_Info(full);
1426     if (full) {
1427     cout << " Id Coordinates" << endl;
1428     cout.precision(15);
1429     cout.setf(ios::scientific, ios::floatfield);
1430 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
1431     pair<double,double> ydy = getFirstCoordAndSpacing(1);
1432     pair<double,double> zdz = getFirstCoordAndSpacing(2);
1433 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
1434     cout << " " << setw(5) << m_nodeId[i]
1435 caltinay 3698 << " " << xdx.first+(i%m_N0)*xdx.second
1436     << " " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
1437     << " " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
1438 caltinay 3691 }
1439     }
1440     }
1441    
1442 caltinay 3698 IndexVector Brick::getNumNodesPerDim() const
1443     {
1444     IndexVector ret;
1445     ret.push_back(m_N0);
1446     ret.push_back(m_N1);
1447     ret.push_back(m_N2);
1448     return ret;
1449     }
1450    
1451     IndexVector Brick::getNumElementsPerDim() const
1452     {
1453     IndexVector ret;
1454     ret.push_back(m_NE0);
1455     ret.push_back(m_NE1);
1456     ret.push_back(m_NE2);
1457     return ret;
1458     }
1459    
1460     IndexVector Brick::getNumFacesPerBoundary() const
1461     {
1462     IndexVector ret(6, 0);
1463     //left
1464     if (m_offset0==0)
1465     ret[0]=m_NE1*m_NE2;
1466     //right
1467     if (m_mpiInfo->rank%m_NX==m_NX-1)
1468     ret[1]=m_NE1*m_NE2;
1469     //bottom
1470     if (m_offset1==0)
1471     ret[2]=m_NE0*m_NE2;
1472     //top
1473     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
1474     ret[3]=m_NE0*m_NE2;
1475     //front
1476     if (m_offset2==0)
1477     ret[4]=m_NE0*m_NE1;
1478     //back
1479     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
1480     ret[5]=m_NE0*m_NE1;
1481     return ret;
1482     }
1483    
1484     pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
1485     {
1486     if (dim==0)
1487     return