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
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
    error("did not find cell containing point")

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]))

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

Let’s run that:

julia> cell = find_cell_containing_point(grid, point)

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")
        N = JuAFEM.value(interpolation, local_guess)

        global_guess = zero(Vec{dim})
        for j in 1:n_basefuncs
            global_guess += N[j] * cell_coordinates[j]
        residual = global_guess - global_coordinate
        if norm(residual) <= tol_norm
        dNdξ = JuAFEM.derivative(interpolation, local_guess)
        J = zero(Tensor{2, 2})
        for j in 1:n_basefuncs
            J += cell_coordinates[j] ⊗ dNdξ[j]
        local_guess -= inv(J) ⋅ residual
    return local_guess

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}:

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)

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


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


It’s been a long time, but as of today, this is a feature we have in Ferrite.jl! :tada:

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

Thanks to everyone who contributed! :slight_smile: