-
Notifications
You must be signed in to change notification settings - Fork 99
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Evaluate PDE solution #345
Comments
I would also really like to know the answer to this problem. |
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. |
Thank you for this fantastic answer. It might be worth mentioning this in the tutorial/documentation - unless I missed it of course. |
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 This way:
Thus, we could evaluate things like this
We also want to implement 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. |
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]) |
How would one implement this function in case of linear Lagrangian? |
It is already implemented. If you use it in your example, you should get the values on the nodes of each cell. |
Ah right, it is available as |
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. |
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 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 |
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 " |
Yes, this is the part that is not implemented now in gridap. There is some WIP in PR #523 |
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. |
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 |
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?
I tried the following, both of which did not work:
The text was updated successfully, but these errors were encountered: