September 2015 - firedrake

PyOP2 branch with logging support
by Eike Mueller 30 Sep '15

30 Sep '15

Hi Lawrence, I pulled your PyOP2 branch with logging support and updated by firedrake modules+dependencies (ffc, ufl, coffee, fiat - all master), and I get the following when I try to run firedrake-clean: eikemueller@Eikes-MacBook-Pro $ ./scripts/firedrake-clean Traceback (most recent call last): File "./scripts/firedrake-clean", line 3, in <module> from firedrake.ffc_interface import clear_cache, FFCKernel File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/__init__.py", line 31, in <module> from pyop2.logger import set_log_level, info_red, info_green, info_blue, log # noqa File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/__init__.py", line 13, in <module> from op2 import * File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/op2.py", line 39, in <module> import base File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/base.py", line 60, in <module> from coffee.visitors import FindInstances, EstimateFlops ImportError: No module named visitors To me it looks like PyOP2 tries to import a non-existant coffee module, but that’s weird, since there is a subdirectory coffee/visitors in my COFFEE installation, and the __init__.py contains: eikemueller@Eikes-MacBook-Pro $ cat __init__.py from __future__ import absolute_import from coffee.visitors.utilities import * from coffee.visitors.inspectors import * and FindInstances is defined in inspectors.py eikemueller@Eikes-MacBook-Pro $ grep FindInstances *py inspectors.py: "SymbolDeclarations", "FindInstances"] inspectors.py:class FindInstances(Visitor): inspectors.py: super(FindInstances, self).__init__() Thanks, Eike -- Dr Eike Hermann Mueller Lecturer in Scientific Computing Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 6470 e.mueller(a)bath.ac.uk<mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/

2 4

Validate exact Jacobian by Finite Differences
by Buesing, Henrik 30 Sep '15

30 Sep '15

Dear Firedrakers, I compared my Jacobian generated by finite differences with the Jacobian generated by firedrake. The firedrake Jacobian seems to be wrong... I generated the FD Jacobian in my Fortran code. Is there a way to simply calculate the FD Jacobian within firedrake to validate the exact Jacobian? Thank you! 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(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 1

Inspecting values of Jacobian matrix
by Buesing, Henrik 30 Sep '15

30 Sep '15

Dear Firedrakers, I would like to inspect the values of my Jacobian matrix derived from my form F. I did a J=derivative(F,u) print assemble(J).M But this did not result in the actual values (see below for output). Can I somehow obtain the actual values? Thank you! Henrik Output of print assemble(J).M OP2 Mat: None_None_None_None_matrix, sparsity (OP2 Sparsity: dsets (MixedDataSet((DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), (1,), 'None_nodes_dset'), DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), (1,), 'None_nodes_dset'))), MixedDataSet((DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), (1,), 'None_nodes_dset'), DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), (1,), 'None_nodes_dset')))), rmaps (MixedMap((DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None), DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None))), MixedMap((DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), None), DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), None)))), cmaps (MixedMap((DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None), DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None))), MixedMap((DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), None), DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), None)))), name None_None_None_None_sparsity), datatype float64 -- 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(a)eonerc.rwth-aachen.de ------------------------------------------------------

3 4

Re: [firedrake] Inspecting values of Jacobian matrix
by Fabio Luporini 30 Sep '15

30 Sep '15

