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

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

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

revision 3702 by caltinay, Fri Dec 2 06:12:32 2011 UTC revision 3735 by caltinay, Fri Dec 9 05:03:08 2011 UTC
# Line 232  const int* Brick::borrowSampleReferenceI Line 232  const int* Brick::borrowSampleReferenceI
232  {  {
233      switch (fsType) {      switch (fsType) {
234          case Nodes:          case Nodes:
235            case DegreesOfFreedom: //FIXME
236            case ReducedDegreesOfFreedom: //FIXME
237              return &m_nodeId[0];              return &m_nodeId[0];
238          case Elements:          case Elements:
239            case ReducedElements:
240              return &m_elementId[0];              return &m_elementId[0];
241          case FaceElements:          case ReducedFaceElements:
242              return &m_faceId[0];              return &m_faceId[0];
243          default:          default:
244              break;              break;
# Line 261  bool Brick::ownSample(int fsCode, index_ Line 264  bool Brick::ownSample(int fsCode, index_
264  #endif  #endif
265  }  }
266    
267    void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const
268    {
269        escript::Data& in = *const_cast<escript::Data*>(&cIn);
270        const dim_t numComp = in.getDataPointSize();
271        const double h0 = m_l0/m_gNE0;
272        const double h1 = m_l1/m_gNE1;
273        const double h2 = m_l1/m_gNE2;
274        const double C0 = .044658198738520451079;
275        const double C1 = .16666666666666666667;
276        const double C2 = .21132486540518711775;
277        const double C3 = .25;
278        const double C4 = .5;
279        const double C5 = .62200846792814621559;
280        const double C6 = .78867513459481288225;
281    
282        if (out.getFunctionSpace().getTypeCode() == Elements) {
283            /*** GENERATOR SNIP_GRAD_ELEMENTS TOP */
284    #pragma omp parallel for
285            for (index_t k2=0; k2 < m_NE2; ++k2) {
286                for (index_t k1=0; k1 < m_NE1; ++k1) {
287                    for (index_t k0=0; k0 < m_NE0; ++k0) {
288                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
289                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
290                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
291                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
292                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
293                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
294                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
295                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
296                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
297                        for (index_t i=0; i < numComp; ++i) {
298                            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;
299                            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;
300                            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;
301                            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;
302                            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;
303                            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;
304                            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;
305                            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;
306                            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;
307                            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;
308                            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;
309                            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;
310                            o[INDEX3(i,0,0,numComp,3)] = V0;
311                            o[INDEX3(i,1,0,numComp,3)] = V4;
312                            o[INDEX3(i,2,0,numComp,3)] = V8;
313                            o[INDEX3(i,0,1,numComp,3)] = V0;
314                            o[INDEX3(i,1,1,numComp,3)] = V5;
315                            o[INDEX3(i,2,1,numComp,3)] = V9;
316                            o[INDEX3(i,0,2,numComp,3)] = V1;
317                            o[INDEX3(i,1,2,numComp,3)] = V4;
318                            o[INDEX3(i,2,2,numComp,3)] = V10;
319                            o[INDEX3(i,0,3,numComp,3)] = V1;
320                            o[INDEX3(i,1,3,numComp,3)] = V5;
321                            o[INDEX3(i,2,3,numComp,3)] = V11;
322                            o[INDEX3(i,0,4,numComp,3)] = V2;
323                            o[INDEX3(i,1,4,numComp,3)] = V6;
324                            o[INDEX3(i,2,4,numComp,3)] = V8;
325                            o[INDEX3(i,0,5,numComp,3)] = V2;
326                            o[INDEX3(i,1,5,numComp,3)] = V7;
327                            o[INDEX3(i,2,5,numComp,3)] = V9;
328                            o[INDEX3(i,0,6,numComp,3)] = V3;
329                            o[INDEX3(i,1,6,numComp,3)] = V6;
330                            o[INDEX3(i,2,6,numComp,3)] = V10;
331                            o[INDEX3(i,0,7,numComp,3)] = V3;
332                            o[INDEX3(i,1,7,numComp,3)] = V7;
333                            o[INDEX3(i,2,7,numComp,3)] = V11;
334                        } /* end of component loop i */
335                    } /* end of k0 loop */
336                } /* end of k1 loop */
337            } /* end of k2 loop */
338            /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
339        } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
340            /*** GENERATOR SNIP_GRAD_REDUCED_ELEMENTS TOP */
341    #pragma omp parallel for
342            for (index_t k2=0; k2 < m_NE2; ++k2) {
343                for (index_t k1=0; k1 < m_NE1; ++k1) {
344                    for (index_t k0=0; k0 < m_NE0; ++k0) {
345                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
346                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
347                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
348                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
349                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
350                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
351                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
352                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
353                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
354                        for (index_t i=0; i < numComp; ++i) {
355                            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;
356                            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;
357                            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;
358                        } /* end of component loop i */
359                    } /* end of k0 loop */
360                } /* end of k1 loop */
361            } /* end of k2 loop */
362            /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
363        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
364            /*** GENERATOR SNIP_GRAD_FACES TOP */
365    #pragma omp parallel
366            {
367                if (m_faceOffset[0] > -1) {
368    #pragma omp for nowait
369                    for (index_t k2=0; k2 < m_NE2; ++k2) {
370                        for (index_t k1=0; k1 < m_NE1; ++k1) {
371                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
372                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
373                            const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
374                            const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
375                            const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
376                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
377                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
378                            const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
379                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
380                            for (index_t i=0; i < numComp; ++i) {
381                                const double V0=((f_010[i]-f_000[i])*C6 + (f_011[i]-f_001[i])*C2) / h1;
382                                const double V1=((f_010[i]-f_000[i])*C2 + (f_011[i]-f_001[i])*C6) / h1;
383                                const double V2=((f_001[i]-f_000[i])*C6 + (f_010[i]-f_011[i])*C2) / h2;
384                                const double V3=((f_001[i]-f_000[i])*C2 + (f_011[i]-f_010[i])*C6) / h2;
385                                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;
386                                o[INDEX3(i,1,0,numComp,3)] = V0;
387                                o[INDEX3(i,2,0,numComp,3)] = V2;
388                                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;
389                                o[INDEX3(i,1,1,numComp,3)] = V0;
390                                o[INDEX3(i,2,1,numComp,3)] = V3;
391                                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;
392                                o[INDEX3(i,1,2,numComp,3)] = V1;
393                                o[INDEX3(i,2,2,numComp,3)] = V2;
394                                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;
395                                o[INDEX3(i,1,3,numComp,3)] = V1;
396                                o[INDEX3(i,2,3,numComp,3)] = V3;
397                            } /* end of component loop i */
398                        } /* end of k1 loop */
399                    } /* end of k2 loop */
400                } /* end of face 0 */
401                if (m_faceOffset[1] > -1) {
402    #pragma omp for nowait
403                    for (index_t k2=0; k2 < m_NE2; ++k2) {
404                        for (index_t k1=0; k1 < m_NE1; ++k1) {
405                            const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
406                            const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
407                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
408                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
409                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
410                            const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
411                            const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
412                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
413                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
414                            for (index_t i=0; i < numComp; ++i) {
415                                const double V0=((f_110[i]-f_100[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
416                                const double V1=((f_110[i]-f_100[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
417                                const double V2=((f_101[i]-f_100[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
418                                const double V3=((f_101[i]-f_100[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
419                                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;
420                                o[INDEX3(i,1,0,numComp,3)] = V0;
421                                o[INDEX3(i,2,0,numComp,3)] = V2;
422                                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;
423                                o[INDEX3(i,1,1,numComp,3)] = V0;
424                                o[INDEX3(i,2,1,numComp,3)] = V3;
425                                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;
426                                o[INDEX3(i,1,2,numComp,3)] = V1;
427                                o[INDEX3(i,2,2,numComp,3)] = V2;
428                                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;
429                                o[INDEX3(i,1,3,numComp,3)] = V1;
430                                o[INDEX3(i,2,3,numComp,3)] = V3;
431                            } /* end of component loop i */
432                        } /* end of k1 loop */
433                    } /* end of k2 loop */
434                } /* end of face 1 */
435                if (m_faceOffset[2] > -1) {
436    #pragma omp for nowait
437                    for (index_t k2=0; k2 < m_NE2; ++k2) {
438                        for (index_t k0=0; k0 < m_NE0; ++k0) {
439                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
440                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
441                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
442                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
443                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
444                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
445                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
446                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
447                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
448                            for (index_t i=0; i < numComp; ++i) {
449                                const double V0=((f_100[i]-f_000[i])*C6 + (f_101[i]-f_001[i])*C2) / h0;
450                                const double V1=((f_001[i]-f_000[i])*C6 + (f_101[i]-f_100[i])*C2) / h2;
451                                const double V2=((f_001[i]-f_000[i])*C2 + (f_101[i]-f_100[i])*C6) / h2;
452                                o[INDEX3(i,0,0,numComp,3)] = V0;
453                                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;
454                                o[INDEX3(i,2,0,numComp,3)] = V1;
455                                o[INDEX3(i,0,1,numComp,3)] = V0;
456                                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;
457                                o[INDEX3(i,2,1,numComp,3)] = V2;
458                                o[INDEX3(i,0,2,numComp,3)] = V0;
459                                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;
460                                o[INDEX3(i,2,2,numComp,3)] = V1;
461                                o[INDEX3(i,0,3,numComp,3)] = V0;
462                                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;
463                                o[INDEX3(i,2,3,numComp,3)] = V2;
464                            } /* end of component loop i */
465                        } /* end of k0 loop */
466                    } /* end of k2 loop */
467                } /* end of face 2 */
468                if (m_faceOffset[3] > -1) {
469    #pragma omp for nowait
470                    for (index_t k2=0; k2 < m_NE2; ++k2) {
471                        for (index_t k0=0; k0 < m_NE0; ++k0) {
472                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
473                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
474                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
475                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
476                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
477                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
478                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
479                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
480                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
481                            for (index_t i=0; i < numComp; ++i) {
482                                const double V0=((f_110[i]-f_010[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
483                                const double V1=((f_110[i]-f_010[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
484                                const double V2=((f_011[i]-f_010[i])*C6 + (f_111[i]-f_110[i])*C2) / h2;
485                                const double V3=((f_011[i]-f_010[i])*C2 + (f_111[i]-f_110[i])*C6) / h2;
486                                o[INDEX3(i,0,0,numComp,3)] = V0;
487                                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;
488                                o[INDEX3(i,2,0,numComp,3)] = V2;
489                                o[INDEX3(i,0,1,numComp,3)] = V0;
490                                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;
491                                o[INDEX3(i,2,1,numComp,3)] = V3;
492                                o[INDEX3(i,0,2,numComp,3)] = V1;
493                                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;
494                                o[INDEX3(i,2,2,numComp,3)] = V2;
495                                o[INDEX3(i,0,3,numComp,3)] = V1;
496                                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;
497                                o[INDEX3(i,2,3,numComp,3)] = V3;
498                            } /* end of component loop i */
499                        } /* end of k0 loop */
500                    } /* end of k2 loop */
501                } /* end of face 3 */
502                if (m_faceOffset[4] > -1) {
503    #pragma omp for nowait
504                    for (index_t k1=0; k1 < m_NE1; ++k1) {
505                        for (index_t k0=0; k0 < m_NE0; ++k0) {
506                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
507                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
508                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
509                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
510                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
511                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
512                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
513                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
514                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
515                            for (index_t i=0; i < numComp; ++i) {
516                                const double V0=((f_100[i]-f_000[i])*C6 + (f_110[i]-f_010[i])*C2) / h0;
517                                const double V1=((f_100[i]-f_000[i])*C2 + (f_110[i]-f_010[i])*C6) / h0;
518                                const double V2=((f_010[i]-f_000[i])*C6 + (f_110[i]-f_100[i])*C2) / h1;
519                                const double V3=((f_010[i]-f_000[i])*C2 + (f_110[i]-f_100[i])*C6) / h1;
520                                o[INDEX3(i,0,0,numComp,3)] = V0;
521                                o[INDEX3(i,1,0,numComp,3)] = V2;
522                                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;
523                                o[INDEX3(i,0,1,numComp,3)] = V0;
524                                o[INDEX3(i,1,1,numComp,3)] = V3;
525                                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;
526                                o[INDEX3(i,0,2,numComp,3)] = V1;
527                                o[INDEX3(i,1,2,numComp,3)] = V2;
528                                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;
529                                o[INDEX3(i,0,3,numComp,3)] = V1;
530                                o[INDEX3(i,1,3,numComp,3)] = V3;
531                                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;
532                            } /* end of component loop i */
533                        } /* end of k0 loop */
534                    } /* end of k1 loop */
535                } /* end of face 4 */
536                if (m_faceOffset[5] > -1) {
537    #pragma omp for nowait
538                    for (index_t k1=0; k1 < m_NE1; ++k1) {
539                        for (index_t k0=0; k0 < m_NE0; ++k0) {
540                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
541                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
542                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
543                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
544                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
545                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
546                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
547                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
548                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
549                            for (index_t i=0; i < numComp; ++i) {
550                                const double V0=((f_101[i]-f_001[i])*C6 + (f_111[i]-f_011[i])*C2) / h0;
551                                const double V1=((f_101[i]-f_001[i])*C2 + (f_111[i]-f_011[i])*C6) / h0;
552                                const double V2=((f_011[i]-f_001[i])*C6 + (f_111[i]-f_101[i])*C2) / h1;
553                                const double V3=((f_011[i]-f_001[i])*C2 + (f_111[i]-f_101[i])*C6) / h1;
554                                o[INDEX3(i,0,0,numComp,3)] = V0;
555                                o[INDEX3(i,1,0,numComp,3)] = V2;
556                                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;
557                                o[INDEX3(i,0,1,numComp,3)] = V0;
558                                o[INDEX3(i,1,1,numComp,3)] = V3;
559                                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;
560                                o[INDEX3(i,0,2,numComp,3)] = V1;
561                                o[INDEX3(i,1,2,numComp,3)] = V2;
562                                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;
563                                o[INDEX3(i,0,3,numComp,3)] = V1;
564                                o[INDEX3(i,1,3,numComp,3)] = V3;
565                                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;
566                            } /* end of component loop i */
567                        } /* end of k0 loop */
568                    } /* end of k1 loop */
569                } /* end of face 5 */
570            } // end of parallel section
571            /* GENERATOR SNIP_GRAD_FACES BOTTOM */
572        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
573            /*** GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
574    #pragma omp parallel
575            {
576                if (m_faceOffset[0] > -1) {
577    #pragma omp for nowait
578                    for (index_t k2=0; k2 < m_NE2; ++k2) {
579                        for (index_t k1=0; k1 < m_NE1; ++k1) {
580                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
581                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
582                            const register double* f_101 = in.getSampleDataRO(INDEX3(1,k1,k2+1, m_N0,m_N1));
583                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
584                            const register double* f_100 = in.getSampleDataRO(INDEX3(1,k1,k2, m_N0,m_N1));
585                            const register double* f_110 = in.getSampleDataRO(INDEX3(1,k1+1,k2, m_N0,m_N1));
586                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
587                            const register double* f_111 = in.getSampleDataRO(INDEX3(1,k1+1,k2+1, m_N0,m_N1));
588                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
589                            for (index_t i=0; i < numComp; ++i) {
590                                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;
591                                o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_011[i]-f_000[i]-f_001[i])*C4 / h1;
592                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_011[i]-f_000[i]-f_010[i])*C4 / h2;
593                            } /* end of component loop i */
594                        } /* end of k1 loop */
595                    } /* end of k2 loop */
596                } /* end of face 0 */
597                if (m_faceOffset[1] > -1) {
598    #pragma omp for nowait
599                    for (index_t k2=0; k2 < m_NE2; ++k2) {
600                        for (index_t k1=0; k1 < m_NE1; ++k1) {
601                            const register double* f_000 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2, m_N0,m_N1));
602                            const register double* f_001 = in.getSampleDataRO(INDEX3(m_N0-2,k1,k2+1, m_N0,m_N1));
603                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
604                            const register double* f_011 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2+1, m_N0,m_N1));
605                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
606                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
607                            const register double* f_010 = in.getSampleDataRO(INDEX3(m_N0-2,k1+1,k2, m_N0,m_N1));
608                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
609                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
610                            for (index_t i=0; i < numComp; ++i) {
611                                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;
612                                o[INDEX3(i,1,0,numComp,3)] = (f_110[i]+f_111[i]-f_100[i]-f_101[i])*C4 / h1;
613                                o[INDEX3(i,2,0,numComp,3)] = (f_101[i]+f_111[i]-f_100[i]-f_110[i])*C4 / h2;
614                            } /* end of component loop i */
615                        } /* end of k1 loop */
616                    } /* end of k2 loop */
617                } /* end of face 1 */
618                if (m_faceOffset[2] > -1) {
619    #pragma omp for nowait
620                    for (index_t k2=0; k2 < m_NE2; ++k2) {
621                        for (index_t k0=0; k0 < m_NE0; ++k0) {
622                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
623                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
624                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
625                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
626                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,1,k2+1, m_N0,m_N1));
627                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,1,k2, m_N0,m_N1));
628                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,1,k2, m_N0,m_N1));
629                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,1,k2+1, m_N0,m_N1));
630                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
631                            for (index_t i=0; i < numComp; ++i) {
632                                o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_101[i]-f_000[i]-f_001[i])*C4 / h0;
633                                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;
634                                o[INDEX3(i,2,0,numComp,3)] = (f_001[i]+f_101[i]-f_000[i]-f_100[i])*C4 / h2;
635                            } /* end of component loop i */
636                        } /* end of k0 loop */
637                    } /* end of k2 loop */
638                } /* end of face 2 */
639                if (m_faceOffset[3] > -1) {
640    #pragma omp for nowait
641                    for (index_t k2=0; k2 < m_NE2; ++k2) {
642                        for (index_t k0=0; k0 < m_NE0; ++k0) {
643                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
644                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
645                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
646                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
647                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2, m_N0,m_N1));
648                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2+1, m_N0,m_N1));
649                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,m_N1-2,k2+1, m_N0,m_N1));
650                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,m_N1-2,k2, m_N0,m_N1));
651                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
652                            for (index_t i=0; i < numComp; ++i) {
653                                o[INDEX3(i,0,0,numComp,3)] = (f_110[i]+f_111[i]-f_010[i]-f_011[i])*C4 / h0;
654                                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;
655                                o[INDEX3(i,2,0,numComp,3)] = (f_011[i]+f_111[i]-f_010[i]-f_110[i])*C4 / h2;
656                            } /* end of component loop i */
657                        } /* end of k0 loop */
658                    } /* end of k2 loop */
659                } /* end of face 3 */
660                if (m_faceOffset[4] > -1) {
661    #pragma omp for nowait
662                    for (index_t k1=0; k1 < m_NE1; ++k1) {
663                        for (index_t k0=0; k0 < m_NE0; ++k0) {
664                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
665                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
666                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
667                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
668                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,1, m_N0,m_N1));
669                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,1, m_N0,m_N1));
670                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,1, m_N0,m_N1));
671                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,1, m_N0,m_N1));
672                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
673                            for (index_t i=0; i < numComp; ++i) {
674                                o[INDEX3(i,0,0,numComp,3)] = (f_100[i]+f_110[i]-f_000[i]-f_010[i])*C4 / h0;
675                                o[INDEX3(i,1,0,numComp,3)] = (f_010[i]+f_110[i]-f_000[i]-f_100[i])*C4 / h1;
676                                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;
677                            } /* end of component loop i */
678                        } /* end of k0 loop */
679                    } /* end of k1 loop */
680                } /* end of face 4 */
681                if (m_faceOffset[5] > -1) {
682    #pragma omp for nowait
683                    for (index_t k1=0; k1 < m_NE1; ++k1) {
684                        for (index_t k0=0; k0 < m_NE0; ++k0) {
685                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
686                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
687                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
688                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
689                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-2, m_N0,m_N1));
690                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-2, m_N0,m_N1));
691                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-2, m_N0,m_N1));
692                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-2, m_N0,m_N1));
693                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
694                            for (index_t i=0; i < numComp; ++i) {
695                                o[INDEX3(i,0,0,numComp,3)] = (f_101[i]+f_111[i]-f_001[i]-f_011[i])*C4 / h0;
696                                o[INDEX3(i,1,0,numComp,3)] = (f_011[i]+f_111[i]-f_001[i]-f_101[i])*C4 / h1;
697                                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;
698                            } /* end of component loop i */
699                        } /* end of k0 loop */
700                    } /* end of k1 loop */
701                } /* end of face 5 */
702            } // end of parallel section
703            /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
704        } else {
705            stringstream msg;
706            msg << "setToGradient() not implemented for "
707                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
708            throw RipleyException(msg.str());
709        }
710    }
711    
712    void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
713    {
714        escript::Data& in = *const_cast<escript::Data*>(&arg);
715        const dim_t numComp = in.getDataPointSize();
716        const double h0 = m_l0/m_gNE0;
717        const double h1 = m_l1/m_gNE1;
718        const double h2 = m_l2/m_gNE2;
719        if (arg.getFunctionSpace().getTypeCode() == Elements) {
720            const double w_0 = h0*h1*h2/8.;
721    #pragma omp parallel
722            {
723                vector<double> int_local(numComp, 0);
724    #pragma omp for nowait
725                for (index_t k2 = 0; k2 < m_NE2; ++k2) {
726                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
727                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
728                            const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
729                            for (index_t i=0; i < numComp; ++i) {
730                                const register double f_0 = f[INDEX2(i,0,numComp)];
731                                const register double f_1 = f[INDEX2(i,1,numComp)];
732                                const register double f_2 = f[INDEX2(i,2,numComp)];
733                                const register double f_3 = f[INDEX2(i,3,numComp)];
734                                const register double f_4 = f[INDEX2(i,4,numComp)];
735                                const register double f_5 = f[INDEX2(i,5,numComp)];
736                                const register double f_6 = f[INDEX2(i,6,numComp)];
737                                const register double f_7 = f[INDEX2(i,7,numComp)];
738                                int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
739                            }  /* end of component loop i */
740                        } /* end of k0 loop */
741                    } /* end of k1 loop */
742                } /* end of k2 loop */
743    
744    #pragma omp critical
745                for (index_t i=0; i<numComp; i++)
746                    integrals[i]+=int_local[i];
747            } // end of parallel section
748        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
749            const double w_0 = h0*h1*h2;
750    #pragma omp parallel
751            {
752                vector<double> int_local(numComp, 0);
753    #pragma omp for nowait
754                for (index_t k2 = 0; k2 < m_NE2; ++k2) {
755                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
756                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
757                            const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
758                            for (index_t i=0; i < numComp; ++i) {
759                                int_local[i]+=f[i]*w_0;
760                            }  /* end of component loop i */
761                        } /* end of k0 loop */
762                    } /* end of k1 loop */
763                } /* end of k2 loop */
764    
765    #pragma omp critical
766                for (index_t i=0; i<numComp; i++)
767                    integrals[i]+=int_local[i];
768            } // end of parallel section
769        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
770            const double w_0 = h1*h2/4.;
771            const double w_1 = h0*h2/4.;
772            const double w_2 = h0*h1/4.;
773    #pragma omp parallel
774            {
775                vector<double> int_local(numComp, 0);
776                if (m_faceOffset[0] > -1) {
777    #pragma omp for nowait
778                    for (index_t k2=0; k2 < m_NE2; ++k2) {
779                        for (index_t k1=0; k1 < m_NE1; ++k1) {
780                            const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
781                            for (index_t i=0; i < numComp; ++i) {
782                                const register double f_0 = f[INDEX2(i,0,numComp)];
783                                const register double f_1 = f[INDEX2(i,1,numComp)];
784                                const register double f_2 = f[INDEX2(i,2,numComp)];
785                                const register double f_3 = f[INDEX2(i,3,numComp)];
786                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
787                            }  /* end of component loop i */
788                        } /* end of k1 loop */
789                    } /* end of k2 loop */
790                }
791    
792                if (m_faceOffset[1] > -1) {
793    #pragma omp for nowait
794                    for (index_t k2=0; k2 < m_NE2; ++k2) {
795                        for (index_t k1=0; k1 < m_NE1; ++k1) {
796                            const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
797                            for (index_t i=0; i < numComp; ++i) {
798                                const register double f_0 = f[INDEX2(i,0,numComp)];
799                                const register double f_1 = f[INDEX2(i,1,numComp)];
800                                const register double f_2 = f[INDEX2(i,2,numComp)];
801                                const register double f_3 = f[INDEX2(i,3,numComp)];
802                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
803                            }  /* end of component loop i */
804                        } /* end of k1 loop */
805                    } /* end of k2 loop */
806                }
807    
808                if (m_faceOffset[2] > -1) {
809    #pragma omp for nowait
810                    for (index_t k2=0; k2 < m_NE2; ++k2) {
811                        for (index_t k0=0; k0 < m_NE0; ++k0) {
812                            const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
813                            for (index_t i=0; i < numComp; ++i) {
814                                const register double f_0 = f[INDEX2(i,0,numComp)];
815                                const register double f_1 = f[INDEX2(i,1,numComp)];
816                                const register double f_2 = f[INDEX2(i,2,numComp)];
817                                const register double f_3 = f[INDEX2(i,3,numComp)];
818                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
819                            }  /* end of component loop i */
820                        } /* end of k1 loop */
821                    } /* end of k2 loop */
822                }
823    
824                if (m_faceOffset[3] > -1) {
825    #pragma omp for nowait
826                    for (index_t k2=0; k2 < m_NE2; ++k2) {
827                        for (index_t k0=0; k0 < m_NE0; ++k0) {
828                            const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
829                            for (index_t i=0; i < numComp; ++i) {
830                                const register double f_0 = f[INDEX2(i,0,numComp)];
831                                const register double f_1 = f[INDEX2(i,1,numComp)];
832                                const register double f_2 = f[INDEX2(i,2,numComp)];
833                                const register double f_3 = f[INDEX2(i,3,numComp)];
834                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
835                            }  /* end of component loop i */
836                        } /* end of k1 loop */
837                    } /* end of k2 loop */
838                }
839    
840                if (m_faceOffset[4] > -1) {
841    #pragma omp for nowait
842                    for (index_t k1=0; k1 < m_NE1; ++k1) {
843                        for (index_t k0=0; k0 < m_NE0; ++k0) {
844                            const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
845                            for (index_t i=0; i < numComp; ++i) {
846                                const register double f_0 = f[INDEX2(i,0,numComp)];
847                                const register double f_1 = f[INDEX2(i,1,numComp)];
848                                const register double f_2 = f[INDEX2(i,2,numComp)];
849                                const register double f_3 = f[INDEX2(i,3,numComp)];
850                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
851                            }  /* end of component loop i */
852                        } /* end of k1 loop */
853                    } /* end of k2 loop */
854                }
855    
856                if (m_faceOffset[5] > -1) {
857    #pragma omp for nowait
858                    for (index_t k1=0; k1 < m_NE1; ++k1) {
859                        for (index_t k0=0; k0 < m_NE0; ++k0) {
860                            const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
861                            for (index_t i=0; i < numComp; ++i) {
862                                const register double f_0 = f[INDEX2(i,0,numComp)];
863                                const register double f_1 = f[INDEX2(i,1,numComp)];
864                                const register double f_2 = f[INDEX2(i,2,numComp)];
865                                const register double f_3 = f[INDEX2(i,3,numComp)];
866                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
867                            }  /* end of component loop i */
868                        } /* end of k1 loop */
869                    } /* end of k2 loop */
870                }
871    
872    #pragma omp critical
873                for (index_t i=0; i<numComp; i++)
874                    integrals[i]+=int_local[i];
875            } // end of parallel section
876    
877        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
878            const double w_0 = h1*h2;
879            const double w_1 = h0*h2;
880            const double w_2 = h0*h1;
881    #pragma omp parallel
882            {
883                vector<double> int_local(numComp, 0);
884                if (m_faceOffset[0] > -1) {
885    #pragma omp for nowait
886                    for (index_t k2=0; k2 < m_NE2; ++k2) {
887                        for (index_t k1=0; k1 < m_NE1; ++k1) {
888                            const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
889                            for (index_t i=0; i < numComp; ++i) {
890                                int_local[i]+=f[i]*w_0;
891                            }  /* end of component loop i */
892                        } /* end of k1 loop */
893                    } /* end of k2 loop */
894                }
895    
896                if (m_faceOffset[1] > -1) {
897    #pragma omp for nowait
898                    for (index_t k2=0; k2 < m_NE2; ++k2) {
899                        for (index_t k1=0; k1 < m_NE1; ++k1) {
900                            const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
901                            for (index_t i=0; i < numComp; ++i) {
902                                int_local[i]+=f[i]*w_0;
903                            }  /* end of component loop i */
904                        } /* end of k1 loop */
905                    } /* end of k2 loop */
906                }
907    
908                if (m_faceOffset[2] > -1) {
909    #pragma omp for nowait
910                    for (index_t k2=0; k2 < m_NE2; ++k2) {
911                        for (index_t k0=0; k0 < m_NE0; ++k0) {
912                            const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
913                            for (index_t i=0; i < numComp; ++i) {
914                                int_local[i]+=f[i]*w_1;
915                            }  /* end of component loop i */
916                        } /* end of k1 loop */
917                    } /* end of k2 loop */
918                }
919    
920                if (m_faceOffset[3] > -1) {
921    #pragma omp for nowait
922                    for (index_t k2=0; k2 < m_NE2; ++k2) {
923                        for (index_t k0=0; k0 < m_NE0; ++k0) {
924                            const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
925                            for (index_t i=0; i < numComp; ++i) {
926                                int_local[i]+=f[i]*w_1;
927                            }  /* end of component loop i */
928                        } /* end of k1 loop */
929                    } /* end of k2 loop */
930                }
931    
932                if (m_faceOffset[4] > -1) {
933    #pragma omp for nowait
934                    for (index_t k1=0; k1 < m_NE1; ++k1) {
935                        for (index_t k0=0; k0 < m_NE0; ++k0) {
936                            const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
937                            for (index_t i=0; i < numComp; ++i) {
938                                int_local[i]+=f[i]*w_2;
939                            }  /* end of component loop i */
940                        } /* end of k1 loop */
941                    } /* end of k2 loop */
942                }
943    
944                if (m_faceOffset[5] > -1) {
945    #pragma omp for nowait
946                    for (index_t k1=0; k1 < m_NE1; ++k1) {
947                        for (index_t k0=0; k0 < m_NE0; ++k0) {
948                            const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
949                            for (index_t i=0; i < numComp; ++i) {
950                                int_local[i]+=f[i]*w_2;
951                            }  /* end of component loop i */
952                        } /* end of k1 loop */
953                    } /* end of k2 loop */
954                }
955    
956    #pragma omp critical
957                for (index_t i=0; i<numComp; i++)
958                    integrals[i]+=int_local[i];
959            } // end of parallel section
960    
961        } else {
962            stringstream msg;
963            msg << "setToIntegrals() not implemented for "
964                << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
965            throw RipleyException(msg.str());
966        }
967    }
968    
969    void Brick::setToNormal(escript::Data& out) const
970    {
971        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
972    #pragma omp parallel
973            {
974                if (m_faceOffset[0] > -1) {
975    #pragma omp for nowait
976                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
977                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
978                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
979                            // set vector at four quadrature points
980                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
981                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
982                            *o++ = -1.; *o++ = 0.; *o++ = 0.;
983                            *o++ = -1.; *o++ = 0.; *o = 0.;
984                        }
985                    }
986                }
987    
988                if (m_faceOffset[1] > -1) {
989    #pragma omp for nowait
990                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
991                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
992                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
993                            // set vector at four quadrature points
994                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
995                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
996                            *o++ = 1.; *o++ = 0.; *o++ = 0.;
997                            *o++ = 1.; *o++ = 0.; *o = 0.;
998                        }
999                    }
1000                }
1001    
1002                if (m_faceOffset[2] > -1) {
1003    #pragma omp for nowait
1004                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1005                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1006                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1007                            // set vector at four quadrature points
1008                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1009                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1010                            *o++ = 0.; *o++ = -1.; *o++ = 0.;
1011                            *o++ = 0.; *o++ = -1.; *o = 0.;
1012                        }
1013                    }
1014                }
1015    
1016                if (m_faceOffset[3] > -1) {
1017    #pragma omp for nowait
1018                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1019                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1020                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1021                            // set vector at four quadrature points
1022                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1023                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1024                            *o++ = 0.; *o++ = 1.; *o++ = 0.;
1025                            *o++ = 0.; *o++ = 1.; *o = 0.;
1026                        }
1027                    }
1028                }
1029    
1030                if (m_faceOffset[4] > -1) {
1031    #pragma omp for nowait
1032                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1033                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1034                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1035                            // set vector at four quadrature points
1036                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1037                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1038                            *o++ = 0.; *o++ = 0.; *o++ = -1.;
1039                            *o++ = 0.; *o++ = 0.; *o = -1.;
1040                        }
1041                    }
1042                }
1043    
1044                if (m_faceOffset[5] > -1) {
1045    #pragma omp for nowait
1046                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1047                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1048                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1049                            // set vector at four quadrature points
1050                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1051                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1052                            *o++ = 0.; *o++ = 0.; *o++ = 1.;
1053                            *o++ = 0.; *o++ = 0.; *o = 1.;
1054                        }
1055                    }
1056                }
1057            } // end of parallel section
1058        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1059    #pragma omp parallel
1060            {
1061                if (m_faceOffset[0] > -1) {
1062    #pragma omp for nowait
1063                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1064                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1065                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1066                            *o++ = -1.;
1067                            *o++ = 0.;
1068                            *o = 0.;
1069                        }
1070                    }
1071                }
1072    
1073                if (m_faceOffset[1] > -1) {
1074    #pragma omp for nowait
1075                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1076                        for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1077                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1078                            *o++ = 1.;
1079                            *o++ = 0.;
1080                            *o = 0.;
1081                        }
1082                    }
1083                }
1084    
1085                if (m_faceOffset[2] > -1) {
1086    #pragma omp for nowait
1087                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1088                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1089                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1090                            *o++ = 0.;
1091                            *o++ = -1.;
1092                            *o = 0.;
1093                        }
1094                    }
1095                }
1096    
1097                if (m_faceOffset[3] > -1) {
1098    #pragma omp for nowait
1099                    for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1100                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1101                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1102                            *o++ = 0.;
1103                            *o++ = 1.;
1104                            *o = 0.;
1105                        }
1106                    }
1107                }
1108    
1109                if (m_faceOffset[4] > -1) {
1110    #pragma omp for nowait
1111                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1112                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1113                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1114                            *o++ = 0.;
1115                            *o++ = 0.;
1116                            *o = -1.;
1117                        }
1118                    }
1119                }
1120    
1121                if (m_faceOffset[5] > -1) {
1122    #pragma omp for nowait
1123                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1124                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1125                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1126                            *o++ = 0.;
1127                            *o++ = 0.;
1128                            *o = 1.;
1129                        }
1130                    }
1131                }
1132            } // end of parallel section
1133    
1134        } else {
1135            stringstream msg;
1136            msg << "setToNormal() not implemented for "
1137                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
1138            throw RipleyException(msg.str());
1139        }
1140    }
1141    
1142  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
1143                                              bool reducedColOrder) const                                              bool reducedColOrder) const
1144  {  {
# Line 540  void Brick::populateSampleIds() Line 1418  void Brick::populateSampleIds()
1418              }              }
1419          }          }
1420      }      }
1421        m_nodeTags.assign(getNumNodes(), 0);
1422        updateTagsInUse(Nodes);
1423    
1424      // elements      // elements
1425      m_elementId.resize(getNumElements());      m_elementId.resize(getNumElements());
# Line 547  void Brick::populateSampleIds() Line 1427  void Brick::populateSampleIds()
1427      for (dim_t k=0; k<getNumElements(); k++) {      for (dim_t k=0; k<getNumElements(); k++) {
1428          m_elementId[k]=k;          m_elementId[k]=k;
1429      }      }
1430        m_elementTags.assign(getNumElements(), 0);
1431        updateTagsInUse(Elements);
1432    
1433      // face elements      // face elements
1434      m_faceId.resize(getNumFaceElements());      m_faceId.resize(getNumFaceElements());
# Line 554  void Brick::populateSampleIds() Line 1436  void Brick::populateSampleIds()
1436      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
1437          m_faceId[k]=k;          m_faceId[k]=k;
1438      }      }
1439    
1440        // generate face offset vector and set face tags
1441        const IndexVector facesPerEdge = getNumFacesPerBoundary();
1442        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
1443        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
1444        m_faceOffset.assign(facesPerEdge.size(), -1);
1445        m_faceTags.clear();
1446        index_t offset=0;
1447        for (size_t i=0; i<facesPerEdge.size(); i++) {
1448            if (facesPerEdge[i]>0) {
1449                m_faceOffset[i]=offset;
1450                offset+=facesPerEdge[i];
1451                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1452            }
1453        }
1454        setTagMap("left", LEFT);
1455        setTagMap("right", RIGHT);
1456        setTagMap("bottom", BOTTOM);
1457        setTagMap("top", TOP);
1458        setTagMap("front", FRONT);
1459        setTagMap("back", BACK);
1460        updateTagsInUse(FaceElements);
1461    }
1462    
1463    //protected
1464    void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,
1465                                           bool reduced) const
1466    {
1467        const dim_t numComp = in.getDataPointSize();
1468        if (reduced) {
1469            /*** GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS TOP */
1470            const double c0 = .125;
1471    #pragma omp parallel for
1472            for (index_t k2=0; k2 < m_NE2; ++k2) {
1473                for (index_t k1=0; k1 < m_NE1; ++k1) {
1474                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1475                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1476                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1477                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1478                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1479                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1480                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1481                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1482                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1483                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1484                        for (index_t i=0; i < numComp; ++i) {
1485                            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]);
1486                        } /* end of component loop i */
1487                    } /* end of k0 loop */
1488                } /* end of k1 loop */
1489            } /* end of k2 loop */
1490            /* GENERATOR SNIP_INTERPOLATE_REDUCED_ELEMENTS BOTTOM */
1491        } else {
1492            /*** GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
1493            const double c0 = .0094373878376559314545;
1494            const double c1 = .035220810900864519624;
1495            const double c2 = .13144585576580214704;
1496            const double c3 = .49056261216234406855;
1497    #pragma omp parallel for
1498            for (index_t k2=0; k2 < m_NE2; ++k2) {
1499                for (index_t k1=0; k1 < m_NE1; ++k1) {
1500                    for (index_t k0=0; k0 < m_NE0; ++k0) {
1501                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
1502                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
1503                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
1504                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
1505                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
1506                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
1507                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
1508                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
1509                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
1510                        for (index_t i=0; i < numComp; ++i) {
1511                            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]);
1512                            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]);
1513                            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]);
1514                            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]);
1515                            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]);
1516                            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]);
1517                            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]);
1518                            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]);
1519                        } /* end of component loop i */
1520                    } /* end of k0 loop */
1521                } /* end of k1 loop */
1522            } /* end of k2 loop */
1523            /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
1524        }
1525    }
1526    
1527    //protected
1528    void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,
1529                                        bool reduced) const
1530    {
1531        const dim_t numComp = in.getDataPointSize();
1532        if (reduced) {
1533            const double c0 = .25;
1534    #pragma omp parallel
1535            {
1536                /*** GENERATOR SNIP_INTERPOLATE_REDUCED_FACES TOP */
1537                if (m_faceOffset[0] > -1) {
1538    #pragma omp for nowait
1539                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1540                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1541                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1542                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1543                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1544                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1545                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1546                            for (index_t i=0; i < numComp; ++i) {
1547                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_010[i] + f_011[i]);
1548                            } /* end of component loop i */
1549                        } /* end of k1 loop */
1550                    } /* end of k2 loop */
1551                } /* end of face 0 */
1552                if (m_faceOffset[1] > -1) {
1553    #pragma omp for nowait
1554                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1555                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1556                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1557                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1558                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1559                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1560                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1561                            for (index_t i=0; i < numComp; ++i) {
1562                                o[INDEX2(i,numComp,0)] = c0*(f_100[i] + f_101[i] + f_110[i] + f_111[i]);
1563                            } /* end of component loop i */
1564                        } /* end of k1 loop */
1565                    } /* end of k2 loop */
1566                } /* end of face 1 */
1567                if (m_faceOffset[2] > -1) {
1568    #pragma omp for nowait
1569                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1570                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1571                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1572                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1573                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1574                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1575                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1576                            for (index_t i=0; i < numComp; ++i) {
1577                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_001[i] + f_100[i] + f_101[i]);
1578                            } /* end of component loop i */
1579                        } /* end of k0 loop */
1580                    } /* end of k2 loop */
1581                } /* end of face 2 */
1582                if (m_faceOffset[3] > -1) {
1583    #pragma omp for nowait
1584                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1585                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1586                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1587                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1588                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1589                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1590                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1591                            for (index_t i=0; i < numComp; ++i) {
1592                                o[INDEX2(i,numComp,0)] = c0*(f_010[i] + f_011[i] + f_110[i] + f_111[i]);
1593                            } /* end of component loop i */
1594                        } /* end of k0 loop */
1595                    } /* end of k2 loop */
1596                } /* end of face 3 */
1597                if (m_faceOffset[4] > -1) {
1598    #pragma omp for nowait
1599                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1600                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1601                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1602                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1603                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1604                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
1605                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1606                            for (index_t i=0; i < numComp; ++i) {
1607                                o[INDEX2(i,numComp,0)] = c0*(f_000[i] + f_010[i] + f_100[i] + f_110[i]);
1608                            } /* end of component loop i */
1609                        } /* end of k0 loop */
1610                    } /* end of k1 loop */
1611                } /* end of face 4 */
1612                if (m_faceOffset[5] > -1) {
1613    #pragma omp for nowait
1614                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1615                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1616                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
1617                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1618                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
1619                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
1620                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1621                            for (index_t i=0; i < numComp; ++i) {
1622                                o[INDEX2(i,numComp,0)] = c0*(f_001[i] + f_011[i] + f_101[i] + f_111[i]);
1623                            } /* end of component loop i */
1624                        } /* end of k0 loop */
1625                    } /* end of k1 loop */
1626                } /* end of face 5 */
1627                /* GENERATOR SNIP_INTERPOLATE_REDUCED_FACES BOTTOM */
1628            } // end of parallel section
1629        } else {
1630            const double c0 = 0.044658198738520451079;
1631            const double c1 = 0.16666666666666666667;
1632            const double c2 = 0.62200846792814621559;
1633    #pragma omp parallel
1634            {
1635                /*** GENERATOR SNIP_INTERPOLATE_FACES TOP */
1636                if (m_faceOffset[0] > -1) {
1637    #pragma omp for nowait
1638                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1639                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1640                            const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
1641                            const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
1642                            const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
1643                            const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
1644                            double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1645                            for (index_t i=0; i < numComp; ++i) {
1646                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_011[i]*c0 + c1*(f_001[i] + f_010[i]);
1647                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_010[i]*c2 + c1*(f_000[i] + f_011[i]);
1648                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_010[i]*c0 + c1*(f_000[i] + f_011[i]);
1649                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_011[i]*c2 + c1*(f_001[i] + f_010[i]);
1650                            } /* end of component loop i */
1651                        } /* end of k1 loop */
1652                    } /* end of k2 loop */
1653                } /* end of face 0 */
1654                if (m_faceOffset[1] > -1) {
1655    #pragma omp for nowait
1656                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1657                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1658                            const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
1659                            const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
1660                            const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
1661                            const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
1662                            double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1663                            for (index_t i=0; i < numComp; ++i) {
1664                                o[INDEX2(i,numComp,0)] = f_100[i]*c2 + f_111[i]*c0 + c1*(f_101[i] + f_110[i]);
1665                                o[INDEX2(i,numComp,1)] = f_101[i]*c0 + f_110[i]*c2 + c1*(f_100[i] + f_111[i]);
1666                                o[INDEX2(i,numComp,2)] = f_101[i]*c2 + f_110[i]*c0 + c1*(f_100[i] + f_111[i]);
1667                                o[INDEX2(i,numComp,3)] = f_100[i]*c0 + f_111[i]*c2 + c1*(f_101[i] + f_110[i]);
1668                            } /* end of component loop i */
1669                        } /* end of k1 loop */
1670                    } /* end of k2 loop */
1671                } /* end of face 1 */
1672                if (m_faceOffset[2] > -1) {
1673    #pragma omp for nowait
1674                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1675                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1676                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
1677                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
1678                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
1679                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
1680                            double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1681                            for (index_t i=0; i < numComp; ++i) {
1682                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_100[i]);
1683                                o[INDEX2(i,numComp,1)] = f_001[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_101[i]);
1684                                o[INDEX2(i,numComp,2)] = f_001[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_101[i]);
1685                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_100[i]);
1686                            } /* end of component loop i */
1687                        } /* end of k0 loop */
1688                    } /* end of k2 loop */
1689                } /* end of face 2 */
1690                if (m_faceOffset[3] > -1) {
1691    #pragma omp for nowait
1692                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1693                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1694                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
1695                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
1696                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
1697                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
1698                            double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1699                            for (index_t i=0; i < numComp; ++i) {
1700                                o[INDEX2(i,numComp,0)] = f_010[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_110[i]);
1701                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_111[i]);
1702                                o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_111[i]);
1703                                o[INDEX2(i,numComp,3)] = f_010[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_110[i]);
1704                            } /* end of component loop i */
1705                        } /* end of k0 loop */
1706                    } /* end of k2 loop */
1707                } /* end of face 3 */
1708                if (m_faceOffset[4] > -1) {
1709    #pragma omp for nowait
1710                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1711                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1712                            const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
1713                            const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
1714                            const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
1715                            const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
1716                            double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1717                            for (index_t i=0; i < numComp; ++i) {
1718                                o[INDEX2(i,numComp,0)] = f_000[i]*c2 + f_110[i]*c0 + c1*(f_010[i] + f_100[i]);
1719                                o[INDEX2(i,numComp,1)] = f_010[i]*c0 + f_100[i]*c2 + c1*(f_000[i] + f_110[i]);
1720                                o[INDEX2(i,numComp,2)] = f_010[i]*c2 + f_100[i]*c0 + c1*(f_000[i] + f_110[i]);
1721                                o[INDEX2(i,numComp,3)] = f_000[i]*c0 + f_110[i]*c2 + c1*(f_010[i] + f_100[i]);
1722                            } /* end of component loop i */
1723                        } /* end of k0 loop */
1724                    } /* end of k1 loop */
1725                } /* end of face 4 */
1726                if (m_faceOffset[5] > -1) {
1727    #pragma omp for nowait
1728                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1729                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1730                            const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
1731                            const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
1732                            const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
1733                            const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
1734                            double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1735                            for (index_t i=0; i < numComp; ++i) {
1736                                o[INDEX2(i,numComp,0)] = f_001[i]*c2 + f_111[i]*c0 + c1*(f_011[i] + f_101[i]);
1737                                o[INDEX2(i,numComp,1)] = f_011[i]*c0 + f_101[i]*c2 + c1*(f_001[i] + f_111[i]);
1738                                o[INDEX2(i,numComp,2)] = f_011[i]*c2 + f_101[i]*c0 + c1*(f_001[i] + f_111[i]);
1739                                o[INDEX2(i,numComp,3)] = f_001[i]*c0 + f_111[i]*c2 + c1*(f_011[i] + f_101[i]);
1740                            } /* end of component loop i */
1741                        } /* end of k0 loop */
1742                    } /* end of k1 loop */
1743                } /* end of face 5 */
1744                /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
1745            } // end of parallel section
1746        }
1747  }  }
1748    
1749  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3702  
changed lines
  Added in v.3735

  ViewVC Help
Powered by ViewVC 1.1.26