/[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 3713 by caltinay, Tue Dec 6 04:43:29 2011 UTC revision 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC
# Line 142  void Rectangle::dump(const string& fileN Line 142  void Rectangle::dump(const string& fileN
142      pair<double,double> ydy = getFirstCoordAndSpacing(1);      pair<double,double> ydy = getFirstCoordAndSpacing(1);
143  #pragma omp parallel  #pragma omp parallel
144      {      {
145  #pragma omp for  #pragma omp for nowait
146          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
147              coords[0][i0]=xdx.first+i0*xdx.second;              coords[0][i0]=xdx.first+i0*xdx.second;
148          }          }
149  #pragma omp for  #pragma omp for nowait
150          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
151              coords[1][i1]=ydy.first+i1*ydy.second;              coords[1][i1]=ydy.first+i1*ydy.second;
152          }          }
# Line 321  void Rectangle::setToGradient(escript::D Line 321  void Rectangle::setToGradient(escript::D
321          } /* end of k1 loop */          } /* end of k1 loop */
322          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */          /* GENERATOR SNIP_GRAD_REDUCED_ELEMENTS BOTTOM */
323      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
324          /* GENERATOR SNIP_GRAD_FACES TOP */  #pragma omp parallel
325          if (m_faceOffset[0] > -1) {          {
326              const double tmp0_0 = 0.78867513459481288225/h0;              /* GENERATOR SNIP_GRAD_FACES TOP */
327              const double tmp0_4 = 0.21132486540518711775/h0;              if (m_faceOffset[0] > -1) {
328              const double tmp0_1 = 0.21132486540518711775/h0;                  const double tmp0_0 = 0.78867513459481288225/h0;
329              const double tmp0_8 = 1.0/h1;                  const double tmp0_4 = 0.21132486540518711775/h0;
330              const double tmp0_5 = 0.78867513459481288225/h0;                  const double tmp0_1 = 0.21132486540518711775/h0;
331              const double tmp0_2 = -0.21132486540518711775/h0;                  const double tmp0_8 = 1.0/h1;
332              const double tmp0_9 = -1/h1;                  const double tmp0_5 = 0.78867513459481288225/h0;
333              const double tmp0_6 = -0.78867513459481288225/h0;                  const double tmp0_2 = -0.21132486540518711775/h0;
334              const double tmp0_3 = -0.78867513459481288225/h0;                  const double tmp0_9 = -1/h1;
335              const double tmp0_7 = -0.21132486540518711775/h0;                  const double tmp0_6 = -0.78867513459481288225/h0;
336  #pragma omp parallel for                  const double tmp0_3 = -0.78867513459481288225/h0;
337              for (index_t k1=0; k1 < m_NE1; ++k1) {                  const double tmp0_7 = -0.21132486540518711775/h0;
338                  const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  #pragma omp for nowait
339                  const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
340                  const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
341                  const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
342                  double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
343                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
344                      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;                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
345                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;                      for (index_t i=0; i < numComp; ++i) {
346                      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;                          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;
347                      o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
348                  } /* end of component loop i */                          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;
349              } /* end of k1 loop */                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*tmp0_9 + f_01[i]*tmp0_8;
350          } /* end of face 0 */                      } /* end of component loop i */
351          if (m_faceOffset[1] > -1) {                  } /* end of k1 loop */
352              const double tmp0_0 = 0.78867513459481288225/h0;              } /* end of face 0 */
353              const double tmp0_4 = 0.21132486540518711775/h0;              if (m_faceOffset[1] > -1) {
354              const double tmp0_1 = 0.21132486540518711775/h0;                  const double tmp0_0 = 0.78867513459481288225/h0;
355              const double tmp0_8 = -1/h1;                  const double tmp0_4 = 0.21132486540518711775/h0;
356              const double tmp0_5 = 0.78867513459481288225/h0;                  const double tmp0_1 = 0.21132486540518711775/h0;
357              const double tmp0_2 = -0.21132486540518711775/h0;                  const double tmp0_8 = -1/h1;
358              const double tmp0_9 = 1.0/h1;                  const double tmp0_5 = 0.78867513459481288225/h0;
359              const double tmp0_6 = -0.78867513459481288225/h0;                  const double tmp0_2 = -0.21132486540518711775/h0;
360              const double tmp0_3 = -0.78867513459481288225/h0;                  const double tmp0_9 = 1.0/h1;
361              const double tmp0_7 = -0.21132486540518711775/h0;                  const double tmp0_6 = -0.78867513459481288225/h0;
362  #pragma omp parallel for                  const double tmp0_3 = -0.78867513459481288225/h0;
363              for (index_t k1=0; k1 < m_NE1; ++k1) {                  const double tmp0_7 = -0.21132486540518711775/h0;
364                  const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  #pragma omp for nowait
365                  const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
366                  const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
367                  const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
368                  double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
369                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
370                      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;                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
371                      o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                      for (index_t i=0; i < numComp; ++i) {
372                      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;                          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;
373                      o[INDEX3(i,1,1,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
374                  } /* end of component loop i */                          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;
375              } /* end of k1 loop */                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*tmp0_8 + f_11[i]*tmp0_9;
376          } /* end of face 1 */                      } /* end of component loop i */
377          if (m_faceOffset[2] > -1) {                  } /* end of k1 loop */
378              const double tmp0_0 = -1/h0;              } /* end of face 1 */
379              const double tmp0_4 = 0.21132486540518711775/h1;              if (m_faceOffset[2] > -1) {
380              const double tmp0_1 = 1.0/h0;                  const double tmp0_0 = -1/h0;
381              const double tmp0_8 = 0.78867513459481288225/h1;                  const double tmp0_4 = 0.21132486540518711775/h1;
382              const double tmp0_5 = 0.78867513459481288225/h1;                  const double tmp0_1 = 1.0/h0;
383              const double tmp0_2 = -0.78867513459481288225/h1;                  const double tmp0_8 = 0.78867513459481288225/h1;
384              const double tmp0_9 = 0.21132486540518711775/h1;                  const double tmp0_5 = 0.78867513459481288225/h1;
385              const double tmp0_6 = -0.21132486540518711775/h1;                  const double tmp0_2 = -0.78867513459481288225/h1;
386              const double tmp0_3 = -0.21132486540518711775/h1;                  const double tmp0_9 = 0.21132486540518711775/h1;
387              const double tmp0_7 = -0.78867513459481288225/h1;                  const double tmp0_6 = -0.21132486540518711775/h1;
388  #pragma omp parallel for                  const double tmp0_3 = -0.21132486540518711775/h1;
389              for (index_t k0=0; k0 < m_NE0; ++k0) {                  const double tmp0_7 = -0.78867513459481288225/h1;
390                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  #pragma omp for nowait
391                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
392                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
393                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
394                  double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
395                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
396                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
397                      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;                      for (index_t i=0; i < numComp; ++i) {
398                      o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
399                      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;                          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;
400                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
401              } /* end of k0 loop */                          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;
402          } /* end of face 2 */                      } /* end of component loop i */
403          if (m_faceOffset[3] > -1) {                  } /* end of k0 loop */
404              const double tmp0_0 = 1.0/h0;              } /* end of face 2 */
405              const double tmp0_4 = -0.21132486540518711775/h1;              if (m_faceOffset[3] > -1) {
406              const double tmp0_1 = -1/h0;                  const double tmp0_0 = 1.0/h0;
407              const double tmp0_8 = -0.78867513459481288225/h1;                  const double tmp0_4 = -0.21132486540518711775/h1;
408              const double tmp0_5 = -0.78867513459481288225/h1;                  const double tmp0_1 = -1/h0;
409              const double tmp0_2 = 0.21132486540518711775/h1;                  const double tmp0_8 = -0.78867513459481288225/h1;
410              const double tmp0_9 = -0.21132486540518711775/h1;                  const double tmp0_5 = -0.78867513459481288225/h1;
411              const double tmp0_6 = 0.78867513459481288225/h1;                  const double tmp0_2 = 0.21132486540518711775/h1;
412              const double tmp0_3 = 0.78867513459481288225/h1;                  const double tmp0_9 = -0.21132486540518711775/h1;
413              const double tmp0_7 = 0.21132486540518711775/h1;                  const double tmp0_6 = 0.78867513459481288225/h1;
414  #pragma omp parallel for                  const double tmp0_3 = 0.78867513459481288225/h1;
415              for (index_t k0=0; k0 < m_NE0; ++k0) {                  const double tmp0_7 = 0.21132486540518711775/h1;
416                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  #pragma omp for nowait
417                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
418                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
419                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
420                  double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
421                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
422                      o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
423                      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;                      for (index_t i=0; i < numComp; ++i) {
424                      o[INDEX3(i,0,1,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
425                      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;                          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;
426                  } /* end of component loop i */                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
427              } /* end of k0 loop */                          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;
428          } /* end of face 3 */                      } /* end of component loop i */
429          /* GENERATOR SNIP_GRAD_FACES BOTTOM */                  } /* end of k0 loop */
430                } /* end of face 3 */
431                /* GENERATOR SNIP_GRAD_FACES BOTTOM */
432            } // end of parallel section
433      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
434          /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */  #pragma omp parallel
435          if (m_faceOffset[0] > -1) {          {
436              const double tmp0_3 = -1/h1;              /* GENERATOR SNIP_GRAD_REDUCED_FACES TOP */
437              const double tmp0_2 = 1.0/h1;              if (m_faceOffset[0] > -1) {
438              const double tmp0_1 = -0.5/h0;                  const double tmp0_3 = -1/h1;
439              const double tmp0_0 = 0.5/h0;                  const double tmp0_2 = 1.0/h1;
440  #pragma omp parallel for                  const double tmp0_1 = -0.5/h0;
441              for (index_t k1=0; k1 < m_NE1; ++k1) {                  const double tmp0_0 = 0.5/h0;
442                  const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));  #pragma omp for nowait
443                  const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
444                  const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(1,k1, m_N0));
445                  const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(1,k1+1, m_N0));
446                  double* o = out.getSampleDataRW(m_faceOffset[0]+k1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(0,k1+1, m_N0));
447                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(0,k1, m_N0));
448                      o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);                      double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
449                      o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2;                      for (index_t i=0; i < numComp; ++i) {
450                  } /* end of component loop i */                          o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);
451              } /* end of k1 loop */                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*tmp0_3 + f_01[i]*tmp0_2;
452          } /* end of face 0 */                      } /* end of component loop i */
453          if (m_faceOffset[1] > -1) {                  } /* end of k1 loop */
454              const double tmp0_3 = 1.0/h1;              } /* end of face 0 */
455              const double tmp0_2 = -1/h1;              if (m_faceOffset[1] > -1) {
456              const double tmp0_1 = -0.5/h0;                  const double tmp0_3 = 1.0/h1;
457              const double tmp0_0 = 0.5/h0;                  const double tmp0_2 = -1/h1;
458  #pragma omp parallel for                  const double tmp0_1 = -0.5/h0;
459              for (index_t k1=0; k1 < m_NE1; ++k1) {                  const double tmp0_0 = 0.5/h0;
460                  const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));  #pragma omp for nowait
461                  const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));                  for (index_t k1=0; k1 < m_NE1; ++k1) {
462                  const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(m_N0-1,k1, m_N0));
463                  const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(m_N0-1,k1+1, m_N0));
464                  double* o = out.getSampleDataRW(m_faceOffset[1]+k1);                      const register double* f_01 = in.getSampleDataRO(INDEX2(m_N0-2,k1+1, m_N0));
465                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(m_N0-2,k1, m_N0));
466                      o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);                      double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
467                      o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_2 + f_11[i]*tmp0_3;                      for (index_t i=0; i < numComp; ++i) {
468                  } /* end of component loop i */                          o[INDEX3(i,0,0,numComp,2)] = tmp0_0*(f_10[i] + f_11[i]) + tmp0_1*(f_00[i] + f_01[i]);
469              } /* end of k1 loop */                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*tmp0_2 + f_11[i]*tmp0_3;
470          } /* end of face 1 */                      } /* end of component loop i */
471          if (m_faceOffset[2] > -1) {                  } /* end of k1 loop */
472              const double tmp0_3 = 0.5/h1;              } /* end of face 1 */
473              const double tmp0_2 = -0.5/h1;              if (m_faceOffset[2] > -1) {
474              const double tmp0_1 = 1.0/h0;                  const double tmp0_3 = 0.5/h1;
475              const double tmp0_0 = -1/h0;                  const double tmp0_2 = -0.5/h1;
476  #pragma omp parallel for                  const double tmp0_1 = 1.0/h0;
477              for (index_t k0=0; k0 < m_NE0; ++k0) {                  const double tmp0_0 = -1/h0;
478                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));  #pragma omp for nowait
479                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
480                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,0, m_N0));
481                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,0, m_N0));
482                  double* o = out.getSampleDataRW(m_faceOffset[2]+k0);                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,1, m_N0));
483                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,1, m_N0));
484                      o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;                      double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
485                      o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_00[i] + f_10[i]) + tmp0_3*(f_01[i] + f_11[i]);                      for (index_t i=0; i < numComp; ++i) {
486                  } /* end of component loop i */                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*tmp0_0 + f_10[i]*tmp0_1;
487              } /* end of k0 loop */                          o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_00[i] + f_10[i]) + tmp0_3*(f_01[i] + f_11[i]);
488          } /* end of face 2 */                      } /* end of component loop i */
489          if (m_faceOffset[3] > -1) {                  } /* end of k0 loop */
490              const double tmp0_3 = -0.5/h1;              } /* end of face 2 */
491              const double tmp0_2 = 0.5/h1;              if (m_faceOffset[3] > -1) {
492              const double tmp0_1 = -1/h0;                  const double tmp0_3 = -0.5/h1;
493              const double tmp0_0 = 1.0/h0;                  const double tmp0_2 = 0.5/h1;
494  #pragma omp parallel for                  const double tmp0_1 = -1/h0;
495              for (index_t k0=0; k0 < m_NE0; ++k0) {                  const double tmp0_0 = 1.0/h0;
496                  const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));  #pragma omp for nowait
497                  const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));                  for (index_t k0=0; k0 < m_NE0; ++k0) {
498                  const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));                      const register double* f_11 = in.getSampleDataRO(INDEX2(k0+1,m_N1-1, m_N0));
499                  const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));                      const register double* f_01 = in.getSampleDataRO(INDEX2(k0,m_N1-1, m_N0));
500                  double* o = out.getSampleDataRW(m_faceOffset[3]+k0);                      const register double* f_10 = in.getSampleDataRO(INDEX2(k0+1,m_N1-2, m_N0));
501                  for (index_t i=0; i < numComp; ++i) {                      const register double* f_00 = in.getSampleDataRO(INDEX2(k0,m_N1-2, m_N0));
502                      o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;                      double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
503                      o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_01[i] + f_11[i]) + tmp0_3*(f_00[i] + f_10[i]);                      for (index_t i=0; i < numComp; ++i) {
504                  } /* end of component loop i */                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*tmp0_1 + f_11[i]*tmp0_0;
505              } /* end of k0 loop */                          o[INDEX3(i,1,0,numComp,2)] = tmp0_2*(f_01[i] + f_11[i]) + tmp0_3*(f_00[i] + f_10[i]);
506          } /* end of face 3 */                      } /* end of component loop i */
507          /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */                  } /* end of k0 loop */
508                } /* end of face 3 */
509                /* GENERATOR SNIP_GRAD_REDUCED_FACES BOTTOM */
510            } // end of parallel section
511      } else {      } else {
512          stringstream msg;          stringstream msg;
513          msg << "setToGradient() not implemented for "          msg << "setToGradient() not implemented for "
# Line 521  void Rectangle::setToIntegrals(vector<do Line 527  void Rectangle::setToIntegrals(vector<do
527  #pragma omp parallel  #pragma omp parallel
528          {          {
529              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
530  #pragma omp for  #pragma omp for nowait
531              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = 0; k1 < m_NE1; ++k1) {
532                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
533                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
# Line 538  void Rectangle::setToIntegrals(vector<do Line 544  void Rectangle::setToIntegrals(vector<do
544  #pragma omp critical  #pragma omp critical
545              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
546                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
547          }          } // end of parallel section
548      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
549          const double w_0 = h0*h1;          const double w_0 = h0*h1;
550  #pragma omp parallel  #pragma omp parallel
551          {          {
552              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
553  #pragma omp for  #pragma omp for nowait
554              for (index_t k1 = 0; k1 < m_NE1; ++k1) {              for (index_t k1 = 0; k1 < m_NE1; ++k1) {
555                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
556                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
# Line 557  void Rectangle::setToIntegrals(vector<do Line 563  void Rectangle::setToIntegrals(vector<do
563  #pragma omp critical  #pragma omp critical
564              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
565                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
566          }          } // end of parallel section
567      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
568          const double w_0 = h0/2.;          const double w_0 = h0/2.;
569          const double w_1 = h1/2.;          const double w_1 = h1/2.;
# Line 565  void Rectangle::setToIntegrals(vector<do Line 571  void Rectangle::setToIntegrals(vector<do
571          {          {
572              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
573              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
574  #pragma omp for  #pragma omp for nowait
575                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
576                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
577                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 577  void Rectangle::setToIntegrals(vector<do Line 583  void Rectangle::setToIntegrals(vector<do
583              }              }
584    
585              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
586  #pragma omp for  #pragma omp for nowait
587                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
588                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
589                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 589  void Rectangle::setToIntegrals(vector<do Line 595  void Rectangle::setToIntegrals(vector<do
595              }              }
596    
597              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
598  #pragma omp for  #pragma omp for nowait
599                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
600                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
601                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 601  void Rectangle::setToIntegrals(vector<do Line 607  void Rectangle::setToIntegrals(vector<do
607              }              }
608    
609              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
610  #pragma omp for  #pragma omp for nowait
611                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
612                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
613                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 615  void Rectangle::setToIntegrals(vector<do Line 621  void Rectangle::setToIntegrals(vector<do
621  #pragma omp critical  #pragma omp critical
622              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
623                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
624          }          } // end of parallel section
625      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
626  #pragma omp parallel  #pragma omp parallel
627          {          {
628              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
629              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
630  #pragma omp for  #pragma omp for nowait
631                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
632                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
633                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 631  void Rectangle::setToIntegrals(vector<do Line 637  void Rectangle::setToIntegrals(vector<do
637              }              }
638    
639              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
640  #pragma omp for  #pragma omp for nowait
641                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {                  for (index_t k1 = 0; k1 < m_NE1; ++k1) {
642                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
643                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 641  void Rectangle::setToIntegrals(vector<do Line 647  void Rectangle::setToIntegrals(vector<do
647              }              }
648    
649              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
650  #pragma omp for  #pragma omp for nowait
651                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
652                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
653                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 651  void Rectangle::setToIntegrals(vector<do Line 657  void Rectangle::setToIntegrals(vector<do
657              }              }
658    
659              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
660  #pragma omp for  #pragma omp for nowait
661                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {                  for (index_t k0 = 0; k0 < m_NE0; ++k0) {
662                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
663                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
# Line 663  void Rectangle::setToIntegrals(vector<do Line 669  void Rectangle::setToIntegrals(vector<do
669  #pragma omp critical  #pragma omp critical
670              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
671                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
672          }          } // end of parallel section
673      } else {      } else {
674          stringstream msg;          stringstream msg;
675          msg << "setToIntegrals() not implemented for "          msg << "setToIntegrals() not implemented for "
# Line 672  void Rectangle::setToIntegrals(vector<do Line 678  void Rectangle::setToIntegrals(vector<do
678      }      }
679  }  }
680    
681    void Rectangle::setToNormal(escript::Data& out) const
682    {
683        if (out.getFunctionSpace().getTypeCode() == FaceElements) {
684    #pragma omp parallel
685            {
686                if (m_faceOffset[0] > -1) {
687    #pragma omp for nowait
688                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
689                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
690                        // set vector at two quadrature points
691                        *o++ = -1.;
692                        *o++ = 0.;
693                        *o++ = -1.;
694                        *o = 0.;
695                    }
696                }
697    
698                if (m_faceOffset[1] > -1) {
699    #pragma omp for nowait
700                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
701                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
702                        // set vector at two quadrature points
703                        *o++ = 1.;
704                        *o++ = 0.;
705                        *o++ = 1.;
706                        *o = 0.;
707                    }
708                }
709    
710                if (m_faceOffset[2] > -1) {
711    #pragma omp for nowait
712                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
713                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
714                        // set vector at two quadrature points
715                        *o++ = 0.;
716                        *o++ = -1.;
717                        *o++ = 0.;
718                        *o = -1.;
719                    }
720                }
721    
722                if (m_faceOffset[3] > -1) {
723    #pragma omp for nowait
724                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
725                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
726                        // set vector at two quadrature points
727                        *o++ = 0.;
728                        *o++ = 1.;
729                        *o++ = 0.;
730                        *o = 1.;
731                    }
732                }
733            } // end of parallel section
734        } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
735    #pragma omp parallel
736            {
737                if (m_faceOffset[0] > -1) {
738    #pragma omp for nowait
739                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
740                        double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
741                        *o++ = -1.;
742                        *o = 0.;
743                    }
744                }
745    
746                if (m_faceOffset[1] > -1) {
747    #pragma omp for nowait
748                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
749                        double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
750                        *o++ = 1.;
751                        *o = 0.;
752                    }
753                }
754    
755                if (m_faceOffset[2] > -1) {
756    #pragma omp for nowait
757                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
758                        double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
759                        *o++ = 0.;
760                        *o = -1.;
761                    }
762                }
763    
764                if (m_faceOffset[3] > -1) {
765    #pragma omp for nowait
766                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
767                        double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
768                        *o++ = 0.;
769                        *o = 1.;
770                    }
771                }
772            } // end of parallel section
773    
774        } else {
775            stringstream msg;
776            msg << "setToNormal() not implemented for "
777                << functionSpaceTypeAsString(out.getFunctionSpace().getTypeCode());
778            throw RipleyException(msg.str());
779        }
780    }
781    
782  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,
783          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
784          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,
# Line 1027  void Rectangle::populateSampleIds() Line 1134  void Rectangle::populateSampleIds()
1134          m_faceId[k]=k;          m_faceId[k]=k;
1135      }      }
1136    
1137      // generate face offset vector      // generate face offset vector and set face tags
1138      const IndexVector facesPerEdge = getNumFacesPerBoundary();      const IndexVector facesPerEdge = getNumFacesPerBoundary();
1139        const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
1140        const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
1141      m_faceOffset.assign(facesPerEdge.size(), -1);      m_faceOffset.assign(facesPerEdge.size(), -1);
1142        m_faceTags.clear();
1143      index_t offset=0;      index_t offset=0;
1144      for (size_t i=0; i<facesPerEdge.size(); i++) {      for (size_t i=0; i<facesPerEdge.size(); i++) {
1145          if (facesPerEdge[i]>0) {          if (facesPerEdge[i]>0) {
1146              m_faceOffset[i]=offset;              m_faceOffset[i]=offset;
1147              offset+=facesPerEdge[i];              offset+=facesPerEdge[i];
1148                m_faceTags.insert(m_faceTags.end(), facesPerEdge[i], faceTag[i]);
1149          }          }
1150      }      }
1151        setTagMap("left", LEFT);
1152        setTagMap("right", RIGHT);
1153        setTagMap("bottom", BOTTOM);
1154        setTagMap("top", TOP);
1155        updateTagsInUse(FaceElements);
1156  }  }
1157    
1158  //private  //private

Legend:
Removed from v.3713  
changed lines
  Added in v.3722

  ViewVC Help
Powered by ViewVC 1.1.26