Skip to content

Commit 18be278

Browse files
frederichechtfrederichecht
committed
correct pb of eps in trunc function, and change to close seach in close and nearest case (more safe).
1 parent ae1bef1 commit 18be278

File tree

2 files changed

+46
-41
lines changed

2 files changed

+46
-41
lines changed

‎src/fflib/lgmesh.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -550,7 +550,7 @@ AnyType Op_trunc_mesh::Op::operator()(Stack stack) const {
550550

551551
*mp=mps;
552552
if (verbosity>1)
553-
cout << " -- Trunc mesh: Nb of Triangle = " << kk << " label=" <<label <<endl;
553+
cout << " -- Trunc mesh(2d): Nb of Triangle = " << kk << " label=" <<label <<endl;
554554
Mesh * pmsh = new Mesh(Th,split,false,label);
555555
if(renum)
556556
pmsh->renum();

‎src/fflib/msh3.cpp‎

Lines changed: 45 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,7 @@ void TestSameVertexMesh3(const Mesh3 &Th3, const double &hseuil, const R3 &Psup,
9797
for (int ii = 0; ii < Th3.nv; ii++) {
9898
const R3 r3vi(Th3.vertices[ii].x, Th3.vertices[ii].y, Th3.vertices[ii].z);
9999
const Vertex3 &vi(r3vi);
100-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
100+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
101101

102102
if (!pvi) {
103103
v[nv_t].x = vi.x;
@@ -130,7 +130,7 @@ void TestSameTetrahedraMesh3(const Mesh3 &Th3, const double &hseuil, const R3 &P
130130
for (int jj = 0; jj < 4; jj++) iv[jj] = Th3.operator( )(K[jj]);
131131

132132
const R3 vi(K(PtHat));
133-
Vertex3 *pvi = gtree_t->ToClose(vi, hseuil);
133+
Vertex3 *pvi = gtree_t->ToClose(vi, hseuil,true);
134134

135135
if (!pvi) {
136136
vt[nt_t].x = vi.x;
@@ -165,7 +165,7 @@ void TestSameTetrahedraMesh3(const Mesh3 &Th3, const double &hseuil, const R3 &P
165165
for (int jj = 0; jj < 4; jj++) iv[jj] = Th3.operator( )(K[jj]);
166166

167167
const R3 vi(K(PtHat));
168-
Vertex3 *pvi = gtree_t->ToClose(vi, hseuil);
168+
Vertex3 *pvi = gtree_t->ToClose(vi, hseuil,true);
169169

170170
if (!pvi) {
171171
vt[nt_t].x = vi.x;
@@ -198,7 +198,7 @@ void TestSameTriangleMesh3(const Mesh3 &Th3, const double &hseuil, const R3 &Psu
198198
for (int jj = 0; jj < 3; jj++) iv[jj] = Th3.operator( )(K[jj]);
199199

200200
const R3 vi(K(PtHat));
201-
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil);
201+
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil,true);
202202

203203
if (!pvi) {
204204
vbe[nbe_t].x = vi.x;
@@ -233,7 +233,7 @@ void TestSameTriangleMesh3(const Mesh3 &Th3, const double &hseuil, const R3 &Psu
233233
for (int jj = 0; jj < 3; jj++) iv[jj] = Th3.operator( )(K[jj]);
234234

235235
const R3 vi(K(PtHat));
236-
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil);
236+
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil,true);
237237

238238
if (!pvi) {
239239
vbe[nbe_t].x = vi.x;
@@ -614,7 +614,7 @@ Mesh3 *TestElementMesh3_patch(const Mesh3 &Th3) {
614614

615615
for (int ii = 0; ii < Th3.nv; ii++) {
616616
const Vertex3 &vi(Th3.vertices[ii]);
617-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
617+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
618618

619619
if (!pvi) {
620620
v[nbv].x = vi.x;
@@ -1883,7 +1883,7 @@ Mesh3 *GluMesh3(listMesh3 const &lst) {
18831883

18841884
for (int ii = 0; ii < Th3.nv; ii++) {
18851885
const Vertex3 &vi(Th3.vertices[ii]);
1886-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
1886+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
18871887

18881888
if (!pvi) {
18891889
v[nbv].x = vi.x;
@@ -1898,7 +1898,7 @@ Mesh3 *GluMesh3(listMesh3 const &lst) {
18981898
for (int k = 0; k < Th3.nt; k++) {
18991899
const Tet &K(Th3.elements[k]);
19001900
int iv[4];
1901-
for (int i = 0; i < 4; i++) iv[i] = gtree->ToClose(K[i], hseuil) - v;
1901+
for (int i = 0; i < 4; i++) iv[i] = gtree->ToClose(K[i], hseuil,true) - v;
19021902
(tt++)->set(v, iv, K.lab);
19031903
}
19041904
}
@@ -1926,7 +1926,7 @@ Mesh3 *GluMesh3(listMesh3 const &lst) {
19261926
for (int i = 0; i < 3; i++) iv[i] = Th3.operator( )(K[i]);
19271927

19281928
const Vertex3 vi(K(PtHat));
1929-
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil_border);
1929+
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil_border,true);
19301930
if (!pvi) {
19311931
becog[nbe].x = vi.x;
19321932
becog[nbe].y = vi.y;
@@ -1936,7 +1936,7 @@ Mesh3 *GluMesh3(listMesh3 const &lst) {
19361936

19371937
int igluv[3];
19381938
for (int i = 0; i < 3; i++)
1939-
igluv[i] = gtree->ToClose(K[i], hseuil) - v; // NumSom[iv[0]+nbv0];
1939+
igluv[i] = gtree->ToClose(K[i], hseuil,true) - v; // NumSom[iv[0]+nbv0];
19401940

19411941
(bb++)->set(v, igluv, K.lab);
19421942
}
@@ -2124,7 +2124,7 @@ MeshS *GluMesh(listMeshT<MeshS> const &lst) {
21242124

21252125
for (int ii = 0; ii < Th.nv; ii++) {
21262126
const V &vi(Th(ii));
2127-
V *pvi = gtree->ToClose(vi, hseuil);
2127+
V *pvi = gtree->ToClose(vi, hseuil,true);
21282128
if (!pvi) {
21292129
v[nbv].x = vi.x;
21302130
v[nbv].y = vi.y;
@@ -2156,7 +2156,7 @@ MeshS *GluMesh(listMeshT<MeshS> const &lst) {
21562156
const T &K(Th[k]);
21572157
const R3 r3vi(K(PtHat1));
21582158
const V &vi(r3vi);
2159-
V *pvi = gtree_e->ToClose(vi, hseuil_border);
2159+
V *pvi = gtree_e->ToClose(vi, hseuil_border,true);
21602160
if (!pvi) {
21612161
becog1[nbt].x = vi.x;
21622162
becog1[nbt].y = vi.y;
@@ -2166,7 +2166,7 @@ MeshS *GluMesh(listMeshT<MeshS> const &lst) {
21662166

21672167
int igluv[T::nv];
21682168
for (int i = 0; i < (T::nv); i++)
2169-
igluv[i] = gtree->ToClose(K[i], hseuil) - v;
2169+
igluv[i] = gtree->ToClose(K[i], hseuil,true) - v;
21702170
(tt++)->set(v, igluv, K.lab);
21712171
}
21722172
}
@@ -2179,7 +2179,7 @@ MeshS *GluMesh(listMeshT<MeshS> const &lst) {
21792179

21802180
const R3 r3vi(K(PtHat2));
21812181
const V &vi(r3vi);
2182-
V *pvi = gtree_be->ToClose(vi, hseuil_border);
2182+
V *pvi = gtree_be->ToClose(vi, hseuil_border,true);
21832183
if (!pvi) {
21842184
becog2[nbe].x = vi.x;
21852185
becog2[nbe].y = vi.y;
@@ -2189,7 +2189,7 @@ MeshS *GluMesh(listMeshT<MeshS> const &lst) {
21892189

21902190
int igluv[B::nv];
21912191
for (int i = 0; i < (B::nv); i++)
2192-
igluv[i] = gtree->ToClose(K[i], hseuil) - v;
2192+
igluv[i] = gtree->ToClose(K[i], hseuil,true) - v;
21932193

21942194
(bb++)->set(v, igluv, K.lab);
21952195
}
@@ -2276,7 +2276,7 @@ MeshL *GluMesh(listMeshT<MeshL> const &lst) {
22762276

22772277
for (int ii = 0; ii < Th.nv; ii++) {
22782278
const V &vi(Th(ii));
2279-
V *pvi = gtree->ToClose(vi, hseuil);
2279+
V *pvi = gtree->ToClose(vi, hseuil,true);
22802280
if (!pvi) {
22812281
v[nbv].x = vi.x;
22822282
v[nbv].y = vi.y;
@@ -2312,7 +2312,7 @@ MeshL *GluMesh(listMeshT<MeshL> const &lst) {
23122312

23132313
const R3 r3vi(K(PtHat1));
23142314
const V &vi(r3vi);
2315-
V *pvi = gtree_e->ToClose(vi, hseuil_border);
2315+
V *pvi = gtree_e->ToClose(vi, hseuil_border,true);
23162316
if (!pvi) {
23172317
becog1[nbt].x = vi.x;
23182318
becog1[nbt].y = vi.y;
@@ -2322,7 +2322,7 @@ MeshL *GluMesh(listMeshT<MeshL> const &lst) {
23222322

23232323
int igluv[T::nv];
23242324
for (int i = 0; i < 2/*T::nv*/; i++)
2325-
igluv[i] = gtree->ToClose(K[i], hseuil) - v;
2325+
igluv[i] = gtree->ToClose(K[i], hseuil,true) - v;
23262326
(tt++)->set(v, igluv, K.lab);
23272327
}
23282328
}
@@ -2334,15 +2334,15 @@ MeshL *GluMesh(listMeshT<MeshL> const &lst) {
23342334
const B &K(Th.be(k));
23352335
const V &vi(Th(Th.operator( )(K[0]) ));
23362336

2337-
V *pvi = gtree_be->ToClose(vi, hseuil_border);
2337+
V *pvi = gtree_be->ToClose(vi, hseuil_border,true);
23382338
if (!pvi) {
23392339
becog2[nbe].x = vi.x;
23402340
becog2[nbe].y = vi.y;
23412341
becog2[nbe].z = vi.z;
23422342
becog2[nbe].lab = vi.lab;
23432343
gtree_be->Add(becog2[nbe++]);
23442344
int igluv[B::nv];
2345-
igluv[0] = gtree->ToClose(K[0], hseuil) - v;
2345+
igluv[0] = gtree->ToClose(K[0], hseuil,true) - v;
23462346
(bb++)->set(v, igluv, K.lab);
23472347
}
23482348
}
@@ -4041,7 +4041,7 @@ void OrderVertexTransfo_hcode_nv_gtree(const int &tab_nv, const R3 &bmin, const
40414041
for (int ii = 0; ii < tab_nv; ii++) {
40424042
const R3 r3vi(tab_XX[ii], tab_YY[ii], tab_ZZ[ii]);
40434043
const Vertex3 &vi(r3vi);
4044-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
4044+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
40454045
if (!pvi) {
40464046
v[nv_t].x = vi.x;
40474047
v[nv_t].y = vi.y;
@@ -4116,7 +4116,7 @@ void PointCommun_hcode_gtree(const int &dim, const int &NbPoints, const int &poi
41164116
for (int ii = 0; ii < NbPoints; ii++) {
41174117
const R3 r3vi(Coord_Point[ii][0], Coord_Point[ii][1], Coord_Point[ii][2]);
41184118
const Vertex3 &vi(r3vi);
4119-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
4119+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
41204120

41214121
if (!pvi) {
41224122
v[np].x = vi.x;
@@ -4145,7 +4145,7 @@ void PointCommun_hcode_gtree(const int &dim, const int &NbPoints, const int &poi
41454145
const R3 r3vi(Coord_Point[ii][0], Coord_Point[ii][1], Coord_Point[ii][2]);
41464146
// int label = label_point[ii];
41474147
const Vertex3 &vi(r3vi);
4148-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
4148+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
41494149
ind_multiple[pvi - v] = ind_multiple[pvi - v] + 1;
41504150
}
41514151

@@ -5348,7 +5348,7 @@ MeshS *truncmesh(const MeshS &Th, const long &kksplit, int *split, bool WithMort
53485348
<< newedge << " edge on border "<< nbeb<< endl;
53495349

53505350

5351-
double hseuil = (hmin / kksplit) / 1000.;
5351+
double hseuil = (hmin / kksplit) / 10.;
53525352
nbe = (nbeb+newedge)*kksplit;
53535353
if (verbosity > 5)
53545354
cout << " - Before trunc, mesh has " << Th.nbe << " new boundary edges, " << nbeb
@@ -5427,7 +5427,7 @@ MeshS *truncmesh(const MeshS &Th, const long &kksplit, int *split, bool WithMort
54275427
// first build old point to keep the numbering order for DDM ...
54285428
for (int i = 0, k = 0; i < Th.nv; i++) {
54295429
if (takevertex[i] >= 0) {
5430-
Vertex3 *pvi = gtree->ToClose(Th(i), hseuil);
5430+
Vertex3 *pvi = gtree->ToClose(Th(i), hseuil,true);
54315431
if (!pvi) {
54325432
(R3 &)vertices[np] = Th(i);
54335433
vertices[np].lab = Th(i).lab;
@@ -5451,7 +5451,7 @@ MeshS *truncmesh(const MeshS &Th, const long &kksplit, int *split, bool WithMort
54515451

54525452
// new points
54535453
for (int iv = 0; iv < nvsub; iv++) {
5454-
Vertex3 *pvi = gtree->ToClose(vertextrisub[iv], hseuil);
5454+
Vertex3 *pvi = gtree->ToClose(vertextrisub[iv], hseuil,true);
54555455

54565456
if (!pvi) {
54575457
vertices[np].x = vertextrisub[iv].x;
@@ -5561,7 +5561,12 @@ MeshS *truncmesh(const MeshS &Th, const long &kksplit, int *split, bool WithMort
55615561

55625562
for (int iv = 0; iv < nv1Dsub; iv++) {
55635563
const Vertex3 &vi(vertexedgesub[iv]);
5564-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
5564+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
5565+
if( ! pvi ) {
5566+
cout << " bug !!!" << endl;
5567+
cout << iv << " " << vi << " " << vertex1Dsub[iv].x << " "<< ivv[0] << " " << ivv[1] << endl;
5568+
ffassert(pvi);
5569+
}
55655570
assert(pvi);
55665571
newindex[iv] = pvi - vertices;
55675572
}
@@ -5651,7 +5656,7 @@ AnyType Op_trunc_meshS::Op::operator( )(Stack stack) const {
56515656

56525657
// *mp=mps;
56535658
if (verbosity > 1) {
5654-
cout << " -- Trunc mesh: Nb of Surface Triangles = " << kk << " label=" << label << endl;
5659+
cout << " -- Trunc meshS: Nb of Surface Triangles = " << kk << " label=" << label << endl;
56555660
}
56565661

56575662
if (pn2o) {
@@ -5803,7 +5808,7 @@ MeshL *truncmesh(const MeshL &Th, const long &kksplit, int *split, bool WithMort
58035808
}
58045809
}
58055810

5806-
double hseuil = (hmin / kksplit) / 1000.;
5811+
double hseuil = (hmin / kksplit) / 10.;
58075812
if (verbosity > 5)
58085813
cout << "Before trunc hseuil=" << hseuil << endl;
58095814

@@ -5868,7 +5873,7 @@ MeshL *truncmesh(const MeshL &Th, const long &kksplit, int *split, bool WithMort
58685873
// first build old point to keep the numbering order for DDM ...
58695874
for (int i = 0, k = 0; i < Th.nv; i++) {
58705875
if (takevertex[i] >= 0) {
5871-
V *pvi = gtree->ToClose(Th(i), hseuil);
5876+
V *pvi = gtree->ToClose(Th(i), hseuil,true);
58725877
if (!pvi) {
58735878
(R3 &)vertices[np] = Th(i);
58745879
vertices[np].lab = Th(i).lab;
@@ -5894,7 +5899,7 @@ MeshL *truncmesh(const MeshL &Th, const long &kksplit, int *split, bool WithMort
58945899

58955900
// new points
58965901
for (int iv = 0; iv < nvsub; iv++) {
5897-
V *pvi = gtree->ToClose(vertexedgesub[iv], hseuil);
5902+
V *pvi = gtree->ToClose(vertexedgesub[iv], hseuil,true);
58985903

58995904
if (!pvi) {
59005905
vertices[np].x = vertexedgesub[iv].x;
@@ -5996,7 +6001,7 @@ AnyType Op_trunc_meshL::Op::operator( )(Stack stack) const {
59966001

59976002
// *mp=mps;
59986003
if (verbosity > 1) {
5999-
cout << " -- Trunc mesh: Nb of line Edges = " << kk << " label=" << label << endl;
6004+
cout << " -- Trunc meshL: Nb of line Edges = " << kk << " label=" << label << endl;
60006005
}
60016006

60026007
if (pn2o) {
@@ -6177,7 +6182,7 @@ Mesh3 *truncmesh(const Mesh3 &Th, const long &kksplit, int *split, bool kk, cons
61776182

61786183
ffassert(nbfi % 2 == 0);
61796184
nface = nbeee + nbfi / 2;
6180-
double hseuil = (hmin / kksplit) / 1000.;
6185+
double hseuil = (hmin / kksplit) / 10.;
61816186
if (verbosity > 5) {
61826187
cout << " number of not intern boundary faces = " << nbeee << ", all faces = " << nbe
61836188
<< ", hseuil=" << hseuil << endl;
@@ -6293,7 +6298,7 @@ Mesh3 *truncmesh(const Mesh3 &Th, const long &kksplit, int *split, bool kk, cons
62936298
if (takevertex[i] >= 0) {
62946299
double heps = hseuil;
62956300
for (int j = 0; j < 3; ++j) {
6296-
Vertex3 *pvi = gtree->ToClose(Th(i), heps);
6301+
Vertex3 *pvi = gtree->ToClose(Th(i), heps,true);
62976302
if (!pvi) {
62986303
(R3 &)v[np] = Th(i);
62996304
v[np].lab = Th(i).lab;
@@ -6650,7 +6655,7 @@ AnyType Op_trunc_mesh3::Op::operator( )(Stack stack) const {
66506655

66516656
// *mp=mps;
66526657
if (verbosity > 1) {
6653-
cout << " -- Trunc mesh: Nb of Tetrahedrons = " << kk << " label=" << label << endl;
6658+
cout << " -- Trunc mesh3: Nb of Tetrahedrons = " << kk << " label=" << label << endl;
66546659
}
66556660

66566661
Mesh3 *Tht = truncmesh(Th, kkksplit, split, false, label, precis_mesh, orientation, cleanmesh, removeduplicate);
@@ -7289,7 +7294,7 @@ Mesh3 *GluMesh3tab(KN< pmesh3 > *const &tab, long const &lab_delete, bool const
72897294

72907295
for (int ii = 0; ii < Th3.nv; ii++) {
72917296
const Vertex3 &vi(Th3.vertices[ii]);
7292-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
7297+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
72937298

72947299
if (!pvi) {
72957300
v[nbv].x = vi.x;
@@ -7304,7 +7309,7 @@ Mesh3 *GluMesh3tab(KN< pmesh3 > *const &tab, long const &lab_delete, bool const
73047309
for (int k = 0; k < Th3.nt; k++) {
73057310
const Tet &K(Th3.elements[k]);
73067311
int iv[Tet::nea];
7307-
for (int iea = 0; iea < Tet::nea; iea++) iv[iea] = gtree->ToClose(K[iea], hseuil) - v;
7312+
for (int iea = 0; iea < Tet::nea; iea++) iv[iea] = gtree->ToClose(K[iea], hseuil,true) - v;
73087313

73097314
(t++)->set(v, iv, K.lab);
73107315
}
@@ -7329,7 +7334,7 @@ Mesh3 *GluMesh3tab(KN< pmesh3 > *const &tab, long const &lab_delete, bool const
73297334
for (int ii = 0; ii < 3; ii++) iv[ii] = Th3.operator( )(K[ii]);
73307335

73317336
const Vertex3 &vi(K(PtHat));
7332-
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil_border);
7337+
Vertex3 *pvi = gtree_be->ToClose(vi, hseuil_border,true);
73337338
if (!pvi) {
73347339
becog[nbe].x = vi.x;
73357340
becog[nbe].y = vi.y;
@@ -7338,7 +7343,7 @@ Mesh3 *GluMesh3tab(KN< pmesh3 > *const &tab, long const &lab_delete, bool const
73387343
gtree_be->Add(becog[nbe++]);
73397344

73407345
int igluv[3];
7341-
for (int i = 0; i < 3; i++) igluv[i] = gtree->ToClose(K[i], hseuil) - v;
7346+
for (int i = 0; i < 3; i++) igluv[i] = gtree->ToClose(K[i], hseuil,true) - v;
73427347

73437348
(bb++)->set(v, igluv, K.lab);
73447349
}
@@ -7377,7 +7382,7 @@ Mesh3 *GluMesh3tab(KN< pmesh3 > *const &tab, long const &lab_delete, bool const
73777382
for (int p = 0; p < 3; p++) {
73787383
const int iv = Th3.operator()(K[p]);
73797384
const Vertex3 &vi(Th3.vertices[iv]);
7380-
Vertex3 *pvi = gtree->ToClose(vi, hseuil);
7385+
Vertex3 *pvi = gtree->ToClose(vi, hseuil,true);
73817386
if (!pvi) {
73827387
v[nbv].x = vi.x;
73837388
v[nbv].y = vi.y;

0 commit comments

Comments
 (0)