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:
- Find the element that the point lies in. For multiple queries it is best to build some tree structure (e.g using NearestNeighbors.jl)
- Compute the corresponding local coordinates of the point in that element. For non-linear shape functions this requires iteration (the Jacobian is the gradient).
- Create a new quadrature rule from the local point inside the element
- Create a
fe_values = CellScalarValuesorfe_values = CellVectorValuewith the interpolation and the just created quadrature rule. -
reinit!thefe_valueson the cell containing the point - Evaluate
function_valueusingfe_valuesand 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
grid = generate_grid(QuadraticTriangle, (3, 3))
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