Re: [firedrake] Cell center on extruded mesh
Can you say a bit more about how you plan to use this? One way to get the coordinates of the cell centre is to make a DG0 VectorFunctionSpace, and 'interpolate' the coordinate field into it: V = VectorFunctionSpace(mesh, "DG", 0) f = Function(V) f.interpolate(Expression(("x[0]", "x[1]", "x[2]"))) This gives a field with one vector-valued degree of freedom per cell, where the value is [or, happens to be] the centre of the cell. A more flexible solution is to write your own Firedrake parloop using the cell coordinates. Again, it would be helpful if you could say how you plan to use this. 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 -----------------------------------------------
________________________________ 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<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> -----------------------------------------------
participants (2)
-
Andrew McRae
-
Buesing, Henrik