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

298
views
0
8 months ago by

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
1
the function "evaluate_basis()" will perform a mini-interpolation over only the basis functions associated with the cell which includes the point you want to evaluate.  Since this has not been supported/implemented, you cannot perform this evaluation.  Perhaps there is a trick around this.

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?
written 8 months ago by pf4d

1
8 months ago by
If you use one quadrature point per cell you can just project onto a DG0 space and you should get the same results:
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 ) )​
Thank you for your recommendation. Sadly, I need higher-order approximations such as third-order Lagrange elements.
written 8 months ago by Philipp O
1
8 months ago by
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)$u=Σ(ui·ϕ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)

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