# Finding the value of a field at a spatial location in JuAFEM

Is there a relatively simple way to extract the value of a field in JuAFEM at a point in space within the grid? Use cases are - visualisation within Julia (not using ParaView); doing inverse problems with point observations

There is no â€śplug and playâ€ť API to use for this. The steps that are needed are the following:

1. Find the element that the point lies in. For multiple queries it is best to build some tree structure (e.g using NearestNeighbors.jl)
2. Compute the corresponding local coordinates of the point in that element. For non-linear shape functions this requires iteration (the Jacobian is the gradient).
3. Create a new quadrature rule from the local point inside the element
4. Create a `fe_values = CellScalarValues` or `fe_values = CellVectorValue` with the interpolation and the just created quadrature rule.
5. `reinit!` the `fe_values` on the cell containing the point
6. Evaluate `function_value` using `fe_values` and the values of the dofs for the field on the cell.

Here is an example but note that the code below is not very pretty. It would be nice to incorporate this into JuAFEM and add a nice API for it:

Letâ€™s say we have the following grid of second-order triangles:

``````using JuAFEM
``````

we want to find the value of some field in say [0.3, 0.3]:

``````point = Vec(0.3, 0.3);
``````

Step 1 is then to find the element that the point lies in. Below is some code that just loops over all element and checks if the point is inside the triangle

``````function find_cell_containing_point(grid, point)
interpolation = JuAFEM.default_interpolation(eltype(grid.cells))
for cell in 1:length(grid.cells)
cell_coords  = getcoordinates(grid, cell)
# [1,2,3] is the indices for the vertices of a quadratic triangle
coords_vertices = cell_coords[[1,2,3]]
if is_point_inside_triangle(point, coords_vertices)
return cell
end
end
error("did not find cell containing point")
end

function is_point_inside_triangle(point, triangle)
# Just transform thing so they fit the "API" for the function from stack overflow
return _is_point_inside_triangle((x = point[1], y = point[2]),
(x = triangle[1][1], y = triangle[1][2]),
(x = triangle[2][1], y = triangle[2][2]),
(x = triangle[3][1], y = triangle[3][2]))
end

# https://stackoverflow.com/a/34093754
function _is_point_inside_triangle(p, p0, p1, p2)
dX = p.x-p2.x
dY = p.y-p2.y
dX21 = p2.x-p1.x
dY12 = p1.y-p2.y
D = dY12*(p0.x-p2.x) + dX21*(p0.y-p2.y)
s = dY12*dX + dX21*dY
t = (p2.y-p0.y)*dX + (p0.x-p2.x)*dY
(D<0) && return s<=0 && t<=0 && s+t>=D
return s>=0 && t>=0 && s+t<=D
end
``````

Letâ€™s run that:

``````julia> cell = find_cell_containing_point(grid, point)
10
``````

So cell number 10 contains the point.

Now we need to find what local coordinates the point corresponds. The code below implements a small Newton-Raphson method to find it. The first argument is the interpolation we use for our field.

``````function find_local_coordinate(interpolation, cell_coordinates, global_coordinate)
dim = length(global_coordinate)
local_guess = zero(Vec{dim})
n_basefuncs = getnbasefunctions(interpolation)
max_iters = 10
tol_norm = 1e-10
for iter in 1:10
if iter == max_iters
error("did not find a local coordinate")
end
N = JuAFEM.value(interpolation, local_guess)

global_guess = zero(Vec{dim})
for j in 1:n_basefuncs
global_guess += N[j] * cell_coordinates[j]
end
residual = global_guess - global_coordinate
if norm(residual) <= tol_norm
break
end
dNdÎľ = JuAFEM.derivative(interpolation, local_guess)
J = zero(Tensor{2, 2})
for j in 1:n_basefuncs
J += cell_coordinates[j] âŠ— dNdÎľ[j]
end
local_guess -= inv(J) â‹… residual
end
return local_guess
end
``````

Letâ€™s run it:

``````julia> cell_coords = getcoordinates(grid, cell);

julia> interpolation = Lagrange{2, RefTetrahedron, 2}(); # second order interpolation for the field as well

julia> local_coordinate = find_local_coordinate(interpolation, cell_coords, point)
2-element Tensor{1,2,Float64,2}:
0.049999999999999926
0.8999999999999996
``````

So the local coordinates for the point is `[0.05, 0.9]`.

Now we create a quadrature rule with a single point (the local coordinate just computed) and the weight 1.0 and the â€ś`FEValues`â€ť:

``````julia> qr = QuadratureRule{2, RefTetrahedron, Float64}([1], [local_coordinate]);

julia> fe_values = CellScalarValues(qr, interpolation)
CellScalarValues{2,Float64,RefTetrahedron} with 6 shape functions and 1 quadrature points
``````

Now, standard `reinit!` on the cell:

``````reinit!(fe_values, cell_coords)
``````

In a real FE simulation we would have a vector of all the dof values and would extract the dof values for this particular cell now. Here, I am just going to put some random values for each dof

``````dof_values = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]; # 6 dof
``````

and now in standard way we can compute the" function value" at any quadrature point in `fe_values` (here we only have one):

``````julia> q_point;

julia> function_value(fe_values, q_point, dof_values)
2.9400000000000044
``````

So the value in the global point [0.3, 0.3] is 2.94.

6 Likes

Thanks @kristoffer.carlsson for a very detailed response! I guess that points 1-3 at least could be precomputed and put in some struct (say `PointEvalHandler` or something like that) that you could make a functor so that calling `peh(u)` where `peh` is the struct and `u` is a vector of dof gives you the values at the points you want - I will try doing something like that and will see if I can make it somewhat efficient.

1 Like

@jbmuir any progress on this? otherwise Iâ€™ll start from Kristofferâ€™s suggestion

I did make some progress but not quite enough to feel it was worth a pull request; maybe I will try to get it up on the github as a WIP and you can look at it

Sounds good. I have some free time and I think this is a good contribution to JuAFEM

It should now be up as a draft pull request - sorry about a bit of a delay over the end-of-year period

2 Likes

Itâ€™s been a long time, but as of today, this is a feature we have in `Ferrite.jl`!

``````ph = PointEvalHandler(grid, points)
point_values = get_point_values(ph, dh, dof_values)
``````

Thanks to everyone who contributed!

4 Likes