gistfile1.py

from firedrake import *
mesh = IcosahedralSphereMesh(radius=1, refinement_level=2)

kernel = op2.Kernel("""void outward_normal (double** coords , double** normal )
{
  double v0[3] ;
  double v1[3] ;
  double n[3] ;
  double x[3] ;
  double dot ;
  double norm ;
  norm = 0.0;
  dot = 0.0;
  int i ;
  int c ;
  int d ;
  for (i = 0; i < 3; ++i)
  {
    v0[i] = coords[1][i] - coords[0][i];
    v1[i] = coords[2][i] - coords[0][i];
    x[i] = 0.0;
  }
  /* Compute normal to triangle (v0 x v1)*/
  n[0] = v0[1] * v1[2] - v0[2] * v1[1];
  n[1] = v0[2] * v1[0] - v0[0] * v1[2];
  n[2] = v0[0] * v1[1] - v0[1] * v1[0];
  /* Compute a reference "outward" direction */
  for (i = 0; i < 3; ++i)
  {
    x[0] += coords[i][0];
    x[1] += coords[i][1];
    x[2] += coords[i][2];
  }
  /* Cell may be back to front, so figure out dot(n, outward) */
  dot += (x[0]) * n[0];
  dot += (x[1]) * n[1];
  dot += (x[2]) * n[2];
  norm += n[0] * n[0];
  norm += n[1] * n[1];
  norm += n[2] * n[2];
  norm = sqrt(norm);
  norm *= (dot < 0 ? -1 : 1);
  /* Write normal to output function */
  for (c = 0; c < 3; ++c)
  {
      normal[0][c] = n[c]/norm;
  }
}""", "outward_normal")


coords = mesh.coordinates

normal_fs = VectorFunctionSpace(mesh, 'DG', 0)
normal = Function(normal_fs)
op2.par_loop(kernel, normal.cell_set,
             coords.dat(op2.READ, coords.cell_node_map()),
             normal.dat(op2.WRITE, normal.cell_node_map()))

out = File('out.vtu')

out << normal