Dear Andrew,
Thank you very much for your very detailed answer. I understand that you can shift the derivative of dot(w,grad(p)) to w. But for my system (see attachment), I already have shifted a derivative on w. So I have sth. like dot(grad(w),grad(p)). I do not know what to do in this case…
Henrik
--
Dipl.-Math. Henrik Büsing
Institute for Applied Geophysics and Geothermal Energy
E.ON Energy Research Center
RWTH Aachen University
------------------------------------------------------
Mathieustr. 10 | Tel +49 (0)241 80 49907
52074 Aachen, Germany | Fax +49 (0)241 80 49889
------------------------------------------------------
http://www.eonerc.rwth-aachen.de/GGE
hbuesing@eonerc.rwth-aachen.de
------------------------------------------------------
Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake-bounces@imperial.ac.uk] Im Auftrag von Andrew McRae
Gesendet: 24 September 2015 16:14
Ah, my first instinct was that the FE way would be exactly equivalent, but of course that's not true.
Anyway, lowest-order FE is similar to FV, so here's the FE approach I was thinking of:
The space (mesh, "NCF", 1) has one degree of freedom per facet (http://femtable.org/, 3rd column, bottom panel). We want to find a function in this space which is somehow the gradient of your p in DG0.
It'd be nice if we could just do
u = TrialFunction(W)
w = TestFunction(W)
a = dot(w, u)*dxL = dot(w, grad(p))*dx
out = Function(W)
solve(a == L, out)
However, this clearly won't give the result you want, since p is in DG0, so grad(p) is 0 everywhere, so the right hand side is 0.
Instead, integrate the right hand side by parts to obtain L = p_bc*dot(w, n)*ds - div(w)*p*dx, where n = FacetNormal(mesh), and p_bc is some function with the values of p that you want on the boundary (this can be just p, or you might choose to put p_bc in a continuous space, for example).
Complete code for the 1D version of this is
---
from firedrake import *
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "DG", 0)
f = Function(V).interpolate(Expression("2*x[0] - 1.0"))
f_bc = Function(FunctionSpace(mesh, "CG", 1)).interpolate(Expression("2*x[0] - 1.0"))
W = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(W)
v = TestFunction(W)
a = u*v*dx
n = FacetNormal(mesh)
L = f_bc*v*n[0]*ds - v.dx(0)*f*dx
g = Function(W)
solve(a == L, g)
print g.dat.data
---Indeed, this verifies that the FE 'weak gradient' of 2x-1 is 2 everywhere. For the 3D version, W needs to be FunctionSpace(mesh, "NCF", 1), and L needs to be as given above.
Andrew
On 24 September 2015 at 07:15, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
--
Dipl.-Math. Henrik Büsing
Institute for Applied Geophysics and Geothermal Energy
E.ON Energy Research Center
RWTH Aachen University
------------------------------------------------------
Mathieustr. 10 | Tel +49 (0)241 80 49907
52074 Aachen, Germany | Fax +49 (0)241 80 49889
------------------------------------------------------
http://www.eonerc.rwth-aachen.de/GGE
hbuesing@eonerc.rwth-aachen.de
------------------------------------------------------
Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake-bounces@imperial.ac.uk] Im Auftrag von Andrew McRae
Gesendet: 22 September 2015 20:42
An: firedrake@imperial.ac.uk
Betreff: Re: [firedrake] Cell center on extruded mesh
For the distances between cell centres, yes.
It feels like you're performing some operations that aren't particularly finite element-y? Would you be using the approximate gradient in a variational form immediately, or are there further calculations you want to carry out?
In the former case, I think there's a nice FE way to get what you want.
[Buesing, Henrik] @Andrew: So, what would be the FE way of doing this?
Henrik
In the latter case, we can probably hack something up if we're smart enough...
On 22 September 2015 at 18:42, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Von: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk]" im Auftrag von "Andrew McRae [A.T.T.McRae@bath.ac.uk]
Gesendet: Dienstag, 22. September 2015 18:59
An: firedrake@imperial.ac.uk
Betreff: Re: [firedrake] Cell center on extruded mesh>>Can you say a bit more about how you plan to use this?
>>I wonder to what end you would like these values?
I would like to approximate \grad(pw) on a DG0 Function space by
(p('+') - p('-'))/|x('+')-x('-')|, where x are the cell centers. On a hexahedral mesh this should just give \Delta x, \Delta y, \Delta z, right?
Henrik
On 22 September 2015 at 17:52, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Dear Firedrakers,
on an extruded mesh:
meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True)
mesh = ExtrudedMesh(meshbase, Nz, Delta_z)
horiz_elt = FiniteElement("DG", quadrilateral, 0)
vert_elt = FiniteElement("DG", interval, 0)
elt = OuterProductElement(horiz_elt, vert_elt)
How can I access the element center and the center of the element faces?
Thank you!
Henrik
--
Dipl.-Math. Henrik Büsing
Applied Geophysics and Geothermal Energy
E.ON Energy Research Center
RWTH Aachen University
-----------------------------------------------
Mathieustr. 10 | Tel +49 (0)241 80 49907
52074 Aachen, Germany | Fax +49 (0)241 80 49889
-----------------------------------------------
http://www.eonerc.rwth-aachen.de/GGE
hbuesing@eonerc.rwth-aachen.de
-----------------------------------------------