Interpolate DG functions between meshes
For functions in CG spaces, I usually do the following to interpolate between meshes: meshc = UnitSquareMesh(2,2) mesh = UnitSquareMesh(4,4) Vc = FunctionSpace(meshc,'CG',1) V = FunctionSpace(mesh,'CG',1) xc,yc = SpatialCoordinate(meshc) fc = 1 + xc + yc Fc = Function(Vc).interpolate(fc) W = VectorFunctionSpace(mesh,'CG',1) X = interpolate(mesh.coordinates,W).dat.data_ro Fc_I = Function(V) Fc_I.dat.data[:] = Fc.at(X) I interpolate the mesh coordinates to the second mesh and use the .at() routine to evaluate the function defined on the coarse mesh on those points. It agrees with directly interpolating the expression onto the second mesh (within machine precision). However, if I am using 'DG' function spaces, the .at() routine is single valued, so if I apply this method to a function Fc that has multiple values at a node, it will force all values at the same mesh node to be identical. How would you suggest interpolating DG functions between meshes? Thanks, Peter
On 04/06/2019 21:57, Sentz, Peter wrote:
How would you suggest interpolating DG functions between meshes?
Dear Peter, You can't interpolate DG functions as they're not regular enough for pointwise evaluation to be defined. Instead, you should compute a projection; to do that, one uses a supermesh (a mesh of the intersections of the elements of the two input meshes). libsupermesh is wired up in firedrake, so the following code runs and computes the projection in the L^2 norm: from firedrake import * meshc = UnitSquareMesh(2,2) mesh = UnitSquareMesh(4,4) Vc = FunctionSpace(meshc,'DG',1) V = FunctionSpace(mesh,'DG',1) xc,yc = SpatialCoordinate(meshc) fc = 1 + xc + yc Fc = Function(Vc).interpolate(fc) f = project(Fc, V) Hope this helps, Patrick
In addition to Patrick's response (which is completely correct), I should point out that your example won't show any discontinuities. This is for two reasons, the first is that fc is affine, so its interpolation into Vc is exact and continuous. Then, V spans Vc so the Galerkin projection from Vc into V is the identity, so f = Fc and there will be no new discontinuities here either. If you want to observe some discontinuities to convince yourself that this is actually working, project a nonlinear function into V, then into Vc: x, y = SpatialCoordinate(mesh) f = project(cos(2*pi*x) * sin(2*pi*y), V) fc = project(f, Vc) Regards, David On 05/06/2019, 07:52, "firedrake-bounces@imperial.ac.uk on behalf of Patrick Farrell" <firedrake-bounces@imperial.ac.uk on behalf of patrick.farrell@maths.ox.ac.uk> wrote: On 04/06/2019 21:57, Sentz, Peter wrote: > How would you suggest interpolating DG functions between meshes? > Dear Peter, You can't interpolate DG functions as they're not regular enough for pointwise evaluation to be defined. Instead, you should compute a projection; to do that, one uses a supermesh (a mesh of the intersections of the elements of the two input meshes). libsupermesh is wired up in firedrake, so the following code runs and computes the projection in the L^2 norm: from firedrake import * meshc = UnitSquareMesh(2,2) mesh = UnitSquareMesh(4,4) Vc = FunctionSpace(meshc,'DG',1) V = FunctionSpace(mesh,'DG',1) xc,yc = SpatialCoordinate(meshc) fc = 1 + xc + yc Fc = Function(Vc).interpolate(fc) f = project(Fc, V) Hope this helps, Patrick _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Ham, David A
-
Patrick Farrell
-
Sentz, Peter