Assemble with dolfin VS dolfinx

dolfinx
1

Hello everyone,

I can’t figure out why I have such big difference between dolfin VS dolfinx with the assemble function.
I am using RT elements. (I have not the problem with CG or DG elements).

Here is a code snippet to reproduce my problem:

import numpy as np
import basix

# Using Dolfin
import dolfin as dlf
mesh_dlf = dlf.BoxMesh(dlf.Point(0.0, 0.0, 0.0),
                dlf.Point(1.0, 1.0, 1.0),
                10,10,10)
V_dlf = dlf.FunctionSpace(mesh_dlf, "RT", 1)
f_dlf = dlf.Function(V_dlf)
f_dlf.vector().set_local(np.ones(V_dlf.dim()))
vol_dlf = dlf.assemble(dlf.Constant(1)  * dlf.Measure('dx', domain=mesh_dlf))
res_dlf = dlf.assemble(dlf.dot(f_dlf,f_dlf)**(1/2) * dlf.dx)
print(vol_dlf)
print(res_dlf)


# Using Dolfinx
import dolfinx as dlfx
from mpi4py import MPI
import ufl
mesh_dlfx = dlfx.mesh.create_unit_cube(MPI.COMM_WORLD,
                                  10,
                                  10,
                                  10,
                                  dlfx.mesh.CellType.tetrahedron)
V_dlfx = dlfx.fem.functionspace(mesh_dlfx,("RT",1))
f_dlfx = dlfx.fem.Function(V_dlfx)
f_dlfx.vector.setArray(np.ones(f_dlfx.vector.size))
vol_dlfx = dlfx.fem.assemble_scalar(dlfx.fem.form(dlfx.fem.Constant(mesh_dlfx, 1.0) * ufl.dx))
res_dlfx = dlfx.fem.assemble_scalar(dlfx.fem.form(ufl.dot(f_dlfx,f_dlfx)**(1/2) * ufl.dx))
print(vol_dlfx)
print(res_dlfx)

I first compute the volume of the domain just to check if everything is ok (and it is indeed: vol_dlf ~ vol_dlfx ~ 1).
But I have a big difference when comparing the norm of the function on the domain:
res_dlf = 141.42
res_dlfx = 312.47

Is it related to different lagrangevariant used ? The default lagrangevariant for RT with dolfinx is “legendre”, but I can’t find the one used with dolfin.

Thanks :slight_smile:

2

This is because setting dof values to a RT array doesn’t really make sense.
You should be quite careful with this, as the degrees of freedom relate to vector basis functions.

Here is a way to illustrate the same, but with an actually well defined function (the constant vector (2,1,3), you could change this to (1,1,1) if you want to)

# # Using Dolfinx

try:
    import dolfinx as dlfx
    from mpi4py import MPI
    import ufl
    import numpy as np

    def f(x):
        values = np.zeros((3, x.shape[1]))
        values[0] = 2
        values[1] = 1
        values[2] = 3
        return values

    mesh_dlfx = dlfx.mesh.create_unit_cube(MPI.COMM_WORLD,
                                           10,
                                           10,
                                           10,
                                           dlfx.mesh.CellType.tetrahedron)
    V_dlfx = dlfx.fem.functionspace(mesh_dlfx, ("RT", 1))
    f_dlfx = dlfx.fem.Function(V_dlfx)
    f_dlfx.interpolate(f)
    vol_dlfx = dlfx.fem.assemble_scalar(dlfx.fem.form(
        dlfx.fem.Constant(mesh_dlfx, 1.0) * ufl.dx))
    res_dlfx = dlfx.fem.assemble_scalar(
        dlfx.fem.form(ufl.dot(f_dlfx, f_dlfx)**(1/2) * ufl.dx))
    print("dolfinx")
    print(f"{vol_dlfx=}")
    print(f"{res_dlfx=}")


