Evaluate_basis() gives strange results

variational formulation
1

Dear all,

I would like to assemble the right hand side of a variational problem which is the integral over a circle of a source term. More precisely, I would like to solve the problem of finding u such that

\int_{\Omega} \nabla u \cdot \nabla v = \int_C{fv},

where \Omega is a square, C is a circle, f is a scalar source term and v is a test function.

To assemble the right hand side, I discretize the circle as a collection of points and, for each one of them, I find the element containing the point. Then I evaluate the (three) basis functions of this element at this particular point using the evaluate_basis function

However, I guess I am doing something wrong because the results I obtain are not correct.

Here is a minimal working example:

from fenics import *
import numpy as np

Nx = 10
Ny = 10
mesh = UnitSquareMesh(10,10)
tree = mesh.bounding_box_tree()
N = 10
theta = np.linspace(0,2*np.pi,N)
xc = 0.5
yc = 0.5
rc = 0.1

V = FunctionSpace(mesh, "CG", 1)
elements = V.element()
dofmap = V.dofmap()

for i in range(N):
    pt = Point(xc+rc*np.cos(theta[i]), yc + rc*np.sin(theta[i]))
    collision = tree.compute_entity_collisions(pt)
    for col in collision:
        cell = Cell(mesh,col)
        coord = elements.tabulate_dof_coordinates(cell)
        dof = dofmap.cell_dofs(cell.index())
        basis = elements.evaluate_basis(dof[0], np.array([pt[0], pt[1]]), coord, cell.orientation())
        print(basis)

And this is the result of the evaluation:

python mwe_evaluate_basis.py
[-3.79570335e+52]
[-3.79570335e+52]
[6.92943466e-310]
[0.]
[1.e-323]
[1.e-323]
[0.]
[0.]
[2.13436318e-314]
[2.13436318e-314]
[0.]
[0.]
[-3.79570335e+52]
[-3.79570335e+52]
[-3.79570335e+52]
[-3.79570335e+52]
[6.92943452e-310]
[0.]
[1.e-323]
[1.e-323]

Any ideas ?

2

I think what’s going on here is that you’re supposed to pass an element-local index (so, for a CG1 triangle, in the range 0–2) for the first argument to elements.evaluate_basis, but you’re passing global DoF indices. If I change dof[0] to 0 and print separators between points, you can see that you get multiple colliding elements with basis functions evaluating to 1 or (machine) zero at the first and last point (since they line up with mesh vertices), but, at each of the intermediate points, one colliding element is found, and the 0^\text{th} basis function evaluates to something between zero and one, as expected.

1 Like
3

Thank you for your response. It was not clear for me what was expected as first argument but I understand now.
An other solution I have found is to use elements.evaluate_basis_all which evaluates the three basis functions of a given element and returns an array.

1 Like