Accessing the basis function of a Finite Element

1

Hello to everybody,

I am working with a problem where I need to directly access the basis functions of a Finite Element. A minimal example would be:

>>> from fenics import *
>>> mesh = UnitSquareMesh(8,8)
>>> V = FunctionSpace(mesh, 'P', 1)
>>> V.dim()

displays 81. Then I’d like to have a Python list of all these 81 functions, and being able to evaluate/use them. The precise reason beyond this is very long to describe, since it involves connections with other mathematical objects I am studying, so for simplicity I prefer stopping here.

This question has been already posted on this website (https://fenicsproject.org/qa/7820/direct-finite-element-functions-specified-functionspace/), but the answer finds another solution not enough for my case (I specifically need the list). My FEniCS knowledge is roughly bases on the official book “tutorial 1”, and I have already tried with help(V), from which I learnt about .dim().

Thanks a lot in advance!

2

Is something like the following what you were thinking of?

from dolfin import *
from numpy import array

# Simple mesh/space for testing:
mesh = UnitIntervalMesh(3)
tree = mesh.bounding_box_tree()
V = FunctionSpace(mesh,"CG",1)

# 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:
list_of_basis_funcs = []
for i in range(0,V.dim()):
    list_of_basis_funcs += [(lambda x : basis_func(i,x)),]

# Try evaluating all basis functions at a point:
for i in range(0,V.dim()):
    print(list_of_basis_funcs[i](array([1/3+1/12,])))
2 Likes
3

Hello kamensky, it looks great! I tested a little bit and seems to do precisely what I need. May I ask for some further reference, in order to better understand the code and adapt to other similar cases?

Thanks a lot

1 Like
4

The API for these lower-level manipulations is not as clearly-documented as the main use cases for FEniCS. I mostly found the functions used in the above answer in old forum posts that came up in Google search results. I also stockpile all of the code snippets from my own forum posts in one directory and run grep -i <relevant string> *.py there when I need to remember how to implement something uncommon.

1 Like
5

Dear Kamensky,

thanks for the grep hint, I like it. Coming back to the original problem, there is still something puzzling me. Using precisely your script, where V.dim() is 4, if I evaluate in a for loop like yours:

point = array([1 / 3 + 1 / 12])
for i in range(0, V.dim()):
\tab print(list_of_basis_funcs [i] (point))

I obtain different values (0, 0.25, 0.75 and 0). But if I evaluate separately:

print(list_of_basis_funcs [0] (point))
print(list_of_basis_funcs [1] (point))
print(list_of_basis_funcs [2] (point))
print(list_of_basis_funcs [3] (point)

the results are always zero. How can it be possible?

6

Interesting! It turns out that I misunderstood what happens with i in the lambda expression. The anonymous function created will take whatever is the value of i in the global scope when it is called, not the one at which it was defined. One trick I found here is to do the following:

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

Looking at this code again, you might also add a [0] at the end of the eval_basis line, to get rid of the array/scalar discrepancy.

2 Likes
7

Great! Thanks so much: not only it works very well, but I took the opportunity for learning this lambda trick which is definitely good to know.

At this point I am wondering if I could go a little bit further: is it maybe possible to obtain a symbolic expression of these functions?

Therefore a feature like: given a finite element F, obtain the list of all the polynomials spanning the space. Thanks a lot in advance!

8

Hi,

This implementation is very interesting! Do you know if it’s possible to obtain the second derivative of the basis functions over a vector field?

Many thanks!

9

How can access basis functions, apply this method in new versions of Fenicsx and Dolfinx ?

10

Consider the following in DOLFINx:

V = fem.functionspace(msh,  ("Lagrange", 1))
basix_element = V.element.basix_element
basis_matrix = basix_element.tabulate(0, x_refs)[0, :, :, 0]

V.element.basix_element gives you a basix_element class, and then you can use it to evaluate the value of basis function with tabulate at specific points x_refs. Rember to pull_back first.

11

Thank you very much for the answer.

In this way, I can access the basis function of the reference triangle (in basix, the element used to create the triangle). How can I access the basis function at any point of a cell (triangle) in any arbitrary geometry?

As an example, for x = np.array([[0.0, 0.0], [0.0,1.0], [1.0, 0.0]]) and x = np.array([[0.2, 0.4], [0.4, 0.9], [0.8 , 0.3]]), I get the same basis function values for mesh/cell/triangles after this process, because it tabulates basix_element = V.element.basix_element.

How can I access the basis functions of different geometries? Theoretically, I know that I need to use the Jacobian matrix and how it is obtained theoretically, and when I create a function space with Fenicsx, I see that a function space is created by adhering to this difference. I want to ask a function or convenient method to reach the basis function values at any point for different triangles?

Thank you very much again.

12

You need to find the reference coordinate first by using pull_back function. See the following for the basis function evaluation process:

and the practical implmentation:

Note it is subject to change depending on the version of DOLFINx you use.

13

If you know what point in the reference element you want to evaluate your function at, you can use dolfinx.fem.Expression, consider the mwe

from mpi4py import MPI

import dolfinx
import numpy as np
import ufl

def eval_trial_function(V, reference_points):
    mesh = V.mesh
    u = ufl.TrialFunction(V)
    expr = dolfinx.fem.Expression(u, reference_points)

    num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
    basis_functions = expr.eval(mesh, np.arange(num_cells, dtype=np.int32)).reshape(num_cells, len(reference_points), -1)
    for i, functions in enumerate(basis_functions):
        for j, point in enumerate(reference_points):    
            print(f"Cell {i}: Point {point} {functions[j, :]}")


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, dolfinx.cpp.mesh.CellType.triangle)
V = dolfinx.fem.functionspace(mesh, ("N1curl", 1))

# Points in reference element
points = np.array([[0.0, 0.0], [1.0, 0.0]])
eval_trial_function(V, points)

yielding:

Cell 0: Point [0. 0.] [ 0.  0.  1.  0.  1. -1.]
Cell 0: Point [1. 0.] [0. 0. 1. 1. 0. 0.]
Cell 1: Point [0. 0.] [-0.  1. -1. -0.  0.  1.]
Cell 1: Point [1. 0.] [-1.  0.  0. -0.  0.  1.]
14

Thank you very much for your answer.

If any point of any arbitrary element I created with this code is mapped to the master element, I correctly find the basis function value of the resulting point (with this code, I correctly find the basis function value of the reference element location after mapping). But when I enter that point directly (that is, in my new arbitrary mesh, location on the arbitrary mesh, not mapped to the reference element), I want to find the basis function value at that point. So if I give an example like this:

For a triangle mesh consisting of only one element, I defined a triangle consisting of points (0,0),(0,0.5),(0.5,0). I want to reach the N1E basis function value at the point (0.1, 0.15) of this element. For example, for the 2nd edge of the triangle (DefElement)

Theoretically, at this point, I am expecting to reach the basis function vector values of Vx = 1.4 and Vy = 0.4 (from the Jacobian transform ), however I reach the basis function values Vx = 1.7 and Vy = 0.2. These values are exactly the expecting basis function values of the coordinates (0.5 and 15/200 (coordinate/2)) to which these points in the reference element (master triangle) correspond, as you said. These are the baiss function values of refernece elements(after mapping of any arbitrary element to reference elements, these input points are the points of reference triangle after mapping)

Is it possible to give real coordinates (any arbitrary triangle coordinates) as input here and get its basis function values? How can I write the code for this? My question is probably very simple but I don’t fully understand it. I’m new to Fenics. I’d be glad if you could help me. Actually, I think my question about to asking a function for isoparametric mapping.

Thank you very much again

15

I had a similar question; given any location for arbirary geometry (not the reference element), how can I reach the basis function value at that point (for arbitrary mesh). Also, similarly, how can the coordinate and basis function isoparametric mapping be established between the reference element and any arbitrary element (for example, for a triangle)?

16

See Point Sources (Redux) - #6 by dokken where the code does this in three steps:

  1. find what cell a physical coordinate is in
  2. pull that coordinate back to reference
  3. create expression on reference element.
  4. eval for said cell