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

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 3706 by caltinay, Mon Dec 5 01:59:08 2011 UTC revision 3707 by caltinay, Mon Dec 5 05:48:25 2011 UTC
# Line 279  void Rectangle::setToGradient(escript::D Line 279  void Rectangle::setToGradient(escript::D
279          const double tmp0_12 = -0.78867513459481288225/h1;          const double tmp0_12 = -0.78867513459481288225/h1;
280          const double tmp0_7 = -0.21132486540518711775/h0;          const double tmp0_7 = -0.21132486540518711775/h0;
281  #pragma omp parallel for  #pragma omp parallel for
282          for (index_t k1 =0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
283              for (index_t k0 =0; k0 < m_NE0; ++k0) {              for (index_t k0=0; k0 < m_NE0; ++k0) {
284                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,k1, m_N0));
285                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
286                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,k1+1, m_N0));
# Line 299  void Rectangle::setToGradient(escript::D Line 299  void Rectangle::setToGradient(escript::D
299              } /* end of k0 loop */              } /* end of k0 loop */
300          } /* end of k1 loop */          } /* end of k1 loop */
301          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */          /* GENERATOR SNIP_GRAD_ELEMENTS BOTTOM */
302        } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
303            /* GENERATOR SNIP_GRAD_FACES TOP */
304            if (m_faceOffset[0] > -1) {
305                const double tmp0_0 = 0.78867513459481288225/h0;
306                const double tmp0_4 = 0.21132486540518711775/h0;
307                const double tmp0_1 = 0.21132486540518711775/h0;
308                const double tmp0_8 = 1.0/h1;
309                const double tmp0_5 = 0.78867513459481288225/h0;
310                const double tmp0_2 = -0.21132486540518711775/h0;
311                const double tmp0_9 = -1/h1;
312                const double tmp0_6 = -0.78867513459481288225/h0;
313                const double tmp0_3 = -0.78867513459481288225/h0;
314                const double tmp0_7 = -0.21132486540518711775/h0;
315    #pragma omp parallel for
316                for (index_t k1=0; k1 < m_NE1; ++k1) {
317                    const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
318                    const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
319                    const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
320                    const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
321                    double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
322                    for (index_t i=0; i < numComp; ++i) {
323                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
324                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
325                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
326                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
327                    } /* end of component loop i */
328                } /* end of k1 loop */
329            } /* end of face 0 */
330            if (m_faceOffset[1] > -1) {
331                const double tmp0_0 = 0.78867513459481288225/h0;
332                const double tmp0_4 = 0.21132486540518711775/h0;
333                const double tmp0_1 = 0.21132486540518711775/h0;
334                const double tmp0_8 = -1/h1;
335                const double tmp0_5 = 0.78867513459481288225/h0;
336                const double tmp0_2 = -0.21132486540518711775/h0;
337                const double tmp0_9 = 1.0/h1;
338                const double tmp0_6 = -0.78867513459481288225/h0;
339                const double tmp0_3 = -0.78867513459481288225/h0;
340                const double tmp0_7 = -0.21132486540518711775/h0;
341    #pragma omp parallel for
342                for (index_t k1=0; k1 < m_NE1; ++k1) {
343                    const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
344                    const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
345                    const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
346                    const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
347                    double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
348                    for (index_t i=0; i < numComp; ++i) {
349                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2 + f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
350                        o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
351                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_7 + f_01[i]*tmp0_6 + f_10[i]*tmp0_4 + f_11[i]*tmp0_5;
352                        o[INDEX3(i,1,1,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
353                    } /* end of component loop i */
354                } /* end of k1 loop */
355            } /* end of face 1 */
356            if (m_faceOffset[2] > -1) {
357                const double tmp0_0 = -1/h0;
358                const double tmp0_4 = 0.21132486540518711775/h1;
359                const double tmp0_1 = 1.0/h0;
360                const double tmp0_8 = 0.78867513459481288225/h1;
361                const double tmp0_5 = 0.78867513459481288225/h1;
362                const double tmp0_2 = -0.78867513459481288225/h1;
363                const double tmp0_9 = 0.21132486540518711775/h1;
364                const double tmp0_6 = -0.21132486540518711775/h1;
365                const double tmp0_3 = -0.21132486540518711775/h1;
366                const double tmp0_7 = -0.78867513459481288225/h1;
367    #pragma omp parallel for
368                for (index_t k0=0; k0 < m_NE0; ++k0) {
369                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
370                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
371                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
372                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
373                    double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
374                    for (index_t i=0; i < numComp; ++i) {
375                        o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
376                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_2 + f_01[i]*tmp0_5 + f_10[i]*tmp0_3 + f_11[i]*tmp0_4;
377                        o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
378                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_6 + f_01[i]*tmp0_9 + f_10[i]*tmp0_7 + f_11[i]*tmp0_8;
379                    } /* end of component loop i */
380                } /* end of k0 loop */
381            } /* end of face 2 */
382            if (m_faceOffset[3] > -1) {
383                const double tmp0_0 = 1.0/h0;
384                const double tmp0_4 = -0.21132486540518711775/h1;
385                const double tmp0_1 = -1/h0;
386                const double tmp0_8 = -0.78867513459481288225/h1;
387                const double tmp0_5 = -0.78867513459481288225/h1;
388                const double tmp0_2 = 0.21132486540518711775/h1;
389                const double tmp0_9 = -0.21132486540518711775/h1;
390                const double tmp0_6 = 0.78867513459481288225/h1;
391                const double tmp0_3 = 0.78867513459481288225/h1;
392                const double tmp0_7 = 0.21132486540518711775/h1;
393    #pragma omp parallel for
394                for (index_t k0=0; k0 < m_NE0; ++k0) {
395                    const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
396                    const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
397                    const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
398                    const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
399                    double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
400                    for (index_t i=0; i < numComp; ++i) {
401                        o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
402                        o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_5 + f_01[i]*tmp0_3 + f_10[i]*tmp0_4 + f_11[i]*tmp0_2;
403                        o[INDEX3(i,0,1,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
404                        o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_7 + f_10[i]*tmp0_8 + f_11[i]*tmp0_6;
405                    } /* end of component loop i */
406                } /* end of k0 loop */
407            } /* end of face 3 */
408            /* GENERATOR SNIP_GRAD_FACES BOTTOM */
409      } else {      } else {
410          throw RipleyException("setToGradient() not implemented");          stringstream msg;
411            msg << "setToGradient() not implemented for "
412                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
413            throw RipleyException(msg.str());
414      }      }
415  }  }
416    
# Line 954  void Rectangle::interpolateNodesOnFaces( Line 1064  void Rectangle::interpolateNodesOnFaces(
1064      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1065      /* GENERATOR SNIP_INTERPOLATE_FACES TOP */      /* GENERATOR SNIP_INTERPOLATE_FACES TOP */
1066      if (m_faceOffset[0] > -1) {      if (m_faceOffset[0] > -1) {
         const index_t k0 = 0;  
1067          const double tmp0_1 = 0.78867513459481288225;          const double tmp0_1 = 0.78867513459481288225;
1068          const double tmp0_0 = 0.21132486540518711775;          const double tmp0_0 = 0.21132486540518711775;
1069  #pragma omp parallel for  #pragma omp parallel for
1070          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1071              const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));              const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
1072              const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));              const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
1073              double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k0,k1,m_NE0));              double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
1074              for (index_t i=0; i < numComp; ++i) {              for (index_t i=0; i < numComp; ++i) {
1075                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_01[i]*tmp0_0;                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_01[i]*tmp0_0;
1076                  o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_01[i]*tmp0_1;                  o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_01[i]*tmp0_1;
# Line 969  void Rectangle::interpolateNodesOnFaces( Line 1078  void Rectangle::interpolateNodesOnFaces(
1078          } /* end of k1 loop */          } /* end of k1 loop */
1079      } /* end of face 0 */      } /* end of face 0 */
1080      if (m_faceOffset[1] > -1) {      if (m_faceOffset[1] > -1) {
         const index_t k0 = 0;  
1081          const double tmp0_1 = 0.21132486540518711775;          const double tmp0_1 = 0.21132486540518711775;
1082          const double tmp0_0 = 0.78867513459481288225;          const double tmp0_0 = 0.78867513459481288225;
1083  #pragma omp parallel for  #pragma omp parallel for
1084          for (index_t k1=0; k1 < m_NE1; ++k1) {          for (index_t k1=0; k1 < m_NE1; ++k1) {
1085              const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));              const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
1086              const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));              const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
1087              double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k0,k1,m_NE0));              double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
1088              for (index_t i=0; i < numComp; ++i) {              for (index_t i=0; i < numComp; ++i) {
1089                  o[INDEX2(i,numComp,0)] = f_10[i]*tmp0_0 + f_11[i]*tmp0_1;                  o[INDEX2(i,numComp,0)] = f_10[i]*tmp0_0 + f_11[i]*tmp0_1;
1090                  o[INDEX2(i,numComp,1)] = f_10[i]*tmp0_1 + f_11[i]*tmp0_0;                  o[INDEX2(i,numComp,1)] = f_10[i]*tmp0_1 + f_11[i]*tmp0_0;
# Line 984  void Rectangle::interpolateNodesOnFaces( Line 1092  void Rectangle::interpolateNodesOnFaces(
1092          } /* end of k1 loop */          } /* end of k1 loop */
1093      } /* end of face 1 */      } /* end of face 1 */
1094      if (m_faceOffset[2] > -1) {      if (m_faceOffset[2] > -1) {
         const index_t k1 = 0;  
1095          const double tmp0_1 = 0.78867513459481288225;          const double tmp0_1 = 0.78867513459481288225;
1096          const double tmp0_0 = 0.21132486540518711775;          const double tmp0_0 = 0.21132486540518711775;
1097  #pragma omp parallel for  #pragma omp parallel for
1098          for (index_t k0=0; k0 < m_NE0; ++k0) {          for (index_t k0=0; k0 < m_NE0; ++k0) {
1099              const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));              const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
1100              const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));              const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
1101              double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k1,m_NE0));              double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
1102              for (index_t i=0; i < numComp; ++i) {              for (index_t i=0; i < numComp; ++i) {
1103                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_10[i]*tmp0_0;                  o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_1 + f_10[i]*tmp0_0;
1104                  o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                  o[INDEX2(i,numComp,1)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
# Line 999  void Rectangle::interpolateNodesOnFaces( Line 1106  void Rectangle::interpolateNodesOnFaces(
1106          } /* end of k0 loop */          } /* end of k0 loop */
1107      } /* end of face 2 */      } /* end of face 2 */
1108      if (m_faceOffset[3] > -1) {      if (m_faceOffset[3] > -1) {
         const index_t k1 = 0;  
1109          const double tmp0_1 = 0.78867513459481288225;          const double tmp0_1 = 0.78867513459481288225;
1110          const double tmp0_0 = 0.21132486540518711775;          const double tmp0_0 = 0.21132486540518711775;
1111  #pragma omp parallel for  #pragma omp parallel for
1112          for (index_t k0=0; k0 < m_NE0; ++k0) {          for (index_t k0=0; k0 < m_NE0; ++k0) {
1113              const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));              const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
1114              const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));              const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
1115              double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k1,m_NE0));              double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
1116              for (index_t i=0; i < numComp; ++i) {              for (index_t i=0; i < numComp; ++i) {
1117                  o[INDEX2(i,numComp,0)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                  o[INDEX2(i,numComp,0)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
1118                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_0 + f_11[i]*tmp0_1;                  o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_0 + f_11[i]*tmp0_1;

Legend:
Removed from v.3706  
changed lines
  Added in v.3707

  ViewVC Help
Powered by ViewVC 1.1.26