pyop2-laplaceknl.cpp

// PyOP2 kernel
static inline void form00_cell_integral_otherwise (double  A[3][3] , const double *const restrict *restrict coords )
{
    double  t19[4] ;
    double  t20[4] ;
    double  t0  = (-1 * coords[0][0]);
    double  t1  = (t0 + coords[1][0]);
    double  t2  = (-1 * coords[0][1]);
    double  t3  = (t2 + coords[2][1]);
    double  t4  = (t0 + coords[2][0]);
    double  t5  = (t2 + coords[1][1]);
    double  t6  = ((t1 * t3) + (-1 * (t4 * t5)));
    double  t7  = (0.5 * fabs(t6));
    double  t8  = (1 / t6);
    double  t9  = ((-1 * t5) * t8);
    double  t10  = (t3 * t8);
    double  t11  = (t1 * t8);
    double  t12  = ((-1 * t4) * t8);
    double  t13  = (t7 * ((t9 * t10) + (t11 * t12)));
    double  t14  = (t7 * ((t10 * t10) + (t12 * t12)));
    double  t15  = (t7 * ((t9 * t9) + (t11 * t11)));
    static const double  t16[4]  = {-1.0, 0.0, 1.0};
    double  t17  = (t7 * ((t10 * t9) + (t12 * t11)));
    static const double  t18[4]  = {-1.0, 1.0, 0.0};
    for (int  k  = 0; k < 3; k += 1)
    {
        t19[k] = (t18[k] * t17) + (t16[k] * t15);
        t20[k] = (t18[k] * t14) + (t16[k] * t13);
    }
    for (int  j  = 0; j < 3; j += 1)
    {
        for (int  k  = 0; k < 3; k += 1)
        {
            #pragma coffee expression
            A[j][k] += (t20[k] * t18[j]) + (t19[k] * t16[j]);
        }
    }
}
// We transalate this to opencl kernels
static inline void form_cell_integral_otherwise (double  A[1] , const double *const restrict *restrict coords , const double *const restrict *restric
t w_0 , const double *const restrict *restrict w_1 )
{
  static const double  t0[3][3]  = {{0.666666666666667, 0.166666666666667, 0.166666666666667},
  {0.166666666666667, 0.166666666666667, 0.666666666666667},
  {0.166666666666667, 0.666666666666667, 0.166666666666667}};
  double  t1  = (-1 * coords[0][0]);
  double  t2  = (-1 * coords[0][1]);
  double  t3  = fabs(((t1 + coords[1][0]) * (t2 + coords[2][1])) + (-1 * ((t1 + coords[2][0]) * (t2 + coords[1][1]))));
  static const double  t4[3]  = {0.166666666666667, 0.166666666666667, 0.166666666666667};
  for (int  ip  = 0; ip < 3; ip += 1)
  {
    #pragma coffee expression
    A[0] += (t4[ip] * t3) * pow((-1 * (((t0[ip][0] * w_1[0][0]) + (t0[ip][1] * w_1[1][0])) + (t0[ip][2] * w_1[2][0]))) + (((t0[ip][0] * w_0[0][0]) + (t0[ip][1] * w_0[1][0])) + (t0[ip][2] * w_0[2][0])), 2);
  }
}