gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
695 stars 95 forks source link

Error with the function `_point_to_cell ` when using a triangular mesh #981

Open ovanvincq opened 7 months ago

ovanvincq commented 7 months ago

Hi,

I encountered a problem with the function evaluate when using a triangular mesh generated by GMSH (the msh file is located at the end of the issue). The geometry is composed of two circles with radius 1 and 2. When I evaluate a field at the point (0.4,0.5), there is an assertion error: Point (0.4,0.5) is not inside any active cell. Of course, this point is indeed in a cell.

using Gridap
using GridapGmsh
model = GmshDiscreteModel("./model.msh");
Ω = Triangulation(model)
test = CellField(x->x[1],Ω)
test(Point(0.4,0.5)) #error

The problem seems to be located in the function _point_to_cell:

using StaticArrays
using NearestNeighbors
cache1,cache2=Gridap.Fields.return_cache(test,Point(0.4,0.5));
searchmethod, kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, table_cache = cache1
x=Point(0.4,0.5)
a=zip(knn(kdtree, SVector(Tuple(x)), searchmethod.num_nearest_vertices, true)...)
id=a.is[1][1]
cells = getindex!(table_cache,vertex_to_cells,id)

By drawing the cells, I get this: bug1 so the point (0.4,0.5) is not in a cell where one of the vertices is the closest point.

The problem is the value of searchmethod.num_nearest_vertices. This variable is set to 1 but must be set to 2:

a=zip(knn(kdtree, SVector(Tuple(x)), 2, true)...)
id=a.is[1][2]
cells2 = getindex!(table_cache,vertex_to_cells,id)

By drawing the cells2, I get this: bug2

Could you fix this issue by setting searchmethod.num_nearest_vertices to 2 for a triangular 2D mesh?

Below is the model.msh file:

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
3 2 2 0
1 0 0 0 1 8 
2 1 0 0 0 
3 2 0 0 0 
1 -1.0000001 -1.0000001 -1e-07 1.0000001 1.0000001 1e-07 1 6 2 2 -2 
2 -2.0000001 -2.0000001 -1e-07 2.0000001 2.0000001 1e-07 1 7 2 3 -3 
1 -2.0000001 -2.0000001 -1e-07 2.0000001 2.0000001 1e-07 1 4 2 2 -1 
2 -1.0000001 -1.0000001 -1e-07 1.0000001 1.0000001 1e-07 1 5 1 1 
$EndEntities
$Nodes
7 15 1 15
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
2 0 0
1 1 0 3
4
5
6
-2.603304941409257e-15 1 0
-1 -1.20980699416795e-15 0
1.592665886326894e-15 -1 0
1 2 0 5
7
8
9
10
11
0.9999999999999976 1.732050807568879 0
-1.000000000000004 1.732050807568875 0
-2 -2.4196139883359e-15 0
-0.9999999999999962 -1.732050807568879 0
1.000000000000002 -1.732050807568876 0
2 1 0 4
12
13
14
15
0.9999999999999987 0.6830127018922197 0
-0.9999999999999986 -0.6830127018922206 0
1.000000000000001 -0.683012701892219 0
-1.000000000000002 0.6830127018922177 0
2 2 0 0
$EndNodes
$Elements
5 33 1 33
0 1 15 1
1 1 
1 1 1 4
2 2 4 
3 4 5 
4 5 6 
5 6 2 
1 2 1 6
6 3 7 
7 7 8 
8 8 9 
9 9 10 
10 10 11 
11 11 3 
2 1 2 18
12 3 12 2 
13 9 13 5 
14 5 15 9 
15 2 14 3 
16 2 12 4 
17 5 13 6 
18 4 15 5 
19 6 14 2 
20 4 12 7 
21 11 14 6 
22 8 15 4 
23 6 13 10 
24 7 12 3 
25 3 14 11 
26 9 15 8 
27 10 13 9 
28 4 7 8 
29 6 10 11 
2 2 2 4
30 1 2 4 
31 6 2 1 
32 1 4 5 
33 1 5 6 
$EndElements