Getting basis functions and evaluating them at a point

1

Howdy,

So per here, back in 2020 (Accessing the basis function of a Finite Element)

a piece of code was provided to not only get the basis functions, but also evaluate them at a point. I tried to do this a few ways today

from dolfin import *
import numpy as np


a=0.0
b=1.0
n=1000

order_of_polynomial = 1
spatial_mesh= IntervalMesh(n,a,b)
#tree = mesh.bounding_box_tree()
tree = BoundingBoxTree()
tree.build()
V = FunctionSpace(spatial_mesh,"CG",order_of_polynomial)
phi_j=TestFunction(V)
phi_i = TrialFunction(V)


a=phi_i*phi_j*dx

M = PETScMatrix()
assemble(a, tensor=M)

Mmat = M.mat()
print(type(Mmat))

Mcopy = Mmat.duplicate()
Mcopy.zeroEntries()

print(Mcopy.getRow(3))
print(Mmat.getRow(3))
print(M.getrow(3))
print(Mmat)

Mmat_dense = Mmat.convert("dense")

print(Mmat_dense.getDenseArray())

y=0.5
point = np.array([1 / 3 + 1 / 12])
list_of_basis_funcs = []
# Evaluate the i-th basis function at point x:
def basis_func(i,x):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = V.dofmap().cell_dofs(cell_index)
    for local_dof in range(0,len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = array([0,])
            return V.element().evaluate_basis(local_dof,x,
                                              cell.get_vertex_coordinates(),
                                              cell.orientation())
    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Build a Python list of all basis functions, as requested:
for i in range(0,V.dim()):
    list_of_basis_funcs += [(lambda i: lambda x : basis_func(i,x))(i),]

# Try evaluating all basis functions at a point:
for i in range(0,V.dim()):
    print(list_of_basis_funcs[i](point))

This then produces the output

Traceback (most recent call last):
  File "fenicsex.py", line 62, in <module>
    print(list_of_basis_funcs[i](point))
  File "fenicsex.py", line 58, in <lambda>
    list_of_basis_funcs += [(lambda i: lambda x : basis_func(i,x))(i),]
  File "fenicsex.py", line 43, in basis_func
    cell_index = tree.compute_first_entity_collision(Point(*x))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***    
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to compute collisions with bounding box tree.
*** Reason:  Bounding box tree has not been built. You need to call tree.build().
*** Where:   This error was encountered inside BoundingBoxTree.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  
*** -------------------------------------------------------------------------

Traceback (most recent call last):
File “fenicsex.py”, line 13, in
tree.build()
TypeError: build(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.geometry.BoundingBoxTree, arg0: dolfin.cpp.mesh.Mesh) → None
2. (self: dolfin.cpp.geometry.BoundingBoxTree, arg0: dolfin.cpp.mesh.Mesh, arg1: int) → None

Invoked with: <dolfin.cpp.geometry.BoundingBoxTree object at 0x7fe460c4ea70>


It would be very helpful if someone could give me some tips on how to fix this code to work as intended, or another simpler way to get all the basis functions and evaluate them at a point.

Thank you.
2

As the error message says, you need to input the mesh to relate the BoundingBoxTree to the mesh, i.e.

tree.build(spatial_mesh)
1 Like