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

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

  ViewVC Help
Powered by ViewVC 1.1.26