except ImportError:
    # Using Dolfin
    import dolfin as dlf
    mesh_dlf = dlf.BoxMesh(dlf.Point(0.0, 0.0, 0.0),
                           dlf.Point(1.0, 1.0, 1.0),
                           10, 10, 10)
    V_dlf = dlf.FunctionSpace(mesh_dlf, "RT", 1)
    f_dlf = dlf.Function(V_dlf)
    f_dlf.interpolate(dlf.Expression(("2", "1", "3"), degree=1))
    vol_dlf = dlf.assemble(dlf.Constant(
        1) * dlf.Measure('dx', domain=mesh_dlf))
    res_dlf = dlf.assemble(dlf.dot(f_dlf, f_dlf)**(1/2) * dlf.dx)
    print("legacy dolfin")
    print(f"{vol_dlf=}")
    print(f"{res_dlf=}")

yields

dolfinx
vol_dlfx=0.9999999999999232
res_dlfx=3.7416573867708562

and

legacy dolfin
vol_dlf=0.9999999999999232
res_dlf=3.74165738677434
1 Like
3

Thank you Dokken for your quick reply !

Then, how can I evaluate a RT basis function ?
I used to do :

V = dlfx.functionspace(mesh,("RT",1))
f = dlfx.Function(V)
values = np.zeros(f.vector.size)
values[5] = 1
f.vector.setArray(values)

to evaluate the 5th basis function for example

4

What do you want to use the basis function for?

I’ve covered how to compute the basis function (and push forward action) for N1curl elements at: Advanced finite elements — FEniCS Workshop
where the covariant piola map is used to map the basis to a given physical cell.

For RT you would need to contravariant piola map.

Note that all of this can be done with

expr = dolfinx.fem.Expression(ufl.TestFunction(V), insert_points_on_reference_element_here)
# select what cell to evaluate in
cells = np.array([1], dtype=np.int32)
print(expr.eval(mesh, cells))

which will give you the basis functions in physical space evaluated at the insert_points_on_reference_element_here pushed forward to the physical cell.

1 Like
5

Interresting documentation, thank’s !

Actually, I want to evaluate an expression for each dof one by one, in order to build a vector G of len=n_dofs, where the ith value correspond to the eval of the expression for the ith dof (dof_i = 1, other_dof = 0)

So I was doing this:

def eval_my_fun(dof_i, vec, f, fxyz):
    vec[dof_i] = 1
    f.vector.setArray(vec)  
    
    # Integrate over all the domain
    dx = ufl.Measure("dx",
                     domain=mesh,
                     metadata={"quadrature_degree": 8})
    res = dlfx.fem.assemble_scalar(dlfx.fem.form((ufl.cross(fxyz,f)[2] / ufl.dot(fxyz,fxyz)**(3/2)) * dx))
    vec[dof_i] = 0
    return res

V = dlfx.fem.functionspace(mesh,("RT",1))
f = dlfx.fem.Function(V)
n_dofs = f.vector.size
u_e = u_expression(c)
f_xyz = dlfx.fem.Function(V)
f_xyz.interpolate(u_e.eval)
vec = np.zeros(n_dofs)
for idof in range(n_dofs):
    G[idof] = eval_my_fun(idof, vec, f, f_xyz)

But since I’ve noticed that there is difference with the assemble function dolfin/dolfinx with RT element, I have some doubts about my way of doing this…

6

Why do you want to do this?
Isn’t this exactly what integrating against a test-function would give you?
If not, please elaborate why?

7

I want to build a Biot&Savart operator, according to:
image

where

  • i is the indice of the magnetic sensor
  • wf_j is the jth face shape function RT (dof_j= 1, other =0)
8

So how would this not be stacking for each i the assembly with a TestFunction? As that would give you a vector with all the contributions were the jth entry of the vector correponds to the contribution of the jth basis function?

9

Sorry, I have difficulties to understand your last reply… Could you show me a code snippet of what you are explaining ? thank you so much for your time

10

Say you want to assemble
L_i=\int_{\Omega} \phi_i(x)\cdot u~\mathrm{d}x
where u is some known vector function.
Then this is exactly what

from mpi4py import MPI
import dolfinx
import ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("RT", 1))

v = ufl.TestFunction(V)
u = dolfinx.fem.Constant(mesh, ((1.0, -3.0)))


L = ufl.inner(u, v) * ufl.dx

b = dolfinx.fem.assemble_vector(dolfinx.fem.form(L))
print(b.array)

Does.

Then for you example, you would vary u for for each r_i to get your matrix operator.

1 Like