Re: [firedrake] Projecting CG1 function into DG0 function
Basically: d = TrialFunction(D.function_space()) e = TestFunction(D.function_space()) solve(d * e * dx == div(q) * e * dx, D) One might have to think a little about what value this gives at the boundary, but I think it's OK. On Thu, 12 Nov 2015 at 11:00 Justin Chang <jychang48@gmail.com> wrote:
Hi all,
Perhaps this may be a simple question, but say I have this bilinear and linear form:
a = grad(u)*D*grad(v)*dx L = F*v*dx
where u,v is trial/test function on CG1 space, and D and F are coefficients that live in DG0 space.
Say I have a velocity vector function q (CG1) on the same mesh and want D to be the element-wise divergence of said velocity.
How would I formulate D in Firedrake?
Thanks, Justin
If you're worried about performance, then you can notice that the LHS matrix is diagonal (because D is in DG0) so you can avoid matrix assembly and solve with: D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) (ie just divide by the diagonal entries of the matrix). On Thu, 12 Nov 2015 at 11:08 David Ham <David.Ham@imperial.ac.uk> wrote:
Basically:
d = TrialFunction(D.function_space()) e = TestFunction(D.function_space())
solve(d * e * dx == div(q) * e * dx, D)
One might have to think a little about what value this gives at the boundary, but I think it's OK.
On Thu, 12 Nov 2015 at 11:00 Justin Chang <jychang48@gmail.com> wrote:
Hi all,
Perhaps this may be a simple question, but say I have this bilinear and linear form:
a = grad(u)*D*grad(v)*dx L = F*v*dx
where u,v is trial/test function on CG1 space, and D and F are coefficients that live in DG0 space.
Say I have a velocity vector function q (CG1) on the same mesh and want D to be the element-wise divergence of said velocity.
How would I formulate D in Firedrake?
Thanks, Justin
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? Thanks, Justin On Thursday, November 12, 2015, David Ham <David.Ham@imperial.ac.uk> wrote:
If you're worried about performance, then you can notice that the LHS matrix is diagonal (because D is in DG0) so you can avoid matrix assembly and solve with:
D.assign(assemble(div(q) * e * dx)/assemble(e * dx))
(ie just divide by the diagonal entries of the matrix).
On Thu, 12 Nov 2015 at 11:08 David Ham <David.Ham@imperial.ac.uk <javascript:_e(%7B%7D,'cvml','David.Ham@imperial.ac.uk');>> wrote:
Basically:
d = TrialFunction(D.function_space()) e = TestFunction(D.function_space())
solve(d * e * dx == div(q) * e * dx, D)
One might have to think a little about what value this gives at the boundary, but I think it's OK.
On Thu, 12 Nov 2015 at 11:00 Justin Chang <jychang48@gmail.com <javascript:_e(%7B%7D,'cvml','jychang48@gmail.com');>> wrote:
Hi all,
Perhaps this may be a simple question, but say I have this bilinear and linear form:
a = grad(u)*D*grad(v)*dx L = F*v*dx
where u,v is trial/test function on CG1 space, and D and F are coefficients that live in DG0 space.
Say I have a velocity vector function q (CG1) on the same mesh and want D to be the element-wise divergence of said velocity.
How would I formulate D in Firedrake?
Thanks, Justin
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
Hi guys, So I was attempting to do something similar: L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx)) where v is test function, u is FE solution, and u_exact is the analytical solution. All of which are CG1 space. I get an error saying "fl.log.UFLException: Division by non-scalar is undefined" Know what's up? Thanks, Justin On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
David Ham
-
Justin Chang
-
Lawrence Mitchell