Dear Sarraf, What you’re actually doing here is hand-coding the Galerkin projection into DG0. This is a straight up finite element calculation. In both FEniCS and Firedrake the right way to do this is to let the system do it for you: DG0 = FunctionSpace(mesh, 'DG', 0) Error = project(div(solution), DG0) Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> Date: Wednesday, 27 May 2020 at 02:09 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Calculate element-wise mass balance Hello, I am trying to calculate the element-wise mass balance (i.e., \int_{\Omega_element} div(u) d\Omega) over every cell and project the output to DG0 space. Has this been done or is there a way (e.g., with par_loop() ) to implement this in firedrake? This is how I achieved it in fenics: ------------------------------------------------------------------- cell_domains = CellFunction('size_t', mesh) cell_domains.set_all(0) dx = Measure('dx')[cell_domains] area = [cell.volume() for cell in cells(mesh)] for i in range(mesh.num_cells()): Error.vector()[i] = (1./area[i]) * assemble(div(solution)*dx(i)) ------------------------------------------------------------------- Regards, Sarraf