Hi Sarraf, apologies for the low responses. I have some comments inline below.
On 14 Sep 2020, at 08:01, Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> wrote:
Hello,
I appreciate if you help me with this question.
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 . .
We can't directly integrate like this. You can do a whole-field operation to get the fluxes into a DG0 trace space: Vt = FunctionSpace(mesh, "HDiv Trace", 0) fluxes = Function(Vt) # you have to code this trace projection by hand: v = TestFunction(Vt) F = (fluxes*v)('+')*dS + fluxes*v*ds - avg(inner(grad(sol), n))*v('+')*dS - inner(grad(sol), n))*v*ds solve(F == 0, fluxes) # Now fluxes contains your integrals. Lawrence