Properly calling evaluate basis derivatives all() function in python

1

Howdy,

Per the documentation here: FiniteElement — FEniCS Project

`evaluate_basis_derivatives_all` ()

Evaluate order n derivatives of all basis functions at given point in cell

Parameters: * **int n** (*unsigned*) –
* *** values** (*double*) –
* **double * x** (*const*) –
* **double * coordinate_dofs** (*const*) –
* **cell_orientation** ([*int*](https://docs.python.org/2/library/functions.html#int)) –
Return type: void

I wrote my code here, assuming since the function returned type void that you would pass a values array of length V.dim() to be filled - this doesn’t look like it works.

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(spatial_mesh)
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([.9])
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(spatial_mesh, cell_index)
            values = np.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


def basis_func_derv(x,n):
	cell_index = tree.compute_first_entity_collision(Point(*x))

	cell = Cell(spatial_mesh, cell_index)
	values = np.zeros(V.dim())
	 

	V.element().evaluate_basis_derivatives_all(n,values,x,
	                                      cell.get_vertex_coordinates(),
	                                      cell.orientation())
	return values
# Build a Python list of all basis functions, as requested:
order=1
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()):
    
	if(list_of_basis_funcs[i](point)!=0.0):
		print((i,list_of_basis_funcs[i](point)[0]))

print(basis_func_derv(point,0))

# println(V.element().evaluate_basis_all())

So a few questions, I effectively followed my evaluate_basis code that looks like it is working. What am I getting wrong here? Also I’d assume that if you call this function and pass the order of the derivative to be 0, then it would just return the basis functions themselves at the point i.e the same as just calling evaluate_basis_all()?

This is the error of the output when I try to run my code:

Traceback (most recent call last):
  File "fenicsex.py", line 80, in <module>
    print(basis_func_derv(point,1))
  File "fenicsex.py", line 67, in basis_func_derv
    cell.orientation())
TypeError: evaluate_basis_derivatives_all(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfin.cpp.fem.FiniteElement, arg0: int, arg1: numpy.ndarray[numpy.float64], arg2: numpy.ndarray[numpy.float64], arg3: int) -> numpy.ndarray[numpy.float64]

Invoked with: <dolfin.cpp.fem.FiniteElement object at 0x7fa5a2a0f1f0>, 1, array([0., 0., 0., ..., 0., 0., 0.]), array([0.9]), [0.899, 0.9], 0
2

Have a look at Parallel Sparse PETSc Matrix Generation with dofs as rows and columns
which calls the function in question.

Additionally, this function: Bitbucket is the glue between the Python and C++ function

3

Thanks @dokken

Just for some else if they get stuck. Evaluate all only returns the values at non zeros. So instead I interate through all using evaluate_basis_derivates, returning 0 if we are not in the support.


def basis_func_derv_i(i,x,n):
    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(spatial_mesh, cell_index)
            values = np.array([0,])
            return V.element().evaluate_basis_derivatives(local_dof,n,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

order =1
first_derv = lambda i,x:basis_func_derv_i(i,x,order)

for i in range(0,V.dim()):
    list_of_basis_functions_der += [(lambda i: lambda x : first_derv(i,x))(i),]

# Try evaluating all basis functions at a point:
for i in range(0,V.dim()):
    
	if(list_of_basis_functions_der[i](point)!=0.0):
		print((i,list_of_basis_functions_der[i](point)[0]))
1 Like