1 // ORIG-DATE: Dec 2007
2 // -*- Mode : c++ -*-
3 //
4 // SUMMARY : Model mesh 3dL
5 // USAGE : LGPL
6 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE
7 // AUTHOR : Frederic Hecht
8 // E-MAIL : frederic.hecht@ann.jussieu.fr
9 //
10
11 /*
12
13 This file is part of Freefem++
14
15 Freefem++ is free software; you can redistribute it and/or modify
16 it under the terms of the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 2.1 of the License, or
18 (at your option) any later version.
19
20 Freefem++ is distributed in the hope that it will be useful,
21 but WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 GNU Lesser General Public License for more details.
24
25 You should have received a copy of the GNU Lesser General Public License
26 along with Freefem++; if not, write to the Free Software
27 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
28
29 Thank to the ARN () FF2A3 grant
30 ref:ANR-07-CIS7-002-01
31 */
32
33 #include <fstream>
34 #include <iostream>
35 #include <cstring>
36 #include "libmesh5.h"
37 #include "ufunction.hpp"
38 #include "error.hpp"
39 #include "RNM.hpp"
40 namespace Fem2D
41 {
42 }
43 #include "Mesh2dn.hpp"
44 #include "Mesh3dn.hpp"
45 #include "MeshSn.hpp"
46 #include "MeshLn.hpp"
47 #include "rgraph.hpp"
48 #include "fem.hpp"
49 #include "PlotStream.hpp"
50
51
52
53 namespace Fem2D
54 {
55 static const int nvfaceSeg[1][3] = {{-1,-1,1}};
56 static const int nvedgeSeg[1][2] = { {0,1} };
57 static const int nvadjSeg[2][1] = { {0},{1} };
58
59 // geometry element for segment ( boundary elements in surface mesh, Rd=3 RdHat=1 )
60 template<> const int (* const GenericElement<DataSeg3>::nvface)[3] = 0 ;
61 template<> const int (* const GenericElement<DataSeg3>::nvedge)[2] = nvedgeSeg; //nvedgeTria ;
62 template<> const int (* const GenericElement<DataSeg3>::nvadj)[1] = nvadjSeg ;
63
64
65
66 template<> const int (* const GenericElement<DataPoint3>::nvface)[3] = 0 ;
67 template<> const int (* const GenericElement<DataPoint3>::nvedge)[2] = 0 ;
68 template<> const int (* const GenericElement<DataPoint3>::nvadj)[1] = 0 ;
69
70 template<> int GenericMesh<EdgeL,BoundaryPointL,Vertex3>::kfind=0;
71 template<> int GenericMesh<EdgeL,BoundaryPointL,Vertex3>::kthrough=0;
72
73 static const int onWhatIsVertex[2][3] = { {1,0,0}, // sommet 0
74 {0,1,0}}; // sommet 1
75
76 template<>
77 const int (* const GenericElement<DataSeg3>::onWhatBorder)[3] = onWhatIsVertex ;
78 template<> const int GenericElement<DataSeg3>::nitemdim[4] = {2,1,0,0 } ;
79
80 const string GsbeginL="MeshS::GSave v0",GsendL="end";
GSave(FILE * ff,int offset) const81 void MeshL::GSave(FILE * ff,int offset) const
82 {
83 PlotStream f(ff);
84
85 f << GsbeginL ;
86 f << nv << nt << nbe;
87 for (int k=0; k<nv; k++) {
88 const Vertex & P = this->vertices[k];
89 f << P.x <<P.y << P.z << P.lab ;
90 }
91
92 for (int k=0; k<nt; ++k) {
93 int iv[EdgeL::nv];
94 const Element & K(this->elements[k]);
95 for(int i=0;i<EdgeL::nv;++i)
96 iv[i]=this->operator()(K[i])+offset;
97 int lab=K.lab;
98 f << iv[0] << iv[1] << lab;
99 }
100 for (int k=0; k<nbe; k++) {
101 const BorderElement & K(this->borderelements[k]);
102 int iv=this->operator()(K[0])+offset;
103 int lab=K.lab;
104 f << iv << lab;
105 }
106 f << GsendL;
107 }
108
109
110
read(istream & f)111 void MeshL::read(istream &f)
112 { // read the mesh
113 int i;
114 string str;
115 int err=0;
116 while(!f.eof())
117 {
118 f >> str;
119 //cout << str << endl;
120 if( str== "Vertices")
121 {
122 f >> nv;
123 assert(!this->vertices );
124 if(verbosity>2)
125 cout << " -- Nb of Vertex " << nv << endl;
126 this->vertices = new Vertex[nv];
127 for (i=0;i<nv;i++)
128 {
129 f >> this->vertices[i];
130 assert(f.good());
131 }
132 }
133 else if (str=="Edges")
134 {
135 f >> nt;
136 assert(this->vertices && !this->elements);
137 this->elements = new Element [nt];
138 mes=0;
139 assert(this->elements);
140 if(verbosity>2)
141 cout << " -- Nb of Elements " << nt << endl;
142 for (int i=0;i<nt;i++)
143 {
144 this->t(i).Read1(f,this->vertices,this->nv);
145 if(this->t(i).mesure()<=0) err++; // Modif FH nov 2014
146 mes += this->t(i).mesure();
147 }
148 }
149 // np present for the moment
150 else if (str=="Points")
151 {
152 mesb=0;
153 cout << " -- No boundary points in ff .mesh format " << endl;
154 /*int kmv=0,ij;
155 f >> nbe;
156 assert(vertices);
157 this->borderelements = new BorderElement[nbe];
158 if(verbosity>2)
159 cout << " -- Nb of border Triangles " << nbe << endl;
160 for (i=0;i<nbe;i++)
161 {
162 this->be(i).Read1(f,this->vertices,this->nv);
163 //mesb += this->be(i).mesure();
164 for(int j=0;j<BorderElement::nv;++j)
165 if(!vertices[ij=this->be(i,j)].lab)
166 {
167 vertices[ij].lab=1;
168 kmv++;
169 }
170 }*/
171 }
172 else if(str[0]=='#') {
173 int c;
174 while ( (c=f.get()) != '\n' && c != EOF)
175 ;
176 }
177 }
178 assert( (nt >= 0 || nbe>=0) && nv>0) ;
179 if(err!=0) {
180 cerr << " MeshL::read: sorry bad mesh. Number of negative Edges " << err << endl;
181 this->~MeshL();
182 ffassert(0);
183 }
184 }
185
186 // no points border in this format
readmsh(ifstream & f,int offset)187 void MeshL::readmsh(ifstream & f,int offset)
188 {
189 int err=0;
190 f >> nv >> nt ; //>> nbe;
191 if(verbosity>2)
192 cout << " GRead : nv " << nv << " " << nt << endl; //" " << nbe << endl;
193 this->vertices = new Vertex[nv];
194 this->elements = new Element[nt];
195 //this->borderelements = new BorderElement[nbe];
196 for (int k=0; k<nv; k++) {
197 Vertex & P = this->vertices[k];
198 f >> P.x >>P.y >> P.z >> P.lab ;
199 }
200 mes=0.;
201 mesb=0.;
202
203 if(nt == 0) {
204 cerr << " A meshL type must have elements " << endl;
205 ffassert(0);exit(1);
206
207 }
208
209 for (int k=0; k<nt; k++) {
210 int i[2],lab;
211 Element & K(this->elements[k]);
212 f >> i[0] >> i[1] >> lab;
213 Add(i,2,offset);
214 K.set(this->vertices,i,lab);
215 mes += K.mesure();
216 err += K.mesure() <0;
217
218 }
219 cout << " -- No boundary points in .msh format " << endl;
220 /*for (int k=0; k<nbe; k++) {
221 int i[1],lab;
222 BorderElement & K(this->borderelements[k]);
223 f >> i[0] >> lab;
224 Add(i,1,offset);
225 K.set(this->vertices,i,lab);
226 mesb += K.mesure();
227
228 }*/
229 if(err!=0)
230 {
231 cerr << " MeshL::readmsh : sorry bad mesh. Number of negative Edges " << err << endl;
232 this->~MeshL();
233 ffassert(0);
234 }
235
236 }
237
238
MeshL(const string filename)239 MeshL::MeshL(const string filename)
240 :mapSurf2Curv(0),mapCurv2Surf(0) {
241 int ok=load(filename);
242 if(verbosity) {
243 cout << "read meshL ok " << ok ;
244 cout << "curve Mesh, num element Edges:= " << nt << ", num Vertice:= " << nv << " num boundary Points:= " << nbe << endl;
245 }
246
247 if (ok) {
248 ifstream f(filename.c_str());
249 if(!f) {
250 cerr << " -- MeshL::MeshL Erreur opening " << filename<<endl;ffassert(0);exit(1);}
251 if(verbosity>2)
252 cout << " -- MeshL: Read On file \"" <<filename<<"\""<< endl;
253 if(filename.rfind(".msh")==filename.length()-4)
254 readmsh(f,-1);
255 else
256 read(f);
257 }
258
259 BuildBound();
260 BuildAdj();
261 Buildbnormalv();
262 BuildjElementConteningVertex();
263 if (nbe==0) {
264 if(verbosity>5)
265 cout << " building of boundary " << endl;
266 BuildBdElem();
267 delete [] TheAdjacencesLink;
268 delete [] BoundaryElementHeadLink;
269 TheAdjacencesLink=0;
270 BoundaryElementHeadLink=0;
271 BuildBound();
272 BuildAdj();
273 Buildbnormalv();
274 BuildjElementConteningVertex();
275 }
276
277 if(verbosity>2) cout << " -- End of read: MeshL mesure = " << mes << endl;
278
279 if(verbosity) cout << " -- MeshL : "<<filename << ", space dimension "<< 3 << ", num Edges elts " << nt << ", num Vertice "
280 << nv << " num Bondary Points " << nbe << endl;
281
282
283 ffassert(mes>=0);
284
285 }
286
287
288
289
load(const string & filename)290 int MeshL::load(const string & filename)
291 {
292 int bin;
293 int ver,inm,dim;
294 int lf=filename.size()+20;
295 KN<char> fileb(lf),filef(lf);
296 char *data = new char[filename.size()+1];
297 size_t ssize = filename.size()+1;
298 char *ptr;
299 char *pfile=data;
300 strncpy( data, filename.c_str(),ssize);
301 ptr = strstr(data,".mesh");
302 if( !ptr ){
303 strcpy(filef,filename.c_str());
304 strcpy(fileb,filef);
305 strcat(filef,".mesh");
306 strcat(fileb,".meshb");
307 if( (inm=GmfOpenMesh(pfile=fileb, GmfRead,&ver,&dim)) )
308 bin=true;
309 else if( (inm=GmfOpenMesh(pfile=filef, GmfRead,&ver,&dim)) )
310 bin=false;
311 else
312 if(verbosity>5){
313 cerr << " Erreur ouverture file " << (char *) fileb << " " << (char *) filef <<endl;
314 return 1;
315 }
316 }
317 else{
318 if( !(inm=GmfOpenMesh(data, GmfRead,&ver,&dim)) ){
319 if(verbosity>5)
320 cerr << " Erreur ouverture file " << (char *) data << endl;
321 return 1;
322 }
323 }
324 // data file is readed and the meshes are initilized
325 int nv=-1,nEdge=-1,nPts=-1;
326 nv=GmfStatKwd(inm,GmfVertices); // vertice
327 nEdge=GmfStatKwd(inm,GmfEdges); // segment elements
328 nPts=0; // GmfStatKwd(inm,GmfVertices); // points border element
329 this->set(nv,nEdge,nPts);
330
331 if(nEdge == 0) {
332 cerr << " A meshL type must have elements " << endl;
333 ffassert(0);exit(1);}
334
335 if(verbosity>1)
336 cout << " -- MeshL(load): "<< (char *) data << ", MeshVersionFormatted:= " << ver << ", space dimension:= "<< dim
337 << ", num Edge elts:= " << nEdge << ", num vertice:= " << nv << " num Points boundaries:= " << nPts << endl;
338
339 if(dim != 3) {
340 cerr << "Err dim == " << dim << " !=3 " <<endl;
341 return 2; }
342 if( nv<=0 && (nEdge <=0 || nPts <0) ) {
343 cerr << " missing data "<< endl;
344 return 3;
345 }
346 int iv[3],lab;
347 float cr[3];
348 int mxlab=0, mnlab=0;
349 // read vertices
350 GmfGotoKwd(inm,GmfVertices);
351 for(int i=0;i<this->nv;++i) {
352 if(ver<2) {
353 GmfGetLin(inm,GmfVertices,&cr[0],&cr[1],&cr[2],&lab);
354 vertices[i].x=cr[0];
355 vertices[i].y=cr[1];
356 vertices[i].z=cr[2];}
357 else
358 GmfGetLin(inm,GmfVertices,&vertices[i].x,&vertices[i].y,&vertices[i].z,&lab);
359 vertices[i].lab=lab;
360 mxlab= max(mxlab,lab);
361 mnlab= min(mnlab,lab);
362 }
363 if(mnlab==0 && mxlab==0 ) {
364 int kmv=0;
365 mes=0;
366 GmfGotoKwd(inm,GmfEdges);
367 for(int i=0;i<nEdge;++i) {
368 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
369 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv);
370 for(int j=0;j<2;++j)
371 if(!vertices[iv[j]-1].lab) {
372 vertices[iv[j]-1].lab=1;
373 kmv++;
374 }
375 for (int j=0;j<2;++j) iv[j]--;
376 elements[i].set(vertices,iv,lab);
377 mes += elements[i].mesure();
378 }
379 if(kmv&& verbosity>1) cout << " Aucun label Hack (FH) ??? => 1 sur les triangle frontiere "<<endl;
380 }
381 else {
382 mes=0;
383 GmfGotoKwd(inm,GmfEdges);
384 for(int i=0;i<nEdge;++i) {
385 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
386 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv);
387 for (int j=0;j<2;++j) iv[j]--;
388 elements[i].set(this->vertices,iv,lab);
389 mes += elements[i].mesure();
390 }
391 }
392 // unused loop...not border points in .mesh
393 mesb=0;
394 GmfGotoKwd(inm,GmfVertices);
395 for(int i=0;i<nPts;++i) {
396 GmfGetLin(inm,GmfVertices,&iv[0],&lab);
397 assert( iv[0]>0 && iv[0]<=nv);
398 borderelements[i].set(this->vertices,iv,lab);
399 mesb += this->borderelements[i].mesure();
400 }
401
402
403 if(verbosity>1)
404 cout << " -- MeshL(load): "<< (char *) data << ", MeshVersionFormatted:= " << ver << ", space dimension:= "<< dim
405 << ", Edges elts:= " << nt << ", num vertice:= " << nv << ", num Points boundaries:= " << nbe << endl;
406
407 GmfCloseMesh(inm);
408 delete[] data;
409 return 0; // OK
410 }
411
412
MeshL(const string filename,bool cleanmesh,bool removeduplicate,bool rebuildboundary,int orientation,double precis_mesh,double ridgeangledetection)413 MeshL::MeshL(const string filename, bool cleanmesh, bool removeduplicate, bool rebuildboundary, int orientation, double precis_mesh, double ridgeangledetection)
414 :mapSurf2Curv(0),mapCurv2Surf(0) {
415
416
417 int ok=load(filename);
418 if(verbosity) {
419 cout << "read meshL ok " << ok << endl;
420 cout << ", nt " << nt << ", nv " << nv << " nbe: = " << nbe << endl;
421 }
422 if(ok) {
423 ifstream f(filename.c_str());
424 if(!f)
425 cerr << " -- MeshL Erreur opening " << filename<<endl;ffassert(0);exit(1);
426 if(verbosity>2)
427 cout << " -- MeshL: Read On file \"" <<filename<<"\""<< endl;
428 if(filename.rfind(".msh")==filename.length()-4)
429 readmsh(f,-1);
430 else
431 read(f);
432 }
433
434 if (cleanmesh) {
435 if(verbosity>3)
436 cout << "before clean meshL, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
437 clean_mesh(precis_mesh, nv, nt, nbe, vertices, elements, borderelements, removeduplicate, rebuildboundary, orientation);
438 if(verbosity>3)
439 cout << "after clean meshL, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
440 }
441
442 BuildBound();
443 BuildAdj();
444 //Buildbnormalv();
445 BuildjElementConteningVertex();
446
447 if (nbe==0) {
448 if(verbosity>5)
449 cout << " building of boundary " << endl;
450 BuildBdElem();
451 delete [] TheAdjacencesLink;
452 delete [] BoundaryElementHeadLink;
453 TheAdjacencesLink=0;
454 BoundaryElementHeadLink=0;
455 BuildBound();
456 BuildAdj();
457 Buildbnormalv();
458 BuildjElementConteningVertex();
459 }
460 if(verbosity>2)
461 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
462 if(verbosity)
463 cout << " -- MeshL : "<<filename << ", d "<< 3 << ", n Edges " << nt << ", n Vtx "
464 << nv << " n Border points " << nbe << endl;
465 ffassert(mes>=0); // add F. Hecht sep 2009.
466 }
467
468
469
MeshL(FILE * f,int offset)470 MeshL::MeshL(FILE *f,int offset)
471 :mapSurf2Curv(0),mapCurv2Surf(0) {
472 GRead(f,offset);// remove 1
473 assert( (nt >= 0 || nbe>=0) && nv>0) ;
474 BuildBound();
475 if(verbosity>2)
476 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
477
478 BuildAdj();
479 Buildbnormalv();
480 BuildjElementConteningVertex();
481
482 if (nbe==0) {
483 if(verbosity>5)
484 cout << " building of boundary " << endl;
485 BuildBdElem();
486 delete [] TheAdjacencesLink;
487 delete [] BoundaryElementHeadLink;
488 TheAdjacencesLink=0;
489 BoundaryElementHeadLink=0;
490 BuildBound();
491 BuildAdj();
492 Buildbnormalv();
493 BuildjElementConteningVertex();
494 }
495 if(verbosity>2)
496 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
497
498 if(verbosity>1)
499 cout << " -- MeshL (File *), d "<< 3 << ", n Tri " << nt << ", n Vtx " << nv << " n Bord " << nbe << endl;
500 ffassert(mes>=0); // add F. Hecht sep 2009.
501 }
502
503
hmin() const504 double MeshL::hmin() const {
505 R3 Pinf(1e100,1e100,1e100),Psup(-1e100,-1e100,-1e100); // Extremite de la boite englobante
506 double hmin=1e10;
507
508 for (int ii=0;ii< this->nv;ii++) {
509 R3 P( vertices[ii].x, vertices[ii].y, vertices[ii].z);
510 Pinf=Minc(P,Pinf);
511 Psup=Maxc(P,Psup);
512 }
513
514 for (int k=0;k<this->nt;k++) {
515
516 if( this->elements[k].mesure() < Norme2(Psup-Pinf)/1e9 ) {
517 const EdgeL & K(this->elements[k]);
518 int iv[2];
519 for(int jj=0; jj <2; jj++)
520 iv[jj] = this->operator()(K[jj]);
521 if(verbosity>2)
522 cout << "EdgeL: " << k << " lenght "<< this->elements[k].mesure() << endl;
523 if(verbosity>2) cout << " A triangleS with a very small edge was created " << endl;
524 return 1;
525 }
526 hmin=min(hmin,this->elements[k].mesure()); // calcul de .lenEdge pour un Mesh3
527
528 }
529 ffassert(hmin>Norme2(Psup-Pinf)/1e9);
530 return hmin;
531 }
532
533
534 // brute force method
GRead(FILE * ff,int offset)535 void MeshL::GRead(FILE * ff,int offset)
536 {
537 PlotStream f(ff);
538 string s;
539 f >> s;
540 ffassert( s== GsbeginL);
541 f >> nv >> nt >> nbe;
542 if(verbosity>2)
543 cout << " GRead : nv " << nv << " " << nt << " " << nbe << endl;
544 this->vertices = new Vertex[nv];
545 this->elements = new Element [nt];
546 this->borderelements = new BorderElement[nbe];
547 for (int k=0; k<nv; k++) {
548 Vertex & P = this->vertices[k];
549 f >> P.x >>P.y >> P.z >> P.lab ;
550 }
551 mes=0.;
552 mesb=0.;
553
554 if(nt == 0) {
555 cerr << " A meshL type must have elements " << endl;
556 ffassert(0);exit(1);}
557
558
559 for (int k=0; k<nt; k++) {
560 int i[2],lab;
561 Element & K(this->elements[k]);
562 f >> i[0] >> i[1] >> lab;
563 Add(i,2,offset);
564 K.set(this->vertices,i,lab);
565 mes += K.mesure();
566
567 }
568 for (int k=0; k<nbe; k++) {
569 int i[2],lab;
570 BorderElement & K(this->borderelements[k]);
571 f >> i[0] >> lab;
572 Add(i,1,offset);
573 K.set(this->vertices,i,lab);
574 mesb += K.mesure();
575
576 }
577 f >> s;
578 ffassert( s== GsendL);
579 }
580
581
Find(Rd P,R1 & Phat,bool & outside,const Element * tstart) const582 const MeshL::Element * MeshL::Find( Rd P, R1 & Phat,bool & outside,const Element * tstart) const
583
584 {
585 if(searchMethod != 2)
586 {
587 GenericDataFindBoundary<GMesh> * gdfb=Buildgdfb();
588 if(gdfb )
589 {
590 double l[2];
591 int loutside;// 0 inside, 1 out close, 2, out fare, , -1 inside
592 int itt =gdfb->Find(P,l,loutside);
593 outside=loutside;
594 Phat=R1(l[1]);
595 Element &K=(this->elements)[itt];
596 if( verbosity > 9)
597 cout << " - Find "<< P << " -> " << K(Phat) << " " << loutside << " k= " << itt
598 << " dist =" << (P-K(Phat)).norme() << " :: " << Phat << endl;
599 return itt<0 ? 0: this->elements + itt; // outside
600
601 }
602 }
603 // rewritie FH 31 jan 2020 ..
604 static int count =0;
605 if( verbosity && count++< 5 )
606 cerr << " MeshL::Find warning brute force to day " << endl;
607 // find the neast points ..
608 double dmin2 = 1e200;
609 R1 Phm;
610 bool out = true;
611 int nopt=-1;
612 for (int i=0;i<nt;i++) {
613 kthrough++;
614 const EdgeL & K(this->elements[i]);
615 R3 A(K[0]),B(K[1]), AB(A,B);
616 R3 AP(A,P);
617 double lab2 = AB.norme2();
618 double l = min(1.,max(0.,(AB,AP)/lab2));
619 R1 Ph(l);
620 R3 Pt=K(Ph);
621 double d2=R3(P,Pt).norme2();
622 if(dmin2>d2)
623 {
624 dmin2 = d2;
625 Phm=Ph;
626 nopt =i;
627 out = d2 < lab2*1e-5; // BofBof FH ...
628 }
629
630 }
631 ffassert( nopt>=0);
632 Phat=Phm;
633 return this->elements+nopt; // outside
634 }
635
636
MeshL(int nnv,int nnt,int nnbe,Vertex3 * vv,EdgeL * tt,BoundaryPointL * bb,bool cleanmesh,bool removeduplicate,bool rebuildboundary,int orientation,double precis_mesh,double ridgeangledetection)637 MeshL::MeshL(int nnv, int nnt, int nnbe, Vertex3 *vv, EdgeL *tt, BoundaryPointL *bb, bool cleanmesh, bool removeduplicate, bool rebuildboundary, int orientation, double precis_mesh, double ridgeangledetection)
638 :mapSurf2Curv(0),mapCurv2Surf(0)
639 {
640 nv = nnv;
641 nt = nnt;
642 nbe =nnbe;
643 vertices = vv;
644 elements = tt;
645 borderelements = bb;
646 mes=0.;
647 mesb=0.;
648
649 for (int i=0;i<nt;i++)
650 mes += this->elements[i].mesure();
651 for (int i=0;i<nbe;i++)
652 mesb += this->be(i).mesure();
653
654 if (cleanmesh) {
655 if(verbosity>3)
656 cout << "before clean meshL, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
657 clean_mesh(precis_mesh, nv, nt, nbe, vertices, elements, borderelements, removeduplicate, rebuildboundary, orientation);
658 if(verbosity>3)
659 cout << "after clean meshL, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
660 }
661 // BuildCurvBasis();
662 BuildBound();
663 BuildAdj();
664 Buildbnormalv();
665 BuildjElementConteningVertex();
666 if (nbe==0) {
667 if(verbosity>5)
668 cout << " building of boundary " << endl;
669 BuildBdElem();
670 delete [] TheAdjacencesLink;
671 delete [] BoundaryElementHeadLink;
672 TheAdjacencesLink=0;
673 BoundaryElementHeadLink=0;
674 BuildBound();
675 BuildAdj();
676 Buildbnormalv();
677 BuildjElementConteningVertex();
678 }
679
680 if(verbosity>1)
681 cout << " -- End of read meshL: mesure = " << mes << " border mesure " << mesb << endl;
682
683 assert(mes>=0.);
684 assert(mesb==0.);
685 }
686
687
Save(const string & filename) const688 int MeshL::Save(const string & filename) const
689 {
690 int ver = GmfDouble, outm;
691 if ( !(outm = GmfOpenMesh(filename.c_str(),GmfWrite,ver,3)) ) {
692 cerr <<" -- MeshL**::Save UNABLE TO OPEN :"<< filename << endl;
693 return(1);
694 }
695 float fx,fy,fz;
696 // write vertice (meshL)
697 GmfSetKwd(outm,GmfVertices,nv);
698 for (int k=0; k<nv; k++) {
699 const Vertex & P = vertices[k];
700 GmfSetLin(outm,GmfVertices,fx=P.x,fy=P.y,fz=P.z,P.lab);
701 }
702 // write triangles (meshS)
703 GmfSetKwd(outm,GmfEdges,nt);
704 for (int k=0; k<nt; k++) {
705 const EdgeL & K(elements[k]);
706 int i0=this->operator()(K[0])+1;
707 int i1=this->operator()(K[1])+1;
708 int lab=K.lab;
709 GmfSetLin(outm,GmfEdges,i0,i1,lab);
710 }
711
712
713
714
715 // no boundary points ?
716
717 GmfCloseMesh(outm);
718 return (0);
719 }
720
721
722
723 // determine the bounder points list for meshL
BuildBdElem(double angle)724 void MeshL::BuildBdElem(double angle) {
725
726
727 delete [] borderelements; // to remove the previous pointers
728 borderelements = new BoundaryPointL[2 * nt]; // 2 * nt upper evaluated
729
730 HashTable<SortArray<int, 1>, int> pointI(2 * nt, nt);
731 int* AdjLink = new int[2 * nt];
732
733 int nbeL=0,nbiL=0,nk=0;
734 // Build border points from the edge list
735 for (int i = 0; i < nt; i++)
736 for (int j = 0; j < 2; j++) {
737 int jt = j, it = ElementAdj(i, jt);
738 EdgeL &K(elements[i]); // current element
739 // True border point -> no adjacence / on domain border
740 if ((it == i || it < 0)) {
741 int iv[1];
742 iv[0] = this->operator () (K [EdgeL::nvedge[0][j]]);
743 if(verbosity>15)
744 cout << " the edge " << iv[0] << " is a boundary " << endl;
745 be(nbeL++).set(vertices,iv,K.lab);
746
747 }
748 // internal point -- check angular and no manifold
749 else {
750 EdgeL &K_adj(elements[it]); // adjacence element
751 int iv[1];
752 iv[0] = this->operator () (K [EdgeL::nvedge[0][j]]);
753 SortArray<int, 1> key(iv[0]);
754 typename HashTable<SortArray<int,1>,int>::iterator p= pointI.find(key);
755 if (!p) {
756 //edge element
757 R3 A(K[0]),B(K[1]);
758 R3 E(B-A);
759 E/=E.norme();
760 // adj edge element
761 R3 A_adj(K_adj[0]),B_adj(K_adj[1]);
762 R3 E_adj(B_adj-A_adj);
763 E_adj/=E_adj.norme();
764
765 R pdt = (E,E_adj); // scalar product
766 pdt = acos(pdt); // radian angle (Normal,Normal_adj)
767 if(verbosity>15)
768 cout << "Element num: " << i << " N " << E << " Element adjacent num: " << it << " E_adj " << E_adj << " angle between N N_adj = " << pdt <<endl;
769
770 if(pdt >= angle) {
771 if(verbosity>15)
772 cout << " the border point " <<nbeL <<": [" << iv[0] << " " << "] is a boundary with the angular criteria" << endl;
773 int lab = min(K.lab, K_adj.lab);
774 be(nbeL).set(vertices,iv,lab);
775 pointI.add(key, nbeL++);
776 }
777 }
778 }
779 nk++; // increment the total edge jump --- nt * 2
780
781 }
782 assert(nt*2==nk);
783 delete [] AdjLink;
784 // update the number of border points
785 nbe = nbeL;
786 if (verbosity>5)
787 cout << " Building border point from meshS nbe: "<< nbeL << " nbi: " << nbiL << endl;
788
789 BuildBound();
790 delete []TheAdjacencesLink;
791 delete [] BoundaryElementHeadLink;
792 TheAdjacencesLink=0;
793 BoundaryElementHeadLink=0;
794 BuildAdj();
795 Buildbnormalv();
796 BuildjElementConteningVertex();
797
798 }
799
MeshL(const Serialize & serialized)800 MeshL::MeshL(const Serialize &serialized)
801 :GenericMesh<EdgeL,BoundaryPointL,Vertex3> (serialized),mapSurf2Curv(0),mapCurv2Surf(0) {
802 BuildBound();
803 if(verbosity>1)
804 cout << " -- End of serialized: mesure = " << mes << " border mesure " << mesb << endl;
805
806 if(nt > 0){
807 BuildAdj();
808 Buildbnormalv();
809 BuildjElementConteningVertex();
810 }
811 if(verbosity>1)
812 cout << " -- MeshL (serialized), d "<< 3 << ", n Edges " << nt << ", n Vtx "
813 << nv << " n Bord " << nbe << endl;
814 ffassert(mes>=0);
815 }
816
817
BuildCurvBasis()818 void MeshL::BuildCurvBasis(){
819 KN<R3> gx(nv), gy(nv), gz(nv);
820 for (int i=0 ; i<nv; i++) {
821 gx[i]=R3(0.,0.,0.);
822 gy[i]=R3(0.,0.,0.);
823 gz[i]=R3(0.,0.,0.);
824 }
825
826 for (int it=0 ; it<nt; it++) {
827 const EdgeL &K = elements[it];
828 R3 ie1, ie2;
829 ie1 = vertices[this->operator()(K[0])];
830 ie2 = vertices[this->operator()(K[1])];
831 for (int i = 0; i < EdgeL::nv; i++) {
832 int iiv = this->operator()(K[i]);
833 gx[iiv]+= K.mesure()*K.Edge(0);
834 gz[iiv]+=K.mesure()*K.NFrenet();
835 gy[iiv]+=K.mesure()*( K.NFrenet()^K.Edge(0) ); // sens ???
836 }
837 }
838 for (int i=0 ; i<nv; i++) {
839 gx[i]/=gx[i].norme();
840 gy[i]/=gy[i].norme();
841 gz[i]/=gz[i].norme();
842 if (verbosity>5)
843 cout << "NORMALIZE test covariant basis i: "<<i << " gx= " << gx[i] << " gy= " << gy[i] << " gz= " << gz[i] << endl;}
844 }
845
846
847 }
848