Evaluate_basis_all in dolfinx

dolfinx
1

Hi everyone ! I used to use the evaluate_basis_all function in dolfin to evaluate basis functions at some points but I can not find it in dolfinx. Can someone tell me where to find it ?

Cheers,
Fabien

2

This function does no longer exist.
Your intended usage is to evalute the basis functions at a point of physical coordinates?

3

I need to evaluate basis functions at physical coordinates inside a cell.

4

I have written a C++ function for evaluating basis functions at a point in a given cell in my dolfinx add-on dolfinx_mpc.

1 Like
5

Ok, it looks like what I want. How do I use it with the python interface ?

6

The installation instructions are written in the repos readme
and is:

rm -rf build
mkdir -p build
cd build
cmake -G Ninja -DCMAKE_BUILD_TYPE=Developer ../cpp/
ninja -j3 install
cd ../python
pip3 install . --upgrade

For use-cases of the function, see: https://github.com/jorgensd/dolfinx_mpc/blob/853b4a62bc23aa2637a34dd4ca7408d08dfadce6/python/dolfinx_mpc/contactcondition.py#L115-L116

7

It works fine thanks ! It would be nice to have this functionality directly in doflinx.

8

You can do all the operations in my function in the dolfinx Python layer, it might just suffer a bit on the performance side.

9

Yes of course. I am just wondering why there is no function for direct basis functions evaluation in dolfinx, just like the evaluate_basis_all in dolfin. Are there any plans to add this feature in dolfinx ?

1 Like
10

Hello dokken, do you have something similar for de derivative of de basis functions?? Or is there any change I can implement to you code to compute the derivatives. I’m not a C++ programmer but I would try.

11

The DOLFINx way to do this would be to:

  1. Figure out which cell the point x is in (using dolfinx.compute_colliding_cells, as done in: Implementation — FEniCSx tutorial)
  2. Pull these coordinates back to the reference element, using mesh.geometry.cmap.pull_back (dolfinx/fem.cpp at 48ece0087588f5f274a36aec223abb129fa59e8e · FEniCS/dolfinx · GitHub)
  3. Use basix.FiniteElement(element_of choice).tabulate(1, X_ref) to tabulate at the reference points.

An example of the last step can be found in: asimov-custom-assemblers/custom_assembler.py at main · Wells-Group/asimov-custom-assemblers · GitHub
and for a more verbose and detailed demo on basix see: Creating and tabulating an element — Basix 0.4.1.0 documentation

2 Likes
12

thank you very much!!

13

One more question: is there a transformation needed??
something like: \nabla_x\phi(x)=(\nabla_yT)^{-T}\nabla_y\hat{\phi}(y) with y\in\Omega_{ref}.
or is it included in the function??

14

It is not included in the function, as you only evaluate on the reference element.
As it is trivial to compute the Jacobian once you have the point on the reference element (for affine geometries you can compute it at any point of the reference element, but for quadrilaterals/hexahedra, you would need to do it at each point in the reference geometry).
A python code doing this can be found at: asimov-custom-assemblers/kernels.py at 6086833d160afbbe00823ca00cc96abca1b2ea8b · Wells-Group/asimov-custom-assemblers · GitHub
where c_tab is the output from the coordinate-element’s tabulate at these points.

What is the particular application where you need the gradient of test functions?

15

Thank you. I guess I’ll have to compute de transformation.
I’m trying to compute a vector g for the RHS that looks like:
g_k=\nabla_x \phi_k|_{x_{ref}}\cdot\vec{n}=\nabla_y\phi_k|_{x_{ref}} (\frac{dx}{dy})^{-1}\cdot \vec{n}

16

What is the relation between y and x here? In your previous post it seemed like you wanted the gradient of a test function on the physical cell, Which can be computed using ufl

17

Totally forgot to specify
x are the coordinates of the point and y are the coordinates of the point in the ref. domain so that:
y\in \Omega_{ref} & x\in\Omega. The therm (\frac{dx}{dy}) would be the Jacobian of the transformation.

18

My question is then, how is this different from

g = ufl.inner(ufl.grad(ufl.TestFunction(V)), ufl.FacetNormal(mesh))*ufl.dx)
g_k = dolfinx.fem.assemble_vector(dolfinx.fem.form(g))
19

Because no integration is made. The complete RHS looks as follows:
RHS=\int v \;dx \nabla u|_{x_{ref} }\cdot \vec{n}
With u as a trial and v as a test function.
The first therm is a “normal” load vector and the second is the term I’m trying to compute.

20

I think you can achieve what you want with the following:

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)
tdim = domain.topology.dim
V = dolfinx.fem.FunctionSpace(domain, ("CG", 2))

point = np.array([[0.2, 0.3, 0]])
bb = dolfinx.geometry.BoundingBoxTree(domain, tdim)
bbox_collisions = dolfinx.geometry.compute_collisions(bb, point)
cells = dolfinx.geometry.compute_colliding_cells(domain, bbox_collisions, point)
cell = cells.links(0)[0]

cmap = domain.geometry.cmap
geom_dofs = domain.geometry.dofmap.links(cell)
x_ref = cmap.pull_back(point, domain.geometry.x[geom_dofs])

n = ufl.as_vector([0.1, 0.2])
expr = ufl.inner(ufl.grad(ufl.TestFunction(V)), n)
E = dolfinx.fem.Expression(expr, x_ref)
grad_u_at_x_ref = E.eval([cell])
print(f"Cell: {cell}, Point: {point}, x_ref: {x_ref}, grad(u)*n {grad_u_at_x_ref}")
3 Likes