double approx::rhs_elem (int k)
{

  int i, j;
  int n = msh->n_mesh;
  int res = 0;
  k2ij (n, i, j, k);
  bool a = 0, b = 0, c = 0, d = 0, e = 0, f = 0;
  double D = msh->box_d;
  double l = msh->box_l;
  double dx = msh->dx;
  double dy = msh->dy;

  if (! (msh->ptempl[j * (n + 1) + i]))
    {
      return 0;
    }

  // +0, +1
  if (j == msh->n_mesh) res += 0;
  else if (msh->ptempl[(j + 1) * (n + 1) + i])
    {
      a = 1;
    }

  // +1, +1
  if (j == n) res += 0;
  else if (i == n) res += 0;
  else if (msh->ptempl[(j + 1) * (n + 1) + i + 1])
    {
      b = 1;
    }

  // +1, +0
  if (i == n) res += 0;
  else if (msh->ptempl[j * (n + 1) + i + 1])
    {
      c = 1;
    }

  // +0, -1
  if (j == 0) res += 0;
  else if (msh->ptempl[(j - 1) * (n + 1) + i])
    {
      d = 1;
    }

  // -1, -1
  if (j == 0) res += 0;
  else if (i == 0) res += 0;
  else if (msh->ptempl[(j - 1) * (n + 1) + i - 1])
    {
      e = 1;
    }

  //-1, 0
  if (i == 0) res += 0;
  else if (msh->ptempl[j * (n + 1) + i - 1])
    {
      f = 1;
    }


  double r = 0;

  if (a && b)
    {
      r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
      r += 1./12. * 1./4. * 1./2. * (F (l + i * dx, D + j * dy + dy/2) + F (l + i * dx + dx/2, D + j * dy + dy/2));
    }

  if (c && b)
    {
      r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
      r += 1./12. * 1./4. * 1./2. * (F (l + i * dx + dx/2, D + j * dy) + F (l + i * dx + dx/2, D + j * dy + dy/2));
    }

  if (c && d)
    {
      if (msh->ptempl[(j - 1) * (n + 1) + i + 1])
        {
          r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
          r += 1./12. * 1./4. * 1./2. * (F (l + i * dx + dx/2, D + j * dy) + F (l + i * dx, D + j * dy - dy/2));
        }
    }

  if (e && d)
    {
      r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
      r += 1./12. * 1./4. * 1./2. * (F (l + i * dx - dx/2, D + j * dy - dy/2) + F (l + i * dx, D + j * dy - dy/2));
    }

  if (e && f)
    {
      r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
      r += 1./12. * 1./4. * 1./2. * (F (l + i * dx - dx/2, D + j * dy - dy/2) + F (l + i * dx - dx/2, D + j * dy));
    }

  if (a && f)
    {
      if (msh->ptempl[(j + 1) * (n + 1) + i - 1])
        {
          r += 1./12. * 1./4. * F (l + i * dx, D + j * dy);
          r += 1./12. * 1./4. * 1./2. * (F (l + i * dx, D + j * dy + dy/2) + F (l + i * dx - dx/2, D + j * dy));
        }
    }
  (void)res;
  return r;
}