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