El.evaluate_basis_all() gives different results

1

I’m having a hard time evaluating basis functions here. The following code creates a mesh of a triangle and evaluates the basis functions at a particular point inside the triangle.

import numpy
from dolfin import BoundingBoxTree, Cell, FunctionSpace, Mesh, MeshEditor, Point

import meshzoo


def create_dolfin_mesh(points, cells):
    # https://bitbucket.org/fenics-project/dolfin/issues/845/initialize-mesh-from-vertices
    editor = MeshEditor()
    mesh = Mesh()
    editor.open(mesh, "triangle", 2, 2)
    editor.init_vertices(points.shape[0])
    editor.init_cells(cells.shape[0])
    for k, point in enumerate(points):
        editor.add_vertex(k, point)
    for k, cell in enumerate(cells):
        editor.add_cell(k, cell)
    editor.close()
    return mesh


points, cells = meshzoo.triangle(
    n=8, corners=numpy.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
)
mesh = create_dolfin_mesh(points, cells)
V = FunctionSpace(mesh, "CG", 1)

bbt = BoundingBoxTree()
bbt.build(mesh)
dofmap = V.dofmap()
el = V.element()
rows = []
cols = []
data = []

xy = [0.38764254, 0.29652301]

cell_id = bbt.compute_first_entity_collision(Point(*xy))
cell = Cell(mesh, cell_id)
coordinate_dofs = cell.get_vertex_coordinates()

print(el.signature())
print(xy, coordinate_dofs, cell_id)
v = el.evaluate_basis_all(xy, coordinate_dofs, cell_id)

print(v)

Everything works fine:

FiniteElement('Lagrange', triangle, 1)
[0.38764254 0.29652301] [0.375, 0.25, 0.5, 0.25, 0.375, 0.375] 31
[0.5266756  0.10114032 0.37218408]

Unfortunately, when I execute the same code embedded in a larger code base, I get

FiniteElement('Lagrange', triangle, 1)
[0.38764254 0.29652301] [0.375, 0.25, 0.5, 0.25, 0.375, 0.375] 31
[-0.72790868  0.10114035  1.62676833]

Everything is the same, except that the basis functions are evaluated wrong. They still add up to 1.0 though.

I have no idea what else I could do to debug. Any hints?

2

That’s odd; the only way I can think of to get negative basis function evaluations with linear triangles is to evaluate the shape functions at a parametric point outside of the reference element, but the xy position and coordinate_dofs are the same in both cases. Do you get the same behavior with the following settings?

fem.parameters["form_compiler"]["representation"] = "quadrature"

and/or

fem.parameters["form_compiler"]["representation"] = "tsfc"

Note that the second one requires installing TSFC for FEniCS, which can be done with

pip3 install git+https://github.com/blechta/tsfc.git@2018.1.0
pip3 install git+https://github.com/blechta/COFFEE.git@2018.1.0
pip3 install git+https://github.com/blechta/FInAT.git@2018.1.0
pip3 install singledispatch networkx pulp
3

Thanks @kamensky for the reply! It appears to be a memory issue after all: I found out that the error can be avoided by passing xy.copy(). Wow, yeah.

I’ve opened a bug at https://bitbucket.org/fenics-project/dolfin/issues/1090/evaluate_basis_all-gives-incorrect-result, perhaps that will turn up something.