How to get vertex values in dolfinx?

Hi,

I am transferring fenics to dolfinx but having a trouble now. I was wondering how can I get the array of vertex values coinciding with the nodal order like compute_vertex_values() in fenics? The reason why I need this function is to verify my solutions.

Thanks for your help!

Looking at dolfinx/function.py, compute_point_values seems to be the dolfinx equivalent of compute_vertex_values.

3 Likes

Thanks for your prompt reply! But it seems that compute_point_values still generates the same values as the original ones but in a tuple form. Maybe I need to write a script to fulfill my needs.

Compute point values does the same thing as compute vertex values for a first order nesh. For a second order mesh it computes the values at every node on the mesh geometry.

If your mesh is first order, it should return you the values at each vertex.

3 Likes

Thank you! I figure it out this works and the original result coincides with the nodal order as well.

But I found something very strange, if I read the mesh from the xdmf file, the nodal order of the mesh in dolfinx is actually changed after I comparing with the input meshing file. So I was wondering how could this happen?

Please produce a minimal working code example.

Hi dokken, here is an example as what I mentioned before. In this script, I tried to generate and save a simple mesh with hexahedral elements and then use this mesh file again to check if there is any change.

import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, Constant, BoxMesh, VectorFunctionSpace
from dolfinx.io import XDMFFile

# Create the mesh file
Mesh =BoxMesh(
    MPI.COMM_WORLD, [np.array([-0.5, -5, 0]),
                     np.array([0.5, 5, 0.5])], [2, 2, 2],
    CellType.hexahedron, dolfinx.cpp.mesh.GhostMode.none)
# Save solution in XDMF format
with XDMFFile(MPI.COMM_WORLD, "mesh_hex.xdmf", "w") as file:
    file.write_mesh(Mesh)

def expr(x):
    return x[0] + x[1] + x[2]

# use "Mesh" to generate function u
W = FunctionSpace(Mesh, ("Lagrange", 1))
w = Function(W)
w.interpolate(expr)
vertex_w = w.compute_point_values()
print(vertex_w)    

# use "mesh" to generate function u
with XDMFFile(MPI.COMM_WORLD,"mesh_hex.xdmf", "r") as infile:
    mesh = infile.read_mesh(name = "mesh")
U = FunctionSpace(mesh, ("Lagrange", 1))
u = Function(U)
u.interpolate(expr)
vertex_u = u.compute_point_values()

# compare the difference
print(vertex_u - vertex_w)
print(np.linalg.norm((vertex_u - vertex_w)))


with XDMFFile(MPI.COMM_WORLD, "test_hex.xdmf", "w") as out_file:
    out_file.write_mesh(mesh)
    out_file.write_function(u)

And the difference of the two functions is not zero, which looks weird to me. Thanks in advance for your any comments.

If you write w to file, you will observe that the function is the same as the test_hex.xdmf output.

The reason for the vertex values not matching, is that the mesh geometry is ordered differently in the read-in mesh.
This can be shown by:

for coord, Coord in zip(mesh.geometry.x, Mesh.geometry.x):
    print(coord - Coord)

Thanks. I do observe that the function w and u are exactly the same but in different orders. So I was wondering if there is anyway to generate an array of function u with the same order to w i.e. the read-in mesh, as I expected compute_point_values will do the job but it did not.

If you are not using a built mesh (i.e. loading the original mesh with XDMF) I thing it should Give the same order. You should verify this.

Thanks for your quick reply. Yes I agree with that, but I guess I do need to use a mesh with the help of a XDMF input file. So does that mean I should mesh the geometry directly in dolfinx if I want to keep the same order?

Is there a particular reason for wanting to use multiple meshes of the same mesh in the same file with IO? What is the application?

I do not intend to use multiple meshes. I would like to verify the numerical solution in dolfinx, so I need to use the same mesh file. This works well when I was using fenics, but after transitioning to dolfinx I haven’t figured out how to accomplish that.

So do you want to compare dolfin and dolfinx solutions? Why not work with a manufactured problem to verify the error norms with the exact solutions? Using compute point values is at least not the way to go. If anything, you should tabulate the dof coordinates, and create a map from the dof coordinates in dolfinx to those in dolfin

Thanks dokken. Actually I want to compare with the solution computed in ABAQUS, since there is not a closed form for my problem. But I believe I can do it with “dofmap” as you suggested.

Has the compute_point_values() function been removed? I get the Error 'Function' object has no attribute 'compute_point_values'…

From the change log of the dolfinx-tutorial,

dolfinx.fem.Function.compute_point_values has been deprecated. Interpolation into a CG-1 is now the way of getting vertex values.

2 Likes