Skip to content
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

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

Open
ovanvincq opened this issue Feb 13, 2024 · 0 comments
Open

Comments

@ovanvincq
Copy link

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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant