Hello, I am trying to implement a postprocessing flux-limiting scheme for a DG advection problem to eliminate the under/overshoots that appears in the vicinity of the front. In order to achieve this, I need to iterate over each mesh element and then on each facet of the element calculate the value of numerical flux. Based on the sign of this Flux value, a constant edge-specific alpha is calculated and the new corrected solutions are then reconstructed for the cell. In firedrake, Is there any way to iterate over facets of each cell, and calculate an integral of an expression (e.g, avg(inner(grad(sol,n))) ) on that facet? This is how I achieved it in fenics: ----------------------------------------------------- facet_domains = MeshFunction('size_t',mesh,mesh.topology().dim()-1) dS = Measure('dS', subdomain_data = facet_domains) elnum = 0 for cell in cells(mesh): flux=np.zeros(3) for fct in facets(cell): flux[fct.index] = assemble( (avg(inner(grad(sol) , n) ) )*dS(fct.index()) print(flux) . . do some operations based on flux[fct.index] values . . elnum += 1 -------------------------------------------------------- Thank you, Sarraf