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

407

views

0

I have a python function $f\left(\vec{x}\right)$

[Edit] Removed incorrect part. For clarification of the question, please see below.

`ƒ`(→`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.[Edit] Removed incorrect part. For clarification of the question, please see below.

Community: FEniCS Project

### 3 Answers

2

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
quad_degree = 2 * order
parameters["form_compiler"]["quadrature_degree"] = quad_degree
# build function spaces with quadrature element and CG elements
mesh = UnitSquareMesh( 3, 3 )
QE = FiniteElement( family = 'Quadrature', cell = mesh.ufl_cell(), degree = quad_degree, quad_scheme = "default" )
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

Interpolation would be much more easier to do in your case. Here is an example

So

```
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

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

Please login to add an answer/comment or follow this question.

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

Basically, I try to compute a function $g$

gwithin function space V for which the difference $f-g$ƒ−gis 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)d→x=∫_{[0,1]2}g(→x)Φ(→x)d→xfor 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.Now u-f is orthogonal to V.

This results in output:

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:

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

What seems to work is the following, actually solving the variational problem:

While u can be plotted using

function u cannot be evaluated using:

as this results in error message:

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

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 ] ).

I compared the results of interpolation and projection onto function space V:

This results in output:

If I increase order to 3, the output reads:

If I increase order to 5, the output reads:

If I also change the mesh to UnitSquareMesh( 100, 100 ), the output reads:

I do not understand, why the quality of the approximation using projection is decreasing with order and resolution. Do I miss a point here?

('fun: ', 0.12135279148004793)

('interpolation: ', 0.12133966940972069)

('projection: ', 0.12133966940972071)

The version I obtained the non-convergent results with still contained the command:

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

https://fenicsproject.org/qa/1425/derivatives-at-the-quadrature-points

does that help?