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 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? Thanks a lot, Eike -- Dr Eike Hermann Mueller Research Associate (PostDoc) Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5803 e.mueller@bath.ac.uk http://people.bath.ac.uk/em459/