On 06/11/14 15:24, Eike Mueller wrote:
Dear firedrakers,
when I store the local DG mass matrix in a PyOP2 dat of 3x3 matrices by hand, why do I need to specify the op2.i[0] iteration space as in here
This is part of the PyOP2 kernel API, you can read about it in gory detail here http://op2.github.io/PyOP2/kernels.html#kernel-api. Essentially the kernel that FFC generates writes to a "local iteration space" if you don't indicate this in the map accessor then we pass a pointer instead.
V = FunctionSpace(mesh,'DG',1) V_DG0 = FunctionSpace(mesh,'DG',0)
u = TestFunction(V) v = TrialFunction(V)
mass = u*v*dx
mass_kernel = compile_form(mass, 'mass')[0][6] mass_matrix = Function(V_DG0, val=op2.Dat(V_DG0.node_set**(3*3)))
op2.par_loop(mass_kernel, mass_matrix.cell_set, mass_matrix.dat(op2.INC, mass_matrix.cell_node_map()[op2.i[0]]),
mesh.coordinates.dat(op2.READ,mesh.coordinates.cell_node_map(),flatten=True),
mesh.coordinates.dat(op2.READ,mesh.coordinates.cell_node_map(),flatten=True))
When I print out mass_kernel.code I get
static inline void mass_cell_integral_0_otherwise (double A[3][3] , double **vertex_coordinates , double **w0 ) { […] }
i.e. on each cell the kernel gets passed a 3x3 matrix (= array with 9 elements). If I leave out the "[op2.i[0]]” I only get zeros in mass_matrix. Also, why do the coordinates have to be passed with flatten=True keyword?
Because FFC does not distinguish between mixed and vector finite elements, it always generates code that accesses: x0 x1 x2 y0 y1 y2 z0 z1 z2, etc... I.e. it expects all the x components of the vector to appear before all the y components and so on. However, we store the data as: x0 y0 z0 x1 y1 z1 x2 y2 z2 If you pass flatten=False (the default) the generated code passes a pointer to the first "vector" entry in the Dat. i.e. you can access x0 y0 z0 x1 y1 z1 x2 y2 z2 ^ ^ ^ | | | w0[0] w0[1] w0[2] and then y0 is accessed as w0[0][1] FFC wants y0 to be at w0[1][0] so flatten=True tells PyOP2 to generate code for pointers to every entry. Diagram also on the referenced webpage above. Cheers, Lawrence