/[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 3704 by caltinay, Mon Dec 5 01:59:08 2011 UTC
# Line 261  bool Brick::ownSample(int fsCode, index_ Line 261  bool Brick::ownSample(int fsCode, index_
261  #endif  #endif
262  }  }
263    
264    void Brick::setToGradient(escript::Data& out, const escript::Data& cIn) const
265    {
266        escript::Data& in = *const_cast<escript::Data*>(&cIn);
267        const dim_t numComp = in.getDataPointSize();
268        const double h0 = m_l0/m_gNE0;
269        const double h1 = m_l1/m_gNE1;
270        const double h2 = m_l1/m_gNE2;
271        if (out.getFunctionSpace().getTypeCode() == Elements) {
272            /* GENERATOR SNIP_GRAD_ELEMENTS TOP */
273            const double tmp0_22 = -0.044658198738520451079/h1;
274            const double tmp0_16 = 0.16666666666666666667/h0;
275            const double tmp0_33 = -0.62200846792814621559/h1;
276            const double tmp0_0 = -0.62200846792814621559/h0;
277            const double tmp0_21 = -0.16666666666666666667/h1;
278            const double tmp0_17 = 0.62200846792814621559/h0;
279            const double tmp0_52 = -0.044658198738520451079/h2;
280            const double tmp0_1 = -0.16666666666666666667/h0;
281            const double tmp0_20 = -0.62200846792814621559/h1;
282            const double tmp0_14 = -0.044658198738520451079/h0;
283            const double tmp0_53 = -0.62200846792814621559/h2;
284            const double tmp0_49 = 0.16666666666666666667/h2;
285            const double tmp0_2 = 0.16666666666666666667/h0;
286            const double tmp0_27 = -0.044658198738520451079/h1;
287            const double tmp0_15 = -0.16666666666666666667/h0;
288            const double tmp0_50 = -0.16666666666666666667/h2;
289            const double tmp0_48 = 0.62200846792814621559/h2;
290            const double tmp0_3 = 0.044658198738520451079/h0;
291            const double tmp0_26 = -0.16666666666666666667/h1;
292            const double tmp0_12 = -0.62200846792814621559/h0;
293            const double tmp0_51 = 0.044658198738520451079/h2;
294            const double tmp0_25 = 0.62200846792814621559/h1;
295            const double tmp0_13 = 0.16666666666666666667/h0;
296            const double tmp0_56 = 0.16666666666666666667/h2;
297            const double tmp0_24 = 0.16666666666666666667/h1;
298            const double tmp0_10 = 0.62200846792814621559/h0;
299            const double tmp0_57 = 0.62200846792814621559/h2;
300            const double tmp0_11 = -0.16666666666666666667/h0;
301            const double tmp0_54 = -0.044658198738520451079/h2;
302            const double tmp0_38 = 0.16666666666666666667/h1;
303            const double tmp0_34 = -0.044658198738520451079/h1;
304            const double tmp0_42 = 0.16666666666666666667/h2;
305            const double tmp0_35 = -0.16666666666666666667/h1;
306            const double tmp0_36 = -0.62200846792814621559/h1;
307            const double tmp0_41 = 0.62200846792814621559/h2;
308            const double tmp0_8 = 0.044658198738520451079/h0;
309            const double tmp0_37 = 0.62200846792814621559/h1;
310            const double tmp0_29 = 0.16666666666666666667/h1;
311            const double tmp0_40 = -0.62200846792814621559/h2;
312            const double tmp0_9 = 0.16666666666666666667/h0;
313            const double tmp0_30 = 0.62200846792814621559/h1;
314            const double tmp0_28 = -0.16666666666666666667/h1;
315            const double tmp0_43 = 0.044658198738520451079/h2;
316            const double tmp0_32 = 0.16666666666666666667/h1;
317            const double tmp0_31 = 0.044658198738520451079/h1;
318            const double tmp0_39 = 0.044658198738520451079/h1;
319            const double tmp0_58 = -0.62200846792814621559/h2;
320            const double tmp0_55 = 0.044658198738520451079/h2;
321            const double tmp0_18 = -0.62200846792814621559/h0;
322            const double tmp0_45 = -0.16666666666666666667/h2;
323            const double tmp0_59 = -0.16666666666666666667/h2;
324            const double tmp0_4 = -0.044658198738520451079/h0;
325            const double tmp0_19 = 0.044658198738520451079/h0;
326            const double tmp0_44 = -0.044658198738520451079/h2;
327            const double tmp0_5 = 0.62200846792814621559/h0;
328            const double tmp0_47 = 0.16666666666666666667/h2;
329            const double tmp0_6 = -0.16666666666666666667/h0;
330            const double tmp0_23 = 0.044658198738520451079/h1;
331            const double tmp0_46 = -0.16666666666666666667/h2;
332            const double tmp0_7 = -0.044658198738520451079/h0;
333    #pragma omp parallel for
334            for (index_t k2 =0; k2 < m_NE2; ++k2) {
335                for (index_t k1 =0; k1 < m_NE1; ++k1) {
336                    for (index_t k0 =0; k0 < m_NE0; ++k0) {
337                        const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
338                        const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
339                        const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
340                        const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
341                        const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
342                        const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
343                        const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
344                        const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
345                        double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
346                        for (index_t i=0; i < numComp; ++i) {
347                            o[INDEX3(i,0,0,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
348                            o[INDEX3(i,1,0,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);
349                            o[INDEX3(i,2,0,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);
350                            o[INDEX3(i,0,1,numComp,3)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_4 + f_100[i]*tmp0_5 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i]) + tmp0_2*(f_101[i] + f_110[i]);
351                            o[INDEX3(i,1,1,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
352                            o[INDEX3(i,2,1,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
353                            o[INDEX3(i,0,2,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
354                            o[INDEX3(i,1,2,numComp,3)] = f_000[i]*tmp0_20 + f_010[i]*tmp0_25 + f_101[i]*tmp0_22 + f_111[i]*tmp0_23 + tmp0_21*(f_001[i] + f_100[i]) + tmp0_24*(f_011[i] + f_110[i]);
355                            o[INDEX3(i,2,2,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
356                            o[INDEX3(i,0,3,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_7 + f_010[i]*tmp0_12 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_8 + f_110[i]*tmp0_10 + f_111[i]*tmp0_9;
357                            o[INDEX3(i,1,3,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_27 + f_010[i]*tmp0_32 + f_011[i]*tmp0_31 + f_100[i]*tmp0_33 + f_101[i]*tmp0_28 + f_110[i]*tmp0_30 + f_111[i]*tmp0_29;
358                            o[INDEX3(i,2,3,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);
359                            o[INDEX3(i,0,4,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
360                            o[INDEX3(i,1,4,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;
361                            o[INDEX3(i,2,4,numComp,3)] = f_000[i]*tmp0_40 + f_001[i]*tmp0_41 + f_110[i]*tmp0_44 + f_111[i]*tmp0_43 + tmp0_42*(f_011[i] + f_101[i]) + tmp0_45*(f_010[i] + f_100[i]);
362                            o[INDEX3(i,0,5,numComp,3)] = f_000[i]*tmp0_6 + f_001[i]*tmp0_12 + f_010[i]*tmp0_7 + f_011[i]*tmp0_11 + f_100[i]*tmp0_13 + f_101[i]*tmp0_10 + f_110[i]*tmp0_8 + f_111[i]*tmp0_9;
363                            o[INDEX3(i,1,5,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);
364                            o[INDEX3(i,2,5,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_52 + f_011[i]*tmp0_51 + f_100[i]*tmp0_53 + f_101[i]*tmp0_48 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
365                            o[INDEX3(i,0,6,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
366                            o[INDEX3(i,1,6,numComp,3)] = f_000[i]*tmp0_26 + f_001[i]*tmp0_33 + f_010[i]*tmp0_32 + f_011[i]*tmp0_30 + f_100[i]*tmp0_27 + f_101[i]*tmp0_28 + f_110[i]*tmp0_31 + f_111[i]*tmp0_29;
367                            o[INDEX3(i,2,6,numComp,3)] = f_000[i]*tmp0_46 + f_001[i]*tmp0_47 + f_010[i]*tmp0_53 + f_011[i]*tmp0_48 + f_100[i]*tmp0_52 + f_101[i]*tmp0_51 + f_110[i]*tmp0_50 + f_111[i]*tmp0_49;
368                            o[INDEX3(i,0,7,numComp,3)] = f_000[i]*tmp0_14 + f_011[i]*tmp0_18 + f_100[i]*tmp0_19 + f_111[i]*tmp0_17 + tmp0_15*(f_001[i] + f_010[i]) + tmp0_16*(f_101[i] + f_110[i]);
369                            o[INDEX3(i,1,7,numComp,3)] = f_000[i]*tmp0_34 + f_010[i]*tmp0_39 + f_101[i]*tmp0_36 + f_111[i]*tmp0_37 + tmp0_35*(f_001[i] + f_100[i]) + tmp0_38*(f_011[i] + f_110[i]);
370                            o[INDEX3(i,2,7,numComp,3)] = f_000[i]*tmp0_54 + f_001[i]*tmp0_55 + f_110[i]*tmp0_58 + f_111[i]*tmp0_57 + tmp0_56*(f_011[i] + f_101[i]) + tmp0_59*(f_010[i] + f_100[i]);
371                        } /* end of component loop i */
372                    } /* end of k0 loop */
373                } /* end of k1 loop */
374            } /* end of k2 loop */
375            /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
376        } else {
377            throw RipleyException("setToGradient() not implemented");
378        }
379    }
380    
381  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
382                                              bool reducedColOrder) const                                              bool reducedColOrder) const
383  {  {
# Line 554  void Brick::populateSampleIds() Line 671  void Brick::populateSampleIds()
671      for (dim_t k=0; k<getNumFaceElements(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
672          m_faceId[k]=k;          m_faceId[k]=k;
673      }      }
674    
675        // generate face offset vector
676        const IndexVector facesPerEdge = getNumFacesPerBoundary();
677        m_faceOffset.assign(facesPerEdge.size(), -1);
678        index_t offset=0;
679        for (size_t i=0; i<facesPerEdge.size(); i++) {
680            if (facesPerEdge[i]>0) {
681                m_faceOffset[i]=offset;
682                offset+=facesPerEdge[i];
683            }
684        }
685    }
686    
687    //protected
688    void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
689    {
690        const dim_t numComp = in.getDataPointSize();
691        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
692        const double tmp0_3 = 0.0094373878376559314545;
693        const double tmp0_2 = 0.035220810900864519624;
694        const double tmp0_1 = 0.13144585576580214704;
695        const double tmp0_0 = 0.49056261216234406855;
696    #pragma omp parallel for
697        for (index_t k2=0; k2 < m_NE2; ++k2) {
698            for (index_t k1=0; k1 < m_NE1; ++k1) {
699                for (index_t k0=0; k0 < m_NE0; ++k0) {
700                    const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
701                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
702                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
703                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
704                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
705                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
706                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
707                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
708                    double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
709                    for (index_t i=0; i < numComp; ++i) {
710                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_111[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_100[i]) + tmp0_2*(f_011[i] + f_101[i] + f_110[i]);
711                        o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_3 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_111[i]);
712                        o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_101[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_110[i]) + tmp0_2*(f_001[i] + f_100[i] + f_111[i]);
713                        o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_3 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_101[i]);
714                        o[INDEX2(i,numComp,4)] = f_001[i]*tmp0_0 + f_110[i]*tmp0_3 + tmp0_1*(f_000[i] + f_011[i] + f_101[i]) + tmp0_2*(f_010[i] + f_100[i] + f_111[i]);
715                        o[INDEX2(i,numComp,5)] = f_010[i]*tmp0_3 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i] + f_111[i]) + tmp0_2*(f_000[i] + f_011[i] + f_110[i]);
716                        o[INDEX2(i,numComp,6)] = f_011[i]*tmp0_0 + f_100[i]*tmp0_3 + tmp0_1*(f_001[i] + f_010[i] + f_111[i]) + tmp0_2*(f_000[i] + f_101[i] + f_110[i]);
717                        o[INDEX2(i,numComp,7)] = f_000[i]*tmp0_3 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i] + f_110[i]) + tmp0_2*(f_001[i] + f_010[i] + f_100[i]);
718                    } /* end of component loop i */
719                } /* end of k0 loop */
720            } /* end of k1 loop */
721        } /* end of k2 loop */
722        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
723    }
724    
725    //protected
726    void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
727    {
728        const dim_t numComp = in.getDataPointSize();
729        /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
730        if (m_faceOffset[0] > -1) {
731            const index_t k0 = 0;
732            const double tmp0_2 = 0.044658198738520451079;
733            const double tmp0_1 = 0.16666666666666666667;
734            const double tmp0_0 = 0.62200846792814621559;
735    #pragma omp parallel for
736            for (index_t k2=0; k2 < m_NE2; ++k2) {
737                for (index_t k1=0; k1 < m_NE1; ++k1) {
738                    const register double* f_000 = in.getSampleDataRO(INDEX3(0,k1,k2, m_N0,m_N1));
739                    const register double* f_001 = in.getSampleDataRO(INDEX3(0,k1,k2+1, m_N0,m_N1));
740                    const register double* f_011 = in.getSampleDataRO(INDEX3(0,k1+1,k2+1, m_N0,m_N1));
741                    const register double* f_010 = in.getSampleDataRO(INDEX3(0,k1+1,k2, m_N0,m_N1));
742                    double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
743                    for (index_t i=0; i < numComp; ++i) {
744                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_011[i]*tmp0_2 + tmp0_1*(f_001[i] + f_010[i]);
745                        o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_010[i]*tmp0_0 + tmp0_1*(f_000[i] + f_011[i]);
746                        o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_010[i]*tmp0_2 + tmp0_1*(f_000[i] + f_011[i]);
747                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_011[i]*tmp0_0 + tmp0_1*(f_001[i] + f_010[i]);
748                    } /* end of component loop i */
749                } /* end of k1 loop */
750            } /* end of k2 loop */
751        } /* end of face 0 */
752        if (m_faceOffset[1] > -1) {
753            const index_t k0 = 0;
754            const double tmp0_2 = 0.044658198738520451079;
755            const double tmp0_1 = 0.62200846792814621559;
756            const double tmp0_0 = 0.16666666666666666667;
757    #pragma omp parallel for
758            for (index_t k2=0; k2 < m_NE2; ++k2) {
759                for (index_t k1=0; k1 < m_NE1; ++k1) {
760                    const register double* f_101 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2+1, m_N0,m_N1));
761                    const register double* f_100 = in.getSampleDataRO(INDEX3(m_N0-1,k1,k2, m_N0,m_N1));
762                    const register double* f_110 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2, m_N0,m_N1));
763                    const register double* f_111 = in.getSampleDataRO(INDEX3(m_N0-1,k1+1,k2+1, m_N0,m_N1));
764                    double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
765                    for (index_t i=0; i < numComp; ++i) {
766                        o[INDEX2(i,numComp,0)] = f_100[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_101[i] + f_110[i]);
767                        o[INDEX2(i,numComp,1)] = f_101[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_100[i] + f_111[i]);
768                        o[INDEX2(i,numComp,2)] = f_101[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_100[i] + f_111[i]);
769                        o[INDEX2(i,numComp,3)] = f_100[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_101[i] + f_110[i]);
770                    } /* end of component loop i */
771                } /* end of k1 loop */
772            } /* end of k2 loop */
773        } /* end of face 1 */
774        if (m_faceOffset[2] > -1) {
775            const index_t k1 = 0;
776            const double tmp0_2 = 0.044658198738520451079;
777            const double tmp0_1 = 0.16666666666666666667;
778            const double tmp0_0 = 0.62200846792814621559;
779    #pragma omp parallel for
780            for (index_t k2=0; k2 < m_NE2; ++k2) {
781                for (index_t k0=0; k0 < m_NE0; ++k0) {
782                    const register double* f_000 = in.getSampleDataRO(INDEX3(k0,0,k2, m_N0,m_N1));
783                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,0,k2+1, m_N0,m_N1));
784                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,0,k2+1, m_N0,m_N1));
785                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,0,k2, m_N0,m_N1));
786                    double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
787                    for (index_t i=0; i < numComp; ++i) {
788                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_100[i]);
789                        o[INDEX2(i,numComp,1)] = f_001[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_101[i]);
790                        o[INDEX2(i,numComp,2)] = f_001[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_101[i]);
791                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_100[i]);
792                    } /* end of component loop i */
793                } /* end of k0 loop */
794            } /* end of k2 loop */
795        } /* end of face 2 */
796        if (m_faceOffset[3] > -1) {
797            const index_t k1 = 0;
798            const double tmp0_2 = 0.044658198738520451079;
799            const double tmp0_1 = 0.62200846792814621559;
800            const double tmp0_0 = 0.16666666666666666667;
801    #pragma omp parallel for
802            for (index_t k2=0; k2 < m_NE2; ++k2) {
803                for (index_t k0=0; k0 < m_NE0; ++k0) {
804                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2, m_N0,m_N1));
805                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2+1, m_N0,m_N1));
806                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,m_N1-1,k2, m_N0,m_N1));
807                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,m_N1-1,k2+1, m_N0,m_N1));
808                    double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
809                    for (index_t i=0; i < numComp; ++i) {
810                        o[INDEX2(i,numComp,0)] = f_010[i]*tmp0_1 + f_111[i]*tmp0_2 + tmp0_0*(f_011[i] + f_110[i]);
811                        o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_110[i]*tmp0_1 + tmp0_0*(f_010[i] + f_111[i]);
812                        o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_1 + f_110[i]*tmp0_2 + tmp0_0*(f_010[i] + f_111[i]);
813                        o[INDEX2(i,numComp,3)] = f_010[i]*tmp0_2 + f_111[i]*tmp0_1 + tmp0_0*(f_011[i] + f_110[i]);
814                    } /* end of component loop i */
815                } /* end of k0 loop */
816            } /* end of k2 loop */
817        } /* end of face 3 */
818        if (m_faceOffset[4] > -1) {
819            const index_t k2 = 0;
820            const double tmp0_2 = 0.044658198738520451079;
821            const double tmp0_1 = 0.16666666666666666667;
822            const double tmp0_0 = 0.62200846792814621559;
823    #pragma omp parallel for
824            for (index_t k1=0; k1 < m_NE1; ++k1) {
825                for (index_t k0=0; k0 < m_NE0; ++k0) {
826                    const register double* f_000 = in.getSampleDataRO(INDEX3(k0,k1,0, m_N0,m_N1));
827                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,0, m_N0,m_N1));
828                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,0, m_N0,m_N1));
829                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,0, m_N0,m_N1));
830                    double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
831                    for (index_t i=0; i < numComp; ++i) {
832                        o[INDEX2(i,numComp,0)] = f_000[i]*tmp0_0 + f_110[i]*tmp0_2 + tmp0_1*(f_010[i] + f_100[i]);
833                        o[INDEX2(i,numComp,1)] = f_010[i]*tmp0_2 + f_100[i]*tmp0_0 + tmp0_1*(f_000[i] + f_110[i]);
834                        o[INDEX2(i,numComp,2)] = f_010[i]*tmp0_0 + f_100[i]*tmp0_2 + tmp0_1*(f_000[i] + f_110[i]);
835                        o[INDEX2(i,numComp,3)] = f_000[i]*tmp0_2 + f_110[i]*tmp0_0 + tmp0_1*(f_010[i] + f_100[i]);
836                    } /* end of component loop i */
837                } /* end of k0 loop */
838            } /* end of k1 loop */
839        } /* end of face 4 */
840        if (m_faceOffset[5] > -1) {
841            const index_t k2 = 0;
842            const double tmp0_2 = 0.044658198738520451079;
843            const double tmp0_1 = 0.16666666666666666667;
844            const double tmp0_0 = 0.62200846792814621559;
845    #pragma omp parallel for
846            for (index_t k1=0; k1 < m_NE1; ++k1) {
847                for (index_t k0=0; k0 < m_NE0; ++k0) {
848                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,m_N2-1, m_N0,m_N1));
849                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,m_N2-1, m_N0,m_N1));
850                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,m_N2-1, m_N0,m_N1));
851                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,m_N2-1, m_N0,m_N1));
852                    double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX3(k0,k1,k2,m_NE0,m_NE1));
853                    for (index_t i=0; i < numComp; ++i) {
854                        o[INDEX2(i,numComp,0)] = f_001[i]*tmp0_0 + f_111[i]*tmp0_2 + tmp0_1*(f_011[i] + f_101[i]);
855                        o[INDEX2(i,numComp,1)] = f_011[i]*tmp0_2 + f_101[i]*tmp0_0 + tmp0_1*(f_001[i] + f_111[i]);
856                        o[INDEX2(i,numComp,2)] = f_011[i]*tmp0_0 + f_101[i]*tmp0_2 + tmp0_1*(f_001[i] + f_111[i]);
857                        o[INDEX2(i,numComp,3)] = f_001[i]*tmp0_2 + f_111[i]*tmp0_0 + tmp0_1*(f_011[i] + f_101[i]);
858                    } /* end of component loop i */
859                } /* end of k0 loop */
860            } /* end of k1 loop */
861        } /* end of face 5 */
862        /* GENERATOR SNIP_INTERPOLATE_FACES BOTTOM */
863  }  }
864    
865  } // end of namespace ripley  } // end of namespace ripley

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

  ViewVC Help
Powered by ViewVC 1.1.26