Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake-bounces@imperial.ac.uk] Im Auftrag von Andrew McRae Gesendet: 28 September 2015 09:31 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Cell center on extruded mesh Hi Henrik, I'm afraid I have no suggestions, and I'd rather not get into the sport of finding discretisations for PDEs that I know nothing about! [Buesing, Henrik] True! :-) Have you tried looking at the literature for FE discretisations of your problem? [Buesing, Henrik] Well, there is a book<http://www.amazon.de/gp/product/089871656X?psc=1&redirect=true&ref_=oh_aui_detailpage_o00_s00> on DG methods, which has a chapter on porous media flow. I ordered that one. Then people often perform a decoupling of the two equations into a pressure equation (solved implicitly) and a saturation equation. There are a lot of discretisations for this one. But I’d rather keep the full system, since I’m also interested into slightly different equations, which govern miscibility of the two phases. Best, Andrew On 28 September 2015 at 08:17, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de<mailto:HBuesing@eonerc.rwth-aachen.de>> wrote: 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<tel:%2B49%20%280%29241%2080%2049907> 52074 Aachen, Germany | Fax +49 (0)241 80 49889<tel:%2B49%20%280%29241%2080%2049889> ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de<mailto:hbuesing@eonerc.rwth-aachen.de> ------------------------------------------------------ Von: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [mailto:firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] Im Auftrag von Andrew McRae Gesendet: 24 September 2015 16:14 An: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Betreff: Re: [firedrake] Cell center on extruded mesh 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)*dx L = 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<mailto: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<tel:%2B49%20%280%29241%2080%2049907> 52074 Aachen, Germany | Fax +49 (0)241 80 49889<tel:%2B49%20%280%29241%2080%2049889> ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de<mailto:hbuesing@eonerc.rwth-aachen.de> ------------------------------------------------------ Von: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [mailto: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<mailto: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<mailto:HBuesing@eonerc.rwth-aachen.de>> wrote: ________________________________ Von: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>]" im Auftrag von "Andrew McRae [A.T.T.McRae@bath.ac.uk<mailto:A.T.T.McRae@bath.ac.uk>] Gesendet: Dienstag, 22. September 2015 18:59 An: firedrake@imperial.ac.uk<mailto: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<mailto: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<tel:%2B49%20%280%29241%2080%2049907> 52074 Aachen, Germany | Fax +49 (0)241 80 49889<tel:%2B49%20%280%29241%2080%2049889> ----------------------------------------------- http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de<mailto:hbuesing@eonerc.rwth-aachen.de> -----------------------------------------------