"project" is giving nodal interpolate rather than projection

8 months ago by
I have an FEM function that I'd like to project onto a coarser mesh.  I'd like the L2 projection, not the nodal interpolation.  The project function, however, is giving the nodal interpolate.   Sample code: 

# Test L2 projection in fenics
# Author : Paul Barbone
# LastEdit date : 20 Dec 2017

from dolfin import *

# Create mesh
mesh1 = IntervalMesh(100,0.0, 1.0)
mesh2 = IntervalMesh( 2,0.0, 1.0)

# Create function spaces
V1 = FunctionSpace(mesh1, "CG", 1) # Fine scale function space
V2 = FunctionSpace(mesh2, "CG", 1) # Coarse scale function space

my_dof_to_vertex_map = dof_to_vertex_map(V2)

f = Expression("sin(x[0]*pi)", degree=3)

# Project f expression onto fine scale space:
f1 = project(f, V1)

# Project f expression onto coarse scale space:
f2p = project(f, V2)

# Interpolate f expression onto coarse scale space:
f2i = interpolate(f,V2)

# Project fine-scale f onto coarse scale space:
f12p = project(f1,V2)

# Interpolate fine-scale f onto coarse scale space:
f12i = interpolate(f1,V2)

print 'Nodal values of various projections'

print 'projected expression'
print f2p.vector().array()[my_dof_to_vertex_map]
print 'interpolated expression'
print f2i.vector().array()[my_dof_to_vertex_map]
print 'projected fine mesh function'
print f12p.vector().array()[my_dof_to_vertex_map]
print 'interpolated fine mesh function'
print f12i.vector().array()[my_dof_to_vertex_map]

Nodal values of various projections
projected expression
[ 0.11519238 1.15932667 0.11519238]
interpolated expression
[ 0.00000000e+00 1.00000000e+00 1.22464680e-16]
projected fine mesh function
[ 5.96786968e-07 1.00008225e+00 5.96786968e-07]
interpolated fine mesh function
[ 5.96786967e-07 1.00008225e+00 5.96786968e-07]

Some observations about these results:
  1. First of all, I'd hope that the first and third results (the two projections) would approximately match, but they don't.
  2. I'd expect the second and fourth (the nodal interpolations) to agree, and they do.
  3. The third is supposed to be an L2 projection, but it's either (a) an interpolation, or (b) (equivalently in this case) a projection based on a nodal quadrature rule using the nodes in the coarse scale.   (Is nodal quadrature the default in fenics?)
Thanks for any suggestions.
Community: FEniCS Project

1 Answer

8 months ago by
Unfortunately, if you project a function from the fine mesh onto a coarse mesh, it first gets interpolated onto the coarse mesh and then the corresponding variational problem is solved. I had the same problem and as far as I know there is no nice workaround.

Edit: See here https://www.allanswered.com/post/jrwoj/does-project-not-conserve-integral/
Thanks, Lukas.  Your answer would explain the behavior I'm seeing.  If projecting first interpolates, however, then the variational projection step is unnecessary.  The projection of the interpolant is always the original interpolant.  In that case, the project command is doing a lot of unnecessary computing.
written 7 months ago by Paul Barbone  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »