### How to evaluate a function on a FiniteElement( family = "Quadrature",... )?

298

views

0

I try to evaluate a function defined on a function space based on finite elements using quadrature:

```
from dolfin import *
mesh = UnitSquareMesh( 1, 1 )
FE = FiniteElement( family = "Quadrature", cell = mesh.ufl_cell(), degree = 1, quad_scheme = "default" )
V = FunctionSpace( mesh, FE )
u = Function( V )
print( u( 0.5, 0.5 ) )
```

results in error message:

```
Traceback (most recent call last):
File "testProjection2.py", line 8, in <module>
print( 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.
Abort trap: 6
```

Is there a way of evaluating function u at arbitrary points within the unit square?

Community: FEniCS Project

### 2 Answers

1

If you use one quadrature point per cell you can just project onto a DG0 space and you should get the same results:

Thank you for your recommendation. Sadly, I need higher-order approximations such as third-order Lagrange elements.

```
from dolfin import *
mesh = UnitSquareMesh( 1, 1 )
FE = FiniteElement( family = "Quadrature", cell = mesh.ufl_cell(), degree = 1, quad_scheme = "default" )
V = FunctionSpace( mesh, FE )
u = Function( V )
FE0 = FiniteElement( 'DG', mesh.ufl_cell(), 0 )
V0 = FunctionSpace( mesh, FE0 )
u0 = project(u, V0)
print( u0( 0.5, 0.5 ) )
```

written
8 months ago by
Philipp O

1

Hi there,

Using the function evaluate_basis you can find the shape functions evaluated at the point you choose, then the interpolation is

$u=\Sigma\left(u_i\cdot\phi_i\right)$

A working example to interpolate the value at a point is:

Good Luck

Using the function evaluate_basis you can find the shape functions evaluated at the point you choose, then the interpolation is

$u=\Sigma\left(u_i\cdot\phi_i\right)$

`u`=Σ(`u`_{i}·`ϕ`_{i})A working example to interpolate the value at a point is:

```
from fenics import *
import numpy as np
mesh = UnitSquareMesh(3, 3)
topology = mesh.cells() #stores mesh connectivity
coordinates = mesh.coordinates() #stores mesh coordinates
V = FunctionSpace(mesh, 'CG', 1)
#your function
u_exp = Expression('x[0]+x[1]',degree=1)
u = interpolate(u_exp, V)
el = V.element()
# Point to be interpolated
x = np.array([0.2, 0.4])
# Find the cell with point
x_point = Point(*x)
cell_id = mesh.bounding_box_tree().compute_first_entity_collision(x_point)
cell = Cell(mesh, cell_id)
coordinate_dofs = cell.get_vertex_coordinates()
cell_topology = topology[cell_id] #stores connectivity of found cell
interpolated_value = 0 #to store sumation
# Array for values. Here it is simply a single scalar
values = np.zeros(1, dtype=float)
for i in range(el.space_dimension()):
#values return the shape function evaluated at node i
el.evaluate_basis(i, values, x, coordinate_dofs, cell.orientation())
node_i = cell_topology[i]
node_xy = coordinates[node_i]
u_node = u(node_xy)
#interpolated_value = summation(u_i * shape_function_i)
interpolated_value += u_node * values
print i, node_i, node_xy, values, u_node
print interpolated_value
```

Good Luck

1

Thank you for the detailed answer. Anyhow, this does not solve the problem as I am looking at finite elements of the family quadrature for which evaluate_basis is not implemented.

written
8 months ago by
Philipp O

1

Yes. Also note that you can perform the equivalent analysis by

```
from fenics import *
mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, 'CG', 1)
u_exp = Expression('x[0]+x[1]',degree=1)
u = interpolate(u_exp, V)
interpolated_value = u(0.2, 0.4)
print interpolated_value
```

written
8 months ago by
pf4d

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

Why don't you create a refined mesh or increase the polynomial order of the non-quadratic element basis --- or both --- and evaluate your function from there?