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

**URL:** <https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975>\
**Category:** Modelling & Simulations\
**Created:** [May 7, 2020, 5:47pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975 "2020-05-07T17:47:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![jbmuir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbmuir/32/4798_2.png) [@jbmuir](https://discourse.julialang.org/u/jbmuir)\
**Post date:** [May 7, 2020, 5:47pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/1 "2020-05-07T17:47:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [May 8, 2020, 9:02am UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2 "2020-05-08T09:02:23Z")

</div>

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:

```julia
using JuAFEM
grid = generate_grid(QuadraticTriangle, (3, 3))

```

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

```julia
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

```julia
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
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.

```julia
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
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
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:

```julia
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

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![jbmuir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbmuir/32/4798_2.png) [@jbmuir](https://discourse.julialang.org/u/jbmuir)\
**Post date:** [May 8, 2020, 3:20pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/3 "2020-05-08T15:20:25Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [December 21, 2020, 1:58pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/4 "2020-12-21T13:58:07Z")

</div>

@jbmuir any progress on this? otherwise I’ll start from Kristoffer’s suggestion

---

<div class="post-metadata">

**Author:** ![jbmuir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbmuir/32/4798_2.png) [@jbmuir](https://discourse.julialang.org/u/jbmuir)\
**Post date:** [December 21, 2020, 9:17pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/5 "2020-12-21T21:17:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [December 21, 2020, 9:35pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/6 "2020-12-21T21:35:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jbmuir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbmuir/32/4798_2.png) [@jbmuir](https://discourse.julialang.org/u/jbmuir)\
**Post date:** [January 8, 2021, 10:03pm UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/7 "2021-01-08T22:03:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kimauth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kimauth/32/25362_2.png) [@kimauth](https://discourse.julialang.org/u/kimauth)\
**Post date:** [October 11, 2021, 9:48am UTC](https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/8 "2021-10-11T09:48:12Z")

</div>

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

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

```

Thanks to everyone who contributed! 🙂