Can I improve the accurace of node coordinates

dolfinx
1

Hello,

I want to ask a question about the Mesh in DOLFINX.

When I use DOLFINx, I usually find that the DOFs are not located at what they “should be”, which means if I print their coordinates, they will not equal to the analytical results. I understand that there will always be some round off errors, but in some cases it does not seem to be a round off error.

For example, if I create a square, and there are two DOFs in the same column, but there x-coordinates are not always equal, especially when I use high-order polynimials. Does any one know where this error comes from or how could I control it?

Thanks a lot!

2

We do have a dtype argument to be passed in at mesh creation. Look for instance at the signature of dolfinx.mesh.create_unit_cube.

With that being said, our tests (for instance at

) use either np.float32 or np.float64, with 64 being the default for users.
numpy has 96 and 128 as well, see
https://numpy.org/doc/stable/reference/arrays.scalars.html#numpy.float96
but it’s very likely that we don’t test those very well.

3

Note that tabulating dof coordinates has more floating point error accumulation than say the mesh.geometry.x, as it involves tabulation of basis functions.

Secondly, whenever one works with floating values, you should use np.isclose or np.allclose with relevant tolerances as one should never do a==b if a and b are to numpy.float64.

4

Hello, thanks for replying!
I have tested the code a bit, now I can feel that “mesh.geometry.x” gives more accurate dof coordinates.

But I think this only gives the coordinates of all “Nodes” instead of the “Dofs”, right? So when I work with high-order FEM, say for example P3, then the two DoFs inside the element I still need to use tabulating dof coordinates, wihch is not what I want.

Also, when I interpolate the function f(x,y) in DOLFINx, the x and y I got are also given by tabulating dof coordinates, right? Is there any way to get rid of it? Thanks!

5

Yes, the mesh.geometry.x is more accurate when the nodes of the mesh coincide with the degrees of freedom, which is only the case in iso-parametric FEM methods.

As I have already stated, there will always be floating error in tabulation (usually of somewhere around sqrt(floating_precision)), which you cannot avoid. It rarely causes any problems as long as you use double-precision floats.

Yes and no, this depends on your function space.
What happens is that the function is evaluated at the interpolation points of the finite element (which is the dof coordinates for Lagrange, but not for Nedelec and RT), and these are then pushed forward to the physical element for evaluation.
For Nedelec/RT etc then a linear combination and a pull-back routine is used to account for the global orientation vs local orientation.