How to project a function using values at quadrature points only?

407
views
0
8 months ago by
I have a python function    $f\left(\vec{x}\right)$ƒ (x)  , which I can evaluate at discrete points but for which I cannot give a closed form. Now I want to approximate this function on a function space using Lagrangian basis functions. It is also important that the resulting function can be evaluated at arbitrary points within the domain.

 Removed incorrect part. For clarification of the question, please see below.
Community: FEniCS Project
Can you clarify the meaning of "project[ing] a function using values at quadrature points only"?
written 8 months ago by pf4d

Sorry, I was not very clear on what I want to achieve.

Basically, I try to compute a function   $g$g  within function space V for which the difference   $f-g$ƒ g  is orthogonal to function space V, i.e.   $\int_{\left[0,1\right]^{^2}}f\left(\vec{x}\right)\Phi\left(\vec{x}\right)d\vec{x}=\int_{\left[0,1\right]^{^2}}g\left(\vec{x}\right)\Phi\left(\vec{x}\right)d\vec{x}$[0,1]2ƒ (x)Φ(x)dx=[0,1]2g(x)Φ(x)dx for all test functions  $\Phi\left(\vec{x}\right)\in V$Φ(x)V . While I do not have a closed form of  $f$ƒ , I can provide values for  $f$ƒ at given points, for example quadrature points.

written 8 months ago by Philipp O
You can define an Expression as shown by @pf4d where you call your own python function inside eval. Then you can use project.
class F(Expression):
def eval(self, value, x):
value[0] = fun(x)

f = F(degree=2)

mesh = UnitSquareMesh(10, 10)
V    = FunctionSpace(mesh, 'CG', 1)
u    = project(f, V)​

Now u-f is orthogonal to V.

written 8 months ago by Praveen C
Thank you for this answer. I tried out this version but it seems that in this case project does the same as interpolate, as can be seen from this version:
from dolfin import *
import numpy as np

fun = lambda x: np.sin(2*np.pi*x[0]) \
* np.sin(2*np.pi*x[1])

class F(Expression):
def eval(self, value, x):
value[0] = fun(x)
print( x )

print( "before f = F" )
f = F(degree=3)

mesh = UnitSquareMesh(1, 1)
V    = FunctionSpace(mesh, 'CG', 3)

print( "before project" )
u    = project(f, V)

print( "before dof coordinates" )
for p in V.tabulate_dof_coordinates().reshape( (-1,2) ):
print( p )​

This results in output:

before f = F
before project
[ 0.  0.]
[ 1.  0.]
[ 1.  1.]
[ 1.          0.33333333]
[ 1.          0.66666667]
[ 0.33333333  0.33333333]
[ 0.66666667  0.66666667]
[ 0.33333333  0.        ]
[ 0.66666667  0.        ]
[ 0.66666667  0.33333333]
[ 0.  0.]
[ 0.  1.]
[ 1.  1.]
[ 0.33333333  1.        ]
[ 0.66666667  1.        ]
[ 0.33333333  0.33333333]
[ 0.66666667  0.66666667]
[ 0.          0.33333333]
[ 0.          0.66666667]
[ 0.33333333  0.66666667]
before dof coordinates
[ 0.  1.]
[ 0.  0.]
[ 1.  1.]
[ 0.33333333  0.33333333]
[ 0.66666667  0.66666667]
[ 0.33333333  1.        ]
[ 0.66666667  1.        ]
[ 0.          0.33333333]
[ 0.          0.66666667]
[ 0.33333333  0.66666667]
[ 1.  0.]
[ 1.          0.33333333]
[ 1.          0.66666667]
[ 0.33333333  0.        ]
[ 0.66666667  0.        ]
[ 0.66666667  0.33333333]

The function f is only evaluated at the coordinates associated with the degrees of freedom. This does not change if I set a higher quadrature degree using:

parameters["form_compiler"]["quadrature_degree"] = 10


Is there a different parameter controlling the way, project works?

What seems to work is the following, actually solving the variational problem:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

fun = lambda x: np.sin(2*np.pi*x[0]) \
* np.sin(2*np.pi*x[1])

