On 12/11/15 11:27, Justin Chang wrote:
David,
So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns cell-wise div(q), what's the denominator "/assemble(e * dx)" for?
It's just the normal FE L2 projection: you want: u = div(q) So you hit both sides with a test function and integrate: u*e*dx = div(q)*e*dx Where u is a trial function in DG0, e is a test function in DG and q is in whereever. But, as David points out, the DG0 mass-matrix is completely diagonal. The values on the diagonal are just obtained by assembling the lone test function: diag = assemble(e*dx) But now, because the matrix is diagonal, you can solve the linear system pointwise by doing a pointwise division: rhs = assemble(div(q)*e*dx) solution = assemble(rhs/diag) Cheers, Lawrence