Evaluate PDE solution

[jw3126]

Thanks a lot for creating Gridap! Its awesome to have such a PDE package in pure julia. So say I solved a PDE using Gridap, how do I actually evaluate the solution?

# solve the trivial PDE
# u = sin

using Gridap
model = CartesianDiscreteModel((0,1),10)
f(pt) = sin(pt[1])
V0 = TestFESpace(
  reffe=:Lagrangian, order=2, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="boundary")
U = TrialFESpace(V0)
trian = Triangulation(model)
quad = CellQuadrature(trian, 2)
A(u, v) = u ⊙ v
b(v) = v*f
t_Ω = AffineFETerm(A,b,trian,quad)
op = AffineFEOperator(U,V0,t_Ω)
u = solve(op)

I tried the following, both of which did not work:

pt = [0.1]
u(pt)
evaluate(u, pt)

[zjwegert]

I would also really like to know the answer to this problem.

[santiagobadia]

What you are asking is possible, but not standard in finite element methods. Usually, in finite element methods, you evaluate the solution in the mesh nodes, and you would do it like this:

get_cell_values(u) # nodal values for scalar Lagrangian at each cell

You could also evaluate them in a given integration quadrature

Qₕ = CellQuadrature(trian,2)
q = get_coordinates(Qₕ)
evaluate(u,q)

or you could choose to evaluate it in the point you choose, e.g., in the middle of each cell (element)

ps = [Point(0.5,)] # An array of points with one point
cps = Fill(ps,num_cells(trian)) # A cell-wise array with one point per cell (constant array for performance)
evaluate(u,cps)

or just write the nodal values in a file that can be open by a visualizer (e.g., paraview)

writevtk(trian,"my solution",cellfields=["u" => u])

Why? (not a short answer)

Finite element functions are piecewise functions. Thus, if you just give a coordinate in your domain, one needs to find in which cell of the mesh this point is located. This is easy for some simple meshes but a non-trivial task for general unstructured meshes.

On the other hand, for performance issues, finite element spaces are defined in a reference space, e.g., a unit d-cube. Thus, if you provide a coordinate in the physical space, one needs to push it back to the reference space. It involves to compute roots of polynomials. There is a unique solution for this problem (for acceptable meshes, no negative Jacobians) but it is something that is not really needed to solve finite element problems. You can compute solutions without this inverse geometrical map. In fact, we have not implemented this in Gridap.

In your test you are using a reference space based element. There is also the possibility in Gridap to use finite element spaces in the physical space, thus working right with physical coordinates but incurring in more cpu cost.

What is more standard in finite elements is to provide an array of points in the reference space (or physical space, depending on where they are defined) and return the cell-wise values of the finite element function in these points (as I did above).

Anyway, we could provide an algorithm that given a global arrays of points in the physical space, it finds the cells containing them, then solves roots of polynomials, and then we use the standard machinery. In the case of Cartesian meshes it is straightforward and could be easily done. In high-order (curved) unstructured meshes, it will be more involved.

By the way, pointwise evaluations of many finite element spaces (discontinuous) are not well-defined.

[zjwegert]

Thank you for this fantastic answer. It might be worth mentioning this in the tutorial/documentation - unless I missed it of course.

[santiagobadia]

Yes, we can do it and add, e.g., one case in which the user defines the array of points.

@fverdugo In order to provide a more general evaluation of FE functions, we could do:

Create methods that return the mesh nodes in the physical space and the reference space.

The result is a new CellPoints struct, an array of arrays of points with a trait DefSpace that can be Reference or Physical.

This way:

• ReferenceFE{Physical} and CellPoints{Physical} we can readily evaluate.

• ReferenceFE{Parametric} and CellPoints{Parametric} we can readily evaluate.

• ReferenceFE{Physical} and CellPoints{Parametric} we map the points to the physical space using the geometrical map

function map_points_to_physical_space(CellPoints{Parametric}

# apply standard map

end

• ReferenceFE{Parametric} and CellPoints{Physical} we map the points to the parametric space using the inverse of the geometrical map (implement it for some simple cases or return a meaningful error if not implemented)

function map_points_to_parametric_space(CellPoints{Physical})

# apply inverse map

end

Thus, we could evaluate things like this

xe = get_cell_coordinates(trian) # in the physical space
evaluate(uh,xe) # cheap if ReferenceFE in the physical space, and can be involved and expensive in the parametric

xe = get_cell_ref_coordinates(trian) # in the parametric space, not implemented yet
evaluate(uh,xe) # cheap if ReferenceFE space defined in the parametric space, and straightforward in the physical

We also want to implement get_cell_ref_coordinates (3 lines of code, but not implemented yet).

On the other hand, we shoud provide methods that given a set of points in the physical space, it generates a CellPoints{Physical} with the nodes located in the cell(s) containing it. That is simple for structured Cartesian and octree meshes, more involved in the most general case.

[jw3126]

Thanks a lot for the comprehensive answer @santiagobadia ! I think being able to evaluate the solution at points (even with bad performance) makes using Gridap extremely lightweight. E.g. one can easily visualize arbitrary slices of the solution without ever leaving julia and without additional API:

using Plots
contour([u([x, y, x+y]) for x in xs, y in ys])

[jw3126]

What you are asking is possible, but not standard in finite element methods. Usually, in finite element methods, you evaluate the solution in the mesh nodes, and you would do it like this:

get_cell_values(u) # nodal values for scalar Lagrangian at each cell

How would one implement this function in case of linear Lagrangian?

[santiagobadia]

It is already implemented. If you use it in your example, you should get the values on the nodes of each cell.

[jw3126]

Ah right, it is available as Gridap.FESpaces.get_cell_values.

[urbainvaes]

What you are asking is possible, but not standard in finite element methods.

I think I would disagree with this. In my opinion, one advantage of the finite element method is precisely that its output is a true function defined everywhere on the domain, unlike e.g. finite difference methods. And for visualization purposes, it would be very nice to be able to evaluate the solution on a finer grid than the element nodes.

[fverdugo]

Hi @urbainvaes ! The FE functions in Gridap can be evaluated at any point you want. However, only a set of practical cases are exposed via the high-level API, but with lower level methods you can implement almost anything you can imagine.

For instance, if you have a Triangulation, trian, and for each cell in this triangulation you have a point or a vector of points cell_to_point, you can create

x = CellPoint(cell_to_point,trian,ReferenceDomain())
# or 
x = CellPoint(cell_to_point,trian,PhysicalDomain())

depending where your points are defined.

then you can evaluate any CellField on these points, e.g., uh(x).

[urbainvaes]

Hi @fverdugo !

Thank you very much for such a quick response! I've only just started using Gridap, and it feels absolutely great (and I say this after years of using FreeFem++ and Fenics)!

I'm not completely sure that I understand your answer. Say I want to evaluate "uh([1, 1])". From what you are saying, it seems that I should first figure out the cell in which my physical point [1, 1] is?

[fverdugo]

it seems that I should first figure out the cell in which my physical point [1, 1] is?

Yes, this is the part that is not implemented now in gridap. There is some WIP in PR #523

[santiagobadia]

Hi @urbainvaes ,

Thanks for the interest in Gridap. In the comment you point out, I was saying that the standard way to evaluate FE functions was cellwise, since FE functions are cellwise. Sure, a FE function can be evaluated anywhere in the domain but you need to know in which cell the point is located to compute its value. In Gridap we can do quite a lot of things without needing the functionality you mention.

In any case, it would be interesting to have the possibility to given a point provide the value of the FE function. It involves a search algorithm as first step. As @fverdugo says, there is a PR in progress about this topic.

It is a FE thing, not particular of Gridap. If you are more familiar with spectral elements, the situation there is different because functions have a global definition.

[urbainvaes]

Hi @fverdugo and @santiagobadia,

Thank you very much for your responses! It's great to see that there is indeed already a PR on this subjet. I agree that the first step and main difficulty is to implement a search algorithm to find the cell in which the point lies. In the past, I have programmed this using a line search (drawing a line between a given element in the mesh and the evaluation point, and moving from element to element along this line until we reach the element that contains the evaluation point). But the associated complexity is O(n) if the domain contains roughly n elements per direction, which is not ideal. I see that you are doing this search with kd trees, which should yield better, logarithmic complexity, I imagine...