mesh = UnitSquareMesh( 2, 2 )
QE   = FiniteElement( family = "Quadrature", cell = mesh.ufl_cell(), degree = 1, quad_scheme = "default" )
Q    = FunctionSpace( mesh, QE )
g    = TrialFunction( Q )
phi  = TestFunction( Q )
f    = Function( Q )
coor = Q.tabulate_dof_coordinates().reshape( (-1,2) )
f.vector().set_local( np.array( [ fun( p ) for p in coor ] ) )
a    = g * phi * dx
L    = f * phi * dx
u    = Function( Q )
solve( a == L, u )​

While u can be plotted using

vtkFile = File( 'test.pvd' )
vtkFile << u 

function u cannot be evaluated using:

u( [ 0.5, 0.5 ] )

as this results in error message:

Traceback (most recent call last):
File "testProjection.py", line 23, in <module>
u( [ 0.5, 0.5 ] )
File "/Users/philipp/anaconda/envs/fenicsproject/lib/python2.7/site-packages/dolfin/functions/function.py", line 598, in __call__
self.eval(values, x)
RuntimeError: evaluate_basis: Function not supported/implemented for QuadratureElement.

Is there a different way of evaluating u at an arbitrary points in the domain?

written 8 months ago by Philipp O
You are printing the x, y coordinates of the nodes of Lagrange polynomials. They of course dont change. You evaluate Function as u(0.5,0.5)
written 8 months ago by Praveen C
In every evaluation of function f, I also print the coordinates at which it is evaluated (between "before project" and "before dof coordinates"). The points, it is evaluated at, coincide with the nodes of the Lagrange polynomials (after "before dof coordinates"). If project used quadrature, function f needed to be evaluated at additional quadrature points. Hence, I suspect that project behaves the same as interpolate in this situation.

For the function u obtained solving the variational problem, u( 0.5, 0.5 ) results in the same runtime error as u( [ 0.5, 0.5 ] ).
written 8 months ago by Philipp O
1
It works like this. During projection, F is first interpolated to f which is degree=3, then the projection is performed. Your V is also of degree=3. So what you observe is correct. Try f=F(degree=6), i.e, use a higher degree than that of V.
written 8 months ago by Praveen C
Thank you for the clarification Praveen.
written 8 months ago by Philipp O
Dear Praveen,

I compared the results of interpolation and projection onto function space V:
from dolfin import *
import numpy as np

order = 1

def fun( x ):
return np.sqrt( ( x[ 0 ] - 0.5 ) * ( x[ 0 ] - 0.5 ) + ( x[ 1 ] - 0.5 ) * ( x[ 1 ] - 0.5 ) ) / np.sqrt( 2 )

class F(Expression):
def eval(self, value, x):
value[0] = fun(x)

f = F( degree = order )

mesh   = UnitSquareMesh( 10, 10 )
V      = FunctionSpace( mesh, 'CG', order )

# fun at ( 0.5, 0.5 )
print( "fun: ", fun( [ 0.342, 0.567 ] ) )

# obtained using interpolation
u_int  = interpolate( f, V )
print( "interpolation: ", u_int( 0.342, 0.567 ) )

# obtained using projection on V
u_proj = project( f, V )
print( "projection: ", u_proj( 0.342, 0.567 ) )
This results in output:
('fun: ', 0.12135279148004793)
('interpolation: ', 0.12819751831041681)
('projection: ', 0.17649268233989812)

If I increase order to 3, the output reads:
('fun: ', 0.12135279148004793)
('interpolation: ', 0.12133966940972071)
('projection: ', -23.182892314863192)
If I increase order to 5, the output reads:
('fun: ', 0.12135279148004793)
('interpolation: ', 0.12135257849628993)
('projection: ', -12.725166601775042)
If I also change the mesh to UnitSquareMesh( 100, 100 ), the output reads:
('fun: ', 0.12135279148004793)
('interpolation: ', 0.12135279147997938)
('projection: ', -97.258135982580441)
I do not understand, why the quality of the approximation using projection is decreasing with order and resolution. Do I miss a point here?
written 8 months ago by Philipp O
1
I dont see any problem when I run your code. With order=3, I get

