/[escript]/trunk/ripley/src/Brick.cpp
ViewVC logotype

Diff of /trunk/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 3703 by caltinay, Sun Dec 4 23:42:52 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 556  void Brick::populateSampleIds() Line 673  void Brick::populateSampleIds()
673      }      }
674  }  }
675    
676    //protected
677    void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
678    {
679        const dim_t numComp = in.getDataPointSize();
680        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS TOP */
681        const double tmp0_3 = 0.0094373878376559314545;
682        const double tmp0_2 = 0.035220810900864519624;
683        const double tmp0_1 = 0.13144585576580214704;
684        const double tmp0_0 = 0.49056261216234406855;
685    #pragma omp parallel for
686        for (index_t k2=0; k2 < m_NE2; ++k2) {
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_000 = in.getSampleDataRO(INDEX3(k0,k1,k2, m_N0,m_N1));
690                    const register double* f_001 = in.getSampleDataRO(INDEX3(k0,k1,k2+1, m_N0,m_N1));
691                    const register double* f_101 = in.getSampleDataRO(INDEX3(k0+1,k1,k2+1, m_N0,m_N1));
692                    const register double* f_011 = in.getSampleDataRO(INDEX3(k0,k1+1,k2+1, m_N0,m_N1));
693                    const register double* f_110 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2, m_N0,m_N1));
694                    const register double* f_010 = in.getSampleDataRO(INDEX3(k0,k1+1,k2, m_N0,m_N1));
695                    const register double* f_100 = in.getSampleDataRO(INDEX3(k0+1,k1,k2, m_N0,m_N1));
696                    const register double* f_111 = in.getSampleDataRO(INDEX3(k0+1,k1+1,k2+1, m_N0,m_N1));
697                    double* o = out.getSampleDataRW(INDEX3(k0,k1,k2,m_NE0,m_NE1));
698                    for (index_t i=0; i < numComp; ++i) {
699                        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]);
700                        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]);
701                        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]);
702                        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]);
703                        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]);
704                        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]);
705                        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]);
706                        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]);
707                    } /* end of component loop i */
708                } /* end of k0 loop */
709            } /* end of k1 loop */
710        } /* end of k2 loop */
711        /* GENERATOR SNIP_INTERPOLATE_ELEMENTS BOTTOM */
712    }
713    
714    //protected
715    void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
716    {
717        throw RipleyException("interpolateNodesOnFaces() not implemented");
718    }
719    
720  } // end of namespace ripley  } // end of namespace ripley
721    

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

  ViewVC Help
Powered by ViewVC 1.1.26