How can I use "evaluate(x, mapping, component, index_values)"?

General
1

Hello,
I am trying to extract the coordinate of mesh by using Spatialcoordinate.

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([x_min, y_min]), np.array([x_max, y_max])], 
                               [nx-1, ny-1], mesh.CellType.triangle)
X = SpatialCoordinate(domain)```

The document says I have to use evaluate to get the information of these coordinates :

*class* `ufl.geometry.` `SpatialCoordinate` (*domain* )

Bases: [`ufl.geometry.GeometricCellQuantity`](https://fenicsproject.org/olddocs/ufl/1.4.0/ufl.html?highlight=spatialcoordinate%20evaluate#ufl.geometry.GeometricCellQuantity)

UFL geometry representation: The coordinate in a domain.

In the context of expression integration, represents the domain coordinate of each quadrature point.

In the context of expression evaluation in a point, represents the value of that point.

`evaluate` (*x*, *mapping*, *component*, *index_values* )

`is_cellwise_constant` ()

Return whether this expression is spatially constant over each cell.

`name` *= 'x'*

`shape` ()

I also found similar questions, but he dealt with object, whereas mine is spatialcoordinate from mesh.

How can I deal with this?
The reason why I am doing this is that I have to extract the calculated FEM solution and coordinate to numpy data.(There is no problem in Pyvista visualiaztion) Thus, to match the coordinate (x,y) and solution f(x,y) I have to find out the indexes of every point.

Thanks for reading.

2

Dear psj1866,

You can find the coordinates of your mesh using (note to import dolfinx)

num_facets_owned_by_proc = domain.topology.index_map(fdim).size_local
geometry_entities = dolfinx.cpp.mesh.entities_to_geometry(domain, fdim, np.arange(num_facets_owned_by_proc, dtype=np.int32), False)
points = domain.geometry.x
print('e, points[entity]')
for e, entity in enumerate(geometry_entities):
    print(e, points[entity])

To find the coordinates corresponding to your function, it depends on what space your function lies in.
For example you could have as a function space

V = FunctionSpace(domain, ("CG", 1))

Then you can find the corresponding coordinates using

dof_coordV = V.tabulate_dof_coordinates()[:,:]
print(dof_coordV)

I hope this helps.

Best,
Robin

1 Like
3

Dear Robin,
I really appreciate your quick, concise, and wonderful solution. :grinning: I can see my coordinates!
However, how can I find some related document?
Thanks a lot.

1 Like
4

Dear psj1886,

You are very welcome :slight_smile: . I could not find some related document unfortunately. If I recall correctly I found these pieces of code before on a different forum post when I was struggling with the same issue. However I can very much recommend the tutorial written by dokken at Implementation — FEniCSx tutorial . In this first chapter you will also find the code for tabulating the DoF for function spaces

Hope this helps!

Best,
Robin

Edit: I also found this https://fenicsproject.org/pub/tutorial/html/._ftut1019.html which could be somewhat helpful

1 Like
5

Dear Robin,
Thanks for your reply.
Seems that I missed the very basics.^^;; I should have read the documents more carefully.
Again, thanks for spending your valuable time on answer.

1 Like
6

Wrt. this function the C++ docs explains the input/output: Mesh (dolfinx::mesh) — DOLFINx 0.3.1 documentation

As @rsmeets says, for certain function spaces you can get the dof coordinates directly with dolfinx.fem.FunctionSpace.tabulate_dof_coordinates. This is not the case for finite elements whose degree of freedom is based on functionals that does integration and not point eval. Then the dofs are associated with a cell,edge or facet, such as https://defelement.com/elements/examples/triangle-N1curl-1.html.

3 Likes