('fun: ', 0.12135279148004793)

('interpolation: ', 0.12133966940972069)

('projection: ', 0.12133966940972071)
written 8 months ago by Praveen C
Thank you Praveen, your help is much appreciated.

The version I obtained the non-convergent results with still contained the command:
parameters["form_compiler"]["quadrature_degree"] = order​

Without this line, I got the same results you did.

written 8 months ago by Philipp O
Not sure if this will help, but I manged to find this really neat bit:

does that help?
written 8 months ago by pf4d

2
8 months ago by
Since Praveen C posted a solution based on project, I wanted to add an additional solution explicitly using quadrature points:
from dolfin import *
import numpy as np

def fun( x ):
return np.sqrt( ( x[ 0 ] - 0.5 ) * ( x[ 0 ] - 0.5 ) + ( x[ 1 ] - 0.5 ) * ( x[ 1 ] - 0.5 ) ) / np.sqrt( 2 )

# order of the CG function space
order = 3

# quadrature needs to be the double to be exact for polynomials up to degree \$order

# build function spaces with quadrature element and CG elements
mesh = UnitSquareMesh( 3, 3 )
Q = FunctionSpace( mesh, QE )
V = FunctionSpace( mesh, 'CG', order )

# trial and test functions are taken from space I want to project on
g = TrialFunction( V )
phi = TestFunction( V )

# fq goes into the RHS and we prescribe function values at quadrature points
fq = Function( Q )
coor = Q.tabulate_dof_coordinates().reshape((-1,2))
coef = np.array( [ fun( p ) for p in coor ] )
fq.vector().set_local( coef )

# Problem ensures Galerkin orthogonality of the result of projection
a = g * phi * dx
L = fq * phi * dx

# solve for a function in the CG function space
u = Function( V )
solve( a == L, u )

# test accuracy of the approach
points = np.random.rand( 20, 2 )
for p in points:
print( "%f\t%f\t%e" % ( p[ 0 ], p[ 1 ], fun( p ) - u( p ) ) )
print( "num of quadrature points: %d" % ( len(V.tabulate_dof_coordinates().reshape((-1,2))) ) )​
Thank you to all that helped by proposing solutions.
2
8 months ago by
Interpolation would be much more easier to do in your case. Here is an example
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

fun = lambda x: np.sin(2*np.pi*x[:,0])*np.sin(2*np.pi*x[:,1])

mesh = UnitSquareMesh(10, 10)
gdim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)
dofs_x = V.tabulate_dof_coordinates().reshape((-1, gdim))
u = Function(V)
u.vector()[:] = fun(dofs_x)
plot(u)
plt.savefig('u.png')​​

So dofs_x[i,:] gives the coordinates of the support point of i'th dof. You can pass this to your Python function.
2
Oh is that what is needed?  you can do this by sub-classing "fenics::Expression()" too:
from dolfin          import *
import numpy             as np
import matplotlib.pyplot as plt

fun = lambda x:  np.sin(2*np.pi*x[0]) \
* np.sin(2*np.pi*x[1])

class F(Expression):
def eval(self, value, x):
value[0] = fun(x)
f = F(degree=2)

mesh = UnitSquareMesh(10, 10)
V    = FunctionSpace(mesh, 'CG', 1)
u    = interpolate(f, V)

plot(u, interactive=True)
plt.savefig('u.png')
written 8 months ago by pf4d
This is a better solution.
written 8 months ago by Praveen C
it does have the advantage of being simpler to implement, but your solution is very illustrative of the underlying representation in memory; nicely done.
written 8 months ago by pf4d
Thank you, this helped me to understand more about interpolation in FEniCS.
written 8 months ago by Philipp O
0
8 months ago by
Not sure about what you are trying to communicate in the second paragraph, but you can evaluate a "fenics::Function()" via the syntax
f   = Function(V)
f_0 = f(x_0, y_0, z_0)​
where $x_0$, $y_0$, and $z_0$ are coordinates in three dimensions.
Sorry, the second paragraph was incorrect.
written 8 months ago by Philipp O