Re: [firedrake] Interpolating/projecting between meshes
Dear Peter, We don’t have nodal evaluation for arbitrary finite elements. The interpolation operator (both interpolate and Function.at) works for spaces with point evaluation nodes i.e. CG and DG so you can rig up a hack by interpolating the RT field into the spanning DG field (which is lossless), interpolating that into a spanning DG field on the new mesh (you can get the coordinates to interpolate at by interpolating the coordinate field into the corresponding vector DG space), and then finally projecting from the DG space down to RT. This is effectively a variant on: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter... However a nicer approach is probably to use Patrick Farrell’s pefarrell/supermesh-mixed-mass-matrix branch which makes project work between different meshes of the same domain. Finally, I notice you are using C string expressions. We’re about to drop that functionality because it’s superceded by UFL expressions. E.g. x = SpatialCoordinate(mesh) interpolate(u, as_vector(x[0] + 1, x[1]) Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Sentz, Peter" <sentz2@illinois.edu> Date: Thursday, 31 January 2019 at 22:41 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Interpolating/projecting between meshes Hello, In Fenics/Dolfin, I can interpolate between different meshes, even with less basic elements like Raviart-Thomas. For example: from fenics import * mesh1 = UnitSquareMesh(1,1) mesh2 = UnitSquareMesh(2,2) V1 = FunctionSpace(mesh1,'RT',1) V2 = FunctionSpace(mesh2,'RT',1) u = Expression(('x[0] + 1','x[1]'),degree = 1) U = interpolate(u,V1) Uf = interpolate(U,V2) Interpolating the function 'U' that is defined on one mesh can be moved to another in a fairly simple matter. Trying something similar in Firedrake with 'project' instead of 'interpolate' results in an error message. A different error appears for nodal FE functions as well, when using 'interpolate'. However, in that case, I have developed a work-around using mesh.coordinates.dat.data. However, I can't see how to extend this to nodal DG or Raviart-Thomas FE functions. Can you access the coordinates of the degrees of freedom of an arbitrary function space? Barring that, is there some mechanism to move between meshes if I use a MeshHierarchy, for example? Thanks, Peter Sentz
Dear Peter,
On 1 Feb 2019, at 09:35, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Peter,
We don’t have nodal evaluation for arbitrary finite elements. The interpolation operator (both interpolate and Function.at) works for spaces with point evaluation nodes i.e. CG and DG so you can rig up a hack by interpolating the RT field into the spanning DG field (which is lossless), interpolating that into a spanning DG field on the new mesh (you can get the coordinates to interpolate at by interpolating the coordinate field into the corresponding vector DG space), and then finally projecting from the DG space down to RT. This is effectively a variant on: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter...
In addition to doing this, if you have a mesh hierarchy, then the interpolation between meshes in the mesh hierarchy is supported for all elements. For those that are not point-evaluation, it is a little ugly right now: we hope to make it simpler in the future. But you can do this: mesh = ... mh = MeshHierarchy(mesh, ...) Vc = FunctionSpace(mh[0], "RT", 1) Vf = FunctionSpace(mh[1], "RT", 1) interpolator = EmbeddedDGTransfer(Vc.ufl_element()) c = Function(Vc) f = Function(Vf) # Move primal field from coarse to fine interpolator.prolong(c, f) # Move primal field from fine to coarse interpolator.inject(f, c) # Move dual field (a residual) from fine to coarse interpolator.restrict(f, c)
However a nicer approach is probably to use Patrick Farrell’s pefarrell/supermesh-mixed-mass-matrix branch which makes project work between different meshes of the same domain.
Note right now this is restricted to different meshes where you can provide an oracle that relates the cells on the coarse mesh to potentially overlapping cells on the fine mesh (A trivial, although not efficient, such oracle is that every coarse cell potentially overlaps with every fine cell). Thanks, Lawrence
On Fri, Feb 1, 2019 at 4:55 AM Lawrence Mitchell <wencel@gmail.com> wrote:
Dear Peter,
On 1 Feb 2019, at 09:35, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Peter,
We don’t have nodal evaluation for arbitrary finite elements. The interpolation operator (both interpolate and Function.at) works for spaces with point evaluation nodes i.e. CG and DG so you can rig up a hack by interpolating the RT field into the spanning DG field (which is lossless), interpolating that into a spanning DG field on the new mesh (you can get the coordinates to interpolate at by interpolating the coordinate field into the corresponding vector DG space), and then finally projecting from the DG space down to RT. This is effectively a variant on: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter...
In addition to doing this, if you have a mesh hierarchy, then the interpolation between meshes in the mesh hierarchy is supported for all elements. For those that are not point-evaluation, it is a little ugly right now: we hope to make it simpler in the future. But you can do this:
Hi Lawrence, Is there something wrong with the way I do interpolation? Namely, I point locate the quadrature points for the dual basis vectors from the target mesh in the source mesh. I admit that this is not conservative like supermesh, but it should work for every element and be first order accurate (I think). I cannot see why its any harder to do this than the DG route since they both require point location between meshes. I know LibMesh does this hierarchical-in-dimension thing from Leszek, but I can never remember exactly why. Thanks, Matt
mesh = ...
mh = MeshHierarchy(mesh, ...)
Vc = FunctionSpace(mh[0], "RT", 1) Vf = FunctionSpace(mh[1], "RT", 1)
interpolator = EmbeddedDGTransfer(Vc.ufl_element())
c = Function(Vc)
f = Function(Vf)
# Move primal field from coarse to fine interpolator.prolong(c, f)
# Move primal field from fine to coarse interpolator.inject(f, c)
# Move dual field (a residual) from fine to coarse interpolator.restrict(f, c)
However a nicer approach is probably to use Patrick Farrell’s pefarrell/supermesh-mixed-mass-matrix branch which makes project work between different meshes of the same domain.
Note right now this is restricted to different meshes where you can provide an oracle that relates the cells on the coarse mesh to potentially overlapping cells on the fine mesh (A trivial, although not efficient, such oracle is that every coarse cell potentially overlaps with every fine cell).
Thanks,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
On 1 Feb 2019, at 15:49, Matthew Knepley <knepley@gmail.com> wrote:
On Fri, Feb 1, 2019 at 4:55 AM Lawrence Mitchell <wencel@gmail.com> wrote: Dear Peter,
On 1 Feb 2019, at 09:35, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Peter,
We don’t have nodal evaluation for arbitrary finite elements. The interpolation operator (both interpolate and Function.at) works for spaces with point evaluation nodes i.e. CG and DG so you can rig up a hack by interpolating the RT field into the spanning DG field (which is lossless), interpolating that into a spanning DG field on the new mesh (you can get the coordinates to interpolate at by interpolating the coordinate field into the corresponding vector DG space), and then finally projecting from the DG space down to RT. This is effectively a variant on: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter...
In addition to doing this, if you have a mesh hierarchy, then the interpolation between meshes in the mesh hierarchy is supported for all elements. For those that are not point-evaluation, it is a little ugly right now: we hope to make it simpler in the future. But you can do this:
Hi Lawrence,
Is there something wrong with the way I do interpolation? Namely, I point locate the quadrature points for the dual basis vectors from the target mesh in the source mesh. I admit that this is not conservative like supermesh, but it should work for every element and be first order accurate (I think). I cannot see why its any harder to do this than the DG route since they both require point location between meshes.
Basically that is exactly what one needs to do: eat the basis functions with the dual basis. However, FIAT only offers an interface for that on the reference cell. This is fine for affine-mapped elements, but for piola-mapped ones you need to add in appropriate pushfowards/pullbacks of the source mesh basis functions into the target mesh. So it's a little uglier. We don't have a good interface to do this cleanly, so have punted for the last while. Lawrence
Hi, thanks to both of you! I tried David's suggestion first; here's a small code snippet: mesh1 = UnitSquareMesh(1,1) RT = FunctionSpace(mesh1,'RT',1) U = project(as_vector([x**2 -1, y + x*y]),RT_c) DG = VectorFunctionSpace(mesh1,'DG',1) Udg = Function(DG).interpolate(U) mesh2 = UnitSquareMesh(2,2) DG2 = VectorFunctionSpace(mesh2,'DG',1) X = interpolate(mesh2.coordinates,DG2) U_interp = Function(DG2) U_interp.dat.data[:] = Udg.at(X.dat.data_ro) However, it doesn't appear the .at() routine deals correctly with DG functions that should be double-valued at a node. Using the above example, I would expect U_interp to be double-valued at the point (0.5,0.5). (6 degrees of freedom at this node, 3 of them correspond to one value, the other 3 correspond to a different value). However, printing out U_interp.dat.data and comparing it with X.dat.data it looks like it is single-valued there. But maybe this is not exactly what you were suggesting. I will try the other two approaches as well, and thanks again. - Peter Sentz ________________________________________ From: Lawrence Mitchell [wencel@gmail.com] Sent: Friday, February 01, 2019 3:54 AM To: David Ham Cc: Sentz, Peter; firedrake Subject: Re: [firedrake] Interpolating/projecting between meshes Dear Peter,
On 1 Feb 2019, at 09:35, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Peter,
We don’t have nodal evaluation for arbitrary finite elements. The interpolation operator (both interpolate and Function.at) works for spaces with point evaluation nodes i.e. CG and DG so you can rig up a hack by interpolating the RT field into the spanning DG field (which is lossless), interpolating that into a spanning DG field on the new mesh (you can get the coordinates to interpolate at by interpolating the coordinate field into the corresponding vector DG space), and then finally projecting from the DG space down to RT. This is effectively a variant on: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter...
In addition to doing this, if you have a mesh hierarchy, then the interpolation between meshes in the mesh hierarchy is supported for all elements. For those that are not point-evaluation, it is a little ugly right now: we hope to make it simpler in the future. But you can do this: mesh = ... mh = MeshHierarchy(mesh, ...) Vc = FunctionSpace(mh[0], "RT", 1) Vf = FunctionSpace(mh[1], "RT", 1) interpolator = EmbeddedDGTransfer(Vc.ufl_element()) c = Function(Vc) f = Function(Vf) # Move primal field from coarse to fine interpolator.prolong(c, f) # Move primal field from fine to coarse interpolator.inject(f, c) # Move dual field (a residual) from fine to coarse interpolator.restrict(f, c)
However a nicer approach is probably to use Patrick Farrell’s pefarrell/supermesh-mixed-mass-matrix branch which makes project work between different meshes of the same domain.
Note right now this is restricted to different meshes where you can provide an oracle that relates the cells on the coarse mesh to potentially overlapping cells on the fine mesh (A trivial, although not efficient, such oracle is that every coarse cell potentially overlaps with every fine cell). Thanks, Lawrence
participants (5)
-
Ham, David A
-
Lawrence Mitchell
-
Lawrence Mitchell
-
Matthew Knepley
-
Sentz, Peter