Problems understanding project function

mathematics
1

I though that I had understood the difference between interpolate and project. The first one looks for a function that whose values coincide at the points of the mesh whereas the second one looks for a function that minimizes some norm. I tried a very simple case:

Lets say that I have a function defined on the interval \left[0,1\right]

f(x)=\begin{cases}\frac{1}{2}\, & x\le \frac{1}{2},\\ x & x\ge \frac{1}{2}\end{cases}

I understand that if I create a mesh with one single cell with vertices \left\lbrace 0, 1\right\rbrace I can construct two different piecewise linear functions. The interpolating one, g would have the values:

x g(x)
0 \frac{1}{2}
1 1

On the other hand, if I project f onto that piecewise linear space, I would expect to obtain a function h whose values at the vertices are

x h(x)
0 \frac{3}{8}=0.375
1 \frac{7}{8}=0.875

as I have computed them (by “hand”) by minimizing the L^2 norm of f-h. However, when I try this at dolfin

from dolfin import IntervalMesh, FunctionSpace, FiniteElement, Function, Expression, project

mesh_1 = IntervalMesh(2,0.,1.)
mesh_2 = IntervalMesh(1,0.,1)

H_1 = FunctionSpace( mesh_1, FiniteElement( family="Lagrange", cell= mesh_1.ufl_cell(), degree=1))
H_2 = FunctionSpace( mesh_2, FiniteElement( family="Lagrange", cell= mesh_2.ufl_cell(), degree=1))


f = Function(H_1)
g = Function(H_2)
h = Function(H_2)

f.interpolate(Expression("x[0]<=0.5 ? 0.5 : x[0]",degree=1))
g.interpolate(f)
h = project(f,V=H_2, mesh=mesh_1)

print(' x   : f(x)')
for (x,y) in zip(mesh_1.coordinates()[:,0],f.compute_vertex_values()):
    print( f'{x:.2f} : {y:.2f}')
print('')
print(' x   : g(x)')
for (x,y) in zip(mesh_2.coordinates()[:,0],g.compute_vertex_values()):
    print( f'{x:.2f} : {y:.2f}')
print('')
print(' x   : h(x)')
for (x,y) in zip(mesh_2.coordinates()[:,0],h.compute_vertex_values()):
    print( f'{x:.2f} : {y:.2f}')

I get:

 x   : f(x)
0.00 : 0.50
0.50 : 0.50
1.00 : 1.00

 x   : g(x)
0.00 : 0.50
1.00 : 1.00

 x   : h(x)
0.00 : 0.40
1.00 : 0.70

What am I missing? ( and what is the meaning of the mesh variable in project?)

Thank you very much!

2

Ok, I think I more or less understand it. As I try to project function f onto the space of function h, first it interpolates f in order to compute the inner products, etc. So that’s why I get the same values for the interpolated and the projected one.

If I want to obtain the correct values, I can do:

h = project(Expression("x[0]<=0.5 ? 0.5 : x[0]",degree=3),H_2)
print('')
print(' x   : h_2(x)')
for (x,y) in zip(mesh_2.coordinates()[:,0],h.compute_vertex_values()):
    print( f'{x:.4f} : {y:.4f}')

which returns:

 x   : h_2(x)
0.0000 : 0.3750
1.0000 : 0.8750

However, I think this is the behaviour one is expecting from project right? I mean, to use the whole cell values, not only the interpolated ones.

3

I would suggest using ufl.conditional with SpatialCoordinate to create something that is evaluated at quadrature points, see for instance: Unable to use conditional module - #2 by dokken

4

I understand that that would be the solution to correctly interpolate an expression onto a given functional space.

My question is. Assuming that I have a fem function, defined on some functional space. How can I get project to use the cuadrature points?

I don’t know if I’m explaining my self. How can I get to project the function f from my example onto the functional space H_2, using the cuadrature points so I get the correct projection (and without using the expression of f, as, in general, I would have no expression )

5

Summarizing:

mesh_1 = IntervalMesh(2,0.,1.)
mesh_2 = IntervalMesh(1,0.,1)
H_1 = FunctionSpace( mesh_1, FiniteElement( family="Lagrange", cell= mesh_1.ufl_cell(), degree=1))
H_2 = FunctionSpace( mesh_2, FiniteElement( family="Lagrange", cell= mesh_2.ufl_cell(), degree=1))
f = Function(H_1)
g = Function(H_2)

such that:

f.compute_vertex_values()

returns

array([0.5, 0.5, 1. ])

Is there any way to project f onto H_2 :

g = project(f,...)

such that, after doing that:

g.compute_vertex_values()

returns

array([0.375, 0.875 ])

( assuming I don’t know the expression of f but only its values on the vertices)

6

I would suggest changing the quadrature order of the projection, as you are trying to project a function of higher regularity (f) onto a space of lower regularity (g).

2 Likes
7

Thanks, I think that is what I need. I don’t know however how to do it. I’ve tried:

h = project(f,V=H_2, form_compiler_parameters = {"quadrature_degree", 2})

which gives me an error:

ValueError: dictionary update sequence element #0 has length 17; 2 is required

Which is the correct syntax? is there a place in the documentation where I can look for the compiler_form_parameters?

I’m using dolfin version 2019.2.0.dev0 (in case it helps)

Thank you very much

8

You would need to write out the projection operator by hand, see for instance: Projection of user expression onto DG1 - #4 by kamensky

1 Like
9

Thank you very much, that explains it