Hi Henrik You should be able to print the values by doing *print assemble(J).M.values*If thematrix is too big a summary will be printed. If you want to print more values (don't know if actually useless, but just in case), you can do import numpy as np np.set_printoptions(threshold=...) and setting the threshold to higher values should do the trick Thanks -- Fabio 2015-09-30 9:06 GMT+01:00 Buesing, Henrik <HBuesing(a)eonerc.rwth-aachen.de>: > Dear Firedrakers, > > > > I would like to inspect the values of my Jacobian matrix derived from my > form F. I did a > > > > > *J=derivative(F,u) print assemble(J).M* > > > > But this did not result in the actual values (see below for output). Can I > somehow obtain the actual values? > Thank you! > > Henrik > > > > Output of *print assemble(J).M* > > > > OP2 Mat: None_None_None_None_matrix, sparsity (OP2 Sparsity: dsets > (MixedDataSet((DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), > (1,), 'None_nodes_dset'), DataSet(ExtrudedSet(Set((4, 4, 4, 4), > 'None_nodes'), 5), (1,), 'None_nodes_dset'))), > MixedDataSet((DataSet(ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), > (1,), 'None_nodes_dset'), DataSet(ExtrudedSet(Set((4, 4, 4, 4), > 'None_nodes'), 5), (1,), 'None_nodes_dset')))), rmaps > (MixedMap((DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), > 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), > 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), > frozenset(['ALL']), frozenset([]), None), > DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, > 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None))), > MixedMap((DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, > 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), > None), DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, > 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), > None)))), cmaps (MixedMap((DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), > 'plexmesh_interior_facets'), 5), ExtrudedSet(Set((4, 4, 4, 4), > 'None_nodes'), 5), 2, None, 'None_interior_facet_node'), > frozenset(['ALL']), frozenset([]), None), > DecoratedMap(Map(ExtrudedSet(Set((0, 0, 0, 0), 'plexmesh_interior_facets'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 2, None, > 'None_interior_facet_node'), frozenset(['ALL']), frozenset([]), None))), > MixedMap((DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, > 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), > None), DecoratedMap(Map(ExtrudedSet(Set((1, 1, 1, 1), 'plexmesh_cells'), > 5), ExtrudedSet(Set((4, 4, 4, 4), 'None_nodes'), 5), 1, None, > 'None_cell_node'), frozenset(['ALL', 'ON_INTERIOR_FACETS']), frozenset([]), > None)))), name None_None_None_None_sparsity), datatype float64 > > > > > > -- > > 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(a)eonerc.rwth-aachen.de > > ------------------------------------------------------ > > >

1 0

Inspecting values of potential function
by Buesing, Henrik 29 Sep '15

29 Sep '15

Dear Firedrakers, I define a potential function on a DG function space (see attached file). Now I would like to inspect the values on my mesh. But when I do a print assemble(psiw).vector().array() I get an error (see below). In contrast evaluating the integral over the facets print assemble(psiw*dS_h) works just fine. How can I inspect the values of the potential function on every grid cell? Thank you! Henrik Error: ValueError: Expression has incompatible function spaces <firedrake.functionspace.FunctionSpace object at 0x7f3c67723390> and <firedrake.functionspace.MixedFunctionSpace object at 0x7f3c67723a50> -- 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(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 2

More questions
by Justin Chang 29 Sep '15

29 Sep '15

Hi everyone, Some general questions about Firedrake: 1) I ran the poisson_mixed.py example and noticed in the Memory Usage section of -log_summary that there are more creations than destructions. In my experience, not properly destroying PETSc objects can inhibit performance. Is this an issue for you guys? As in, are you guys properly taking care of the memory allocation and related issues? 2) How do I read in an auxiliary data from a datafile? In FEniCS I could load data as a numpy array and use CellFunction() to "align" the data with the connectivity/ordering of the cells. 3) Is there any update regarding adding parallel IO? Last I spoke with Matt he said that he and one or two of you Firedrake folks (Michael Lange I think?) were working on this. 4) Do you guys support Hexahedron elements? When I look through the documentation for UnitCubeMesh and BoxMesh, I see no option for non-tetrahedron elements. Thanks, Justin

2 5

Re: [firedrake] Low order DG for pure diffusion?
by Colin Cotter 29 Sep '15

29 Sep '15

Dear Justin, Since we don't support elimination within hybridized methods for DG at present, your best option is to use the interior penalty method since it can be applied in primal form. See for example: Arnold, D. N., Brezzi, F., Cockburn, B., & Marini, L. D. (2002). Unified analysis of discontinuous Galerkin methods for elliptic problems. *SIAM journal on numerical analysis*, *39*(5), 1749-1779. all the best -cjc On 29 September 2015 at 08:48, Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > I am doing an operator-splitting time-stepping scheme for an > advection-diffusion-reaction system. That said, I have two maybe related > questions. > > 1) I see that you guys have a nice low order DG scheme for the advection > equation, and was wondering if a similar discretization can be made for the > transient diffusion equation. > > I have seen plenty of nodally discontinuous galerkin methods like SIPG, > but is it possible to discretize that equation using DG0 elements (i.e., > piece-wise constant)? In other words, I am wondering if it’s possible to > come up with a finite-volume-like weak form of the following: > > du/dt - div(D*grad(u)) = f > > where D is the diffusivity tensor (dim by dim) and we use backward euler > time-stepping. > > 2) Suppose I use the DG0 formulation for advection and CG1 for diffusion, > in operator splitting, the “initial condition” for one operator is the > computed solution from the previous operator. That is, if I first solve the > advection equation for one time step, that solution becomes the initial > condition for the diffusion equation. If I have to pass back and forth > between DG0 and CG1 solution spaces, will I lose a lot of accuracy due to > projecting or “smoothening” out the data? > > Thanks, > Justin > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 2

Low order DG for pure diffusion?
by Justin Chang 29 Sep '15

29 Sep '15

Hi all, I am doing an operator-splitting time-stepping scheme for an advection-diffusion-reaction system. That said, I have two maybe related questions. 1) I see that you guys have a nice low order DG scheme for the advection equation, and was wondering if a similar discretization can be made for the transient diffusion equation. I have seen plenty of nodally discontinuous galerkin methods like SIPG, but is it possible to discretize that equation using DG0 elements (i.e., piece-wise constant)? In other words, I am wondering if it’s possible to come up with a finite-volume-like weak form of the following: du/dt - div(D*grad(u)) = f where D is the diffusivity tensor (dim by dim) and we use backward euler time-stepping. 2) Suppose I use the DG0 formulation for advection and CG1 for diffusion, in operator splitting, the “initial condition” for one operator is the computed solution from the previous operator. That is, if I first solve the advection equation for one time step, that solution becomes the initial condition for the diffusion equation. If I have to pass back and forth between DG0 and CG1 solution spaces, will I lose a lot of accuracy due to projecting or “smoothening” out the data? Thanks, Justin

1 0

Re: [firedrake] Firedrake issue
by David Ham 28 Sep '15

28 Sep '15

