Miklós, thanks for your answer. Le 09/06/16 à 14:28, Homolya, Miklós a écrit :
It would be much easier to help if you could provide a minimal failing example that exposes the bug you are experiencing. Preferably as a GitHub issue.
It would be much easier to debug as well if I could, but when I try to isolate the bug it disappears... What I'm doing: considering a time-dependent equation, my code solves the equation for one time step using Firedrake, calls a third-party software to generate a metric based on the solution, then generates a mesh adapted to this metric using petsc, then interpolates the solution onto the new mesh using ugly loops and point evaluation, then solves the equation for the next time-step and so on a certain number of times. Sometimes during the interpolation step one vertex of the new mesh is found "not in domain" while it obviously is (one vertex every few iterations). To isolate the bug, I have to write the concerned meshes to a file and read them in the minimal code. But then the bug disappears. This could be due to a memory issue (hence valgrind), or rounding errors when writing to the file. I'm going to retry using double precision when writing the files though, to see if it helps. (I checked, mesh or dmplex objects can unfortunately not be serialized.)
So what you're saying is that the point is right on the boundary between two cells? I guess there are two ways the point location can fail. Either the point is not found in any bounding box. Or, the point is found by libspatialindex in a bounding box. But that bucket does not actually contain the cell in question, so then the linear search for point location may fail due to floating point rounding. Which of these two are occurring?
Yes! This is an important distinction.
Let's find out then! But how do I do that ? (cf other email)
If the point is in the bounding box, but found to be outside the cell during the physical coordinates -> reference coordinates transformation, then you could inspect the Newton iteration and the tolerances. Also, which cell types do you have this problem with?
simple triangles.
Please confirm you are not trying to do point evaluations on manifolds.
My domain is a simple unit square, so a manifold I guess, but not the bad ones. -- Nicolas -- Nicolas Barral Dept. of Earth Science and Engineering Imperial College London Royal School of Mines - Office 4.88 London SW7 2AZ