Hi Tomasz, I've gone over your working, and I can't understand how your weak form works. The test function \delta phi seems to simultaneously be defined over the xy free surface and the yz coupling surface. Since each discrete equation corresponds to one matrix row, I don't understand to which matrix rows the contributions are supposed to contribute. Can you enlighten me? Regards, David On Thu, 24 Sep 2015 at 08:06 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > Hi David, > > > Concerning the FD update, it will take a while, as I encounter problems > installing petsc4py and I need to wait for IT guy. So I'm working on the > version of 15/06/15 for the moment. > > > I'm new to github, but this should work: > https://github.com/mmtjs/lin_coupled > <https://github.com/mmtjs/lin_coupled> > mmtjs/lin_coupled > Contribute to lin_coupled development by creating an account on GitHub. > Read more... <https://github.com/mmtjs/lin_coupled> > I had an error in the previous formulation, so I redirived it and at the > moment the solver is partitioned i.e. it doesn't require one matrix for the > whole system, but just mappings of corresponding boundary nodes from one > domain to the other (these work). > > > The issue is in lines 175-179, where the coupling terms are commented out > ... > a_eta = trial * v * ds_t > L_eta = ( eta + dt*dot( n_w, grad(phi) ) ) * v * ds_t #+ dt/rho * phi * v > * ds_v(2) - dt * Pw * v * ds_v(2) > ... > More on it in the previous pdf file. > > Cheers, > Tomasz > > > ------------------------------ > *From:* David Ham <David.Ham(a)imperial.ac.uk> > *Sent:* 23 September 2015 10:27 > *To:* Tomasz Salwa [RPG]; firedrake > *Subject:* Re: Firedrake issue > > Hi Tomasz, > > Welcome back. > > Let's get through a few preliminaries. > > First, if you haven't used firedrake for 2-3 months then it will be quite > out of date. You need to update to the current version. You can do this > using your existing install and update all the packages, or you can use the > new installer to get a new installation which will also have an attached > update script. See http://firedrakeproject.org/download.html > Obtaining Firedrake — Firedrake 0.11.0-658-ge79e125-dirty ... > Running firedrake-install with no arguments will install firedrake in a > python virtualenv created in a firedrake subdirectory of the current > directory. > Read more... <http://firedrakeproject.org/download.html> > > Next, it's much better to have this sort of conversation on the Firedrake > list. This keeps everyone in the loop and enables other people to intervene > when I say something stupid ;). I have taken the liberty of adding you to > the list and I'm replying to the list now so that this conversation ends up > there. > > Now, we need to be talking about code, not only maths, so please push what > you're working on to a git repository on github or bitbucket so we can see > what's going on. > > Finally we can start to think about the actual problem. I think that last > time we talked about this, we discussed setting up a mixed mass matrix > across the two domains and checking that projection works. Does that still > work (and can you post some code to a repository which does this)? > > Once we get past that then we need to talk about how to assemble forms on > the boundary facets. This will involve doing some mucking around with maps, > because at the moment the facets belong to two different meshes. We need to > think a little about how to do that, but we can do so when the first steps > work. > > Cheers, > > David > > > On Tue, 22 Sep 2015 at 12:54 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > >> Dear David, >> >> >> I was busy in the recent months and could not work on the fluid-structure >> interaction problem. Now I returned to it and reformulated it. >> >> >> I am getting the following error: >> >> >> RuntimeError: Nonlinear solve failed to converge after 0 nonlinear >> iterations. >> Reason: >> Inner linear solve failed to converge after 10000 iterations with >> reason: DIVERGED_MAX_IT >> >> >> but the problem is with the proper formulation in Firedrake. >> >> I include the pdf with it. Would you be able to help me with this? >> >> >> Regards, >> >> Tomasz >> >> >> >> ----- >> >> Tomasz Salwa >> >> PhD Research Student >> >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~mmtjs >> >

2 1

Re: [firedrake] Cell center on extruded mesh
by Andrew McRae 28 Sep '15

28 Sep '15

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! Have you tried looking at the literature for FE discretisations of your problem? Best, Andrew On 28 September 2015 at 08:17, Buesing, Henrik < HBuesing(a)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 > > 52074 Aachen, Germany | Fax +49 (0)241 80 49889 > > ------------------------------------------------------ > > http://www.eonerc.rwth-aachen.de/GGE > > hbuesing(a)eonerc.rwth-aachen.de > > ------------------------------------------------------ > > > > *Von:* firedrake-bounces(a)imperial.ac.uk [mailto: > firedrake-bounces(a)imperial.ac.uk] *Im Auftrag von *Andrew McRae > *Gesendet:* 24 September 2015 16:14 > > *An:* firedrake(a)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(a)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(a)eonerc.rwth-aachen.de > > ------------------------------------------------------ > > > > *Von:* firedrake-bounces(a)imperial.ac.uk [mailto: > firedrake-bounces(a)imperial.ac.uk] *Im Auftrag von *Andrew McRae > *Gesendet:* 22 September 2015 20:42 > *An:* firedrake(a)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(a)eonerc.rwth-aachen.de> wrote: > > > ------------------------------ > > *Von:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk]" > im Auftrag von "Andrew McRae [A.T.T.McRae(a)bath.ac.uk] > *Gesendet:* Dienstag, 22. September 2015 18:59 > *An:* firedrake(a)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(a)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(a)eonerc.rwth-aachen.de > ----------------------------------------------- > > > > > > >

2 1