1 // ORIG-DATE: Dec 2007
2 // -*- Mode : c++ -*-
3 //
4 // SUMMARY : Model mesh 2d
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 namespace Fem2D
52 {
53
54 template<> int GenericMesh<TriangleS,BoundaryEdgeS,Vertex3>::kfind=0;
55 template<> int GenericMesh<TriangleS,BoundaryEdgeS,Vertex3>::kthrough=0;
56
57
58 const string GsbeginS="MeshS::GSave v0",GsendS="end";
GSave(FILE * ff,int offset) const59 void MeshS::GSave(FILE * ff,int offset) const
60 {
61 PlotStream f(ff);
62
63 f << GsbeginS ;
64 f << nv << nt << nbe;
65 for (int k=0; k<nv; k++) {
66 const Vertex & P = this->vertices[k];
67 f << P.x <<P.y << P.z << P.lab ;
68 }
69
70 for (int k=0; k<nt; k++) {
71 const Element & K(this->elements[k]);
72 int i0=this->operator()(K[0])+offset;
73 int i1=this->operator()(K[1])+offset;
74 int i2=this->operator()(K[2])+offset;
75 int lab=K.lab;
76 f << i0 << i1 << i2 << lab;
77 }
78
79 for (int k=0; k<nbe; k++) {
80 const BorderElement & K(this->borderelements[k]);
81 int i0=this->operator()(K[0])+offset;
82 int i1=this->operator()(K[1])+offset;
83 int lab=K.lab;
84 f << i0 << i1 << lab;
85 }
86 f << GsendS;
87 }
88
89
90
read(istream & f)91 void MeshS::read(istream &f)
92 { // read the mesh
93 int i;
94 // f >> nv >> nt >> neb ;
95 string str;
96 int err=0;
97 while(!f.eof())
98 {
99 f >> str;
100 //cout << str << endl;
101 if( str== "Vertices")
102 {
103 f >> nv;
104 assert(!this->vertices );
105 if(verbosity>2)
106 cout << " -- Nb of Vertex " << nv << endl;
107 this->vertices = new Vertex[nv];
108 for (i=0;i<nv;i++)
109 {
110 f >> this->vertices[i];
111 assert(f.good());
112 }
113 }
114 else if (str=="Triangles")
115 {
116 f >> nt;
117 assert(this->vertices && !this->elements);
118 this->elements = new Element [nt];
119 mes=0;
120 assert(this->elements);
121 if(verbosity>2)
122 cout << " -- Nb of Elements " << nt << endl;
123 for (int i=0;i<nt;i++)
124 {
125 this->t(i).Read1(f,this->vertices,this->nv);
126 if(this->t(i).mesure()<=0) err++; // Modif FH nov 2014
127 mes += this->t(i).mesure();
128 }
129 }
130 else if (str=="Edges")
131 {
132 mesb=0;
133 int kmv=0,ij;
134 f >> nbe;
135 assert(vertices);
136 this->borderelements = new BorderElement[nbe];
137 if(verbosity>2)
138 cout << " -- Nb of border Triangles " << nbe << endl;
139 for (i=0;i<nbe;i++)
140 {
141 this->be(i).Read1(f,this->vertices,this->nv);
142 mesb += this->be(i).mesure();
143 for(int j=0;j<BorderElement::nv;++j)
144 if(!vertices[ij=this->be(i,j)].lab)
145 {
146 vertices[ij].lab=1;
147 kmv++;
148 }
149 }
150 }
151 else if(str[0]=='#')
152 {// on mange la ligne
153 int c;
154 while ( (c=f.get()) != '\n' && c != EOF)
155 //cout << c << endl;
156 ;
157 }
158 }
159 assert( (nt >= 0 || nbe>=0) && nv>0) ;
160 /* done at up level ...
161 BuildBound();
162
163 if(nt > 0){
164 BuildAdj();
165 Buildbnormalv();
166 BuildjElementConteningVertex();
167 }
168 */
169 if(err!=0)
170 {
171 cerr << " MeshS::read: sorry bad mesh. Number of negative Triangles " << err << endl;
172 this->~MeshS();
173 ffassert(0);
174 }
175 }
176
177
readmsh(ifstream & f,int offset)178 void MeshS::readmsh(ifstream & f,int offset)
179 {
180 int err=0;
181 f >> nv >> nt >> nbe;
182 if(verbosity>2)
183 cout << " GRead : nv " << nv << " " << nt << " " << nbe << endl;
184 this->vertices = new Vertex[nv];
185 this->elements = new Element [nt];
186 this->borderelements = new BorderElement[nbe];
187 for (int k=0; k<nv; k++) {
188 Vertex & P = this->vertices[k];
189 f >> P.x >>P.y >> P.z >> P.lab ;
190 }
191 mes=0.;
192 mesb=0.;
193
194 if(nt == 0) {
195 cerr << " A meshS type must have elements " << endl;
196 ffassert(0);exit(1);
197
198 }
199
200 for (int k=0; k<nt; k++) {
201 int i[3],lab;
202 Element & K(this->elements[k]);
203 f >> i[0] >> i[1] >> i[2] >> lab;
204 Add(i,3,offset);
205 K.set(this->vertices,i,lab);
206 mes += K.mesure();
207 err += K.mesure() <0;
208
209 }
210
211
212
213 for (int k=0; k<nbe; k++) {
214 int i[2],lab;
215 BorderElement & K(this->borderelements[k]);
216 f >> i[0] >> i[1] >> lab;
217 Add(i,2,offset);
218 K.set(this->vertices,i,lab);
219 mesb += K.mesure();
220
221 }
222 if(err!=0)
223 {
224 cerr << " MeshS::readmsh : sorry bad mesh. Number of negative Tri " << err << endl;
225 this->~MeshS();
226 ffassert(0);
227 }
228
229 }
230
231
MeshS(const string filename,double ridgeangledetection)232 MeshS::MeshS(const string filename, double ridgeangledetection)
233 :mapSurf2Vol(0),mapVol2Surf(0),meshL(0) {
234 int ok=load(filename);
235 if(verbosity) {
236 cout << "read meshS ok " << ok ;
237 cout << "surface Mesh, num Triangles:= " << nt << ", num Vertice:= " << nv << " num boundary Edges:= " << nbe << endl;
238 }
239
240 if (ok) {
241 ifstream f(filename.c_str());
242 if(!f) {
243 cerr << " -- MeshS::MeshS Erreur opening " << filename<<endl;ffassert(0);exit(1);}
244 if(verbosity>2)
245 cout << " -- MeshS: Read On file \"" <<filename<<"\""<< endl;
246 if(filename.rfind(".msh")==filename.length()-4)
247 readmsh(f,-1);
248 else
249 read(f);
250 }
251
252 BuildBound();
253 BuildAdj();
254 Buildbnormalv();
255 BuildjElementConteningVertex();
256
257
258 if(verbosity>2) cout << " -- End of read: MeshS mesure = " << mes << " border mesure " << mesb << endl;
259
260 if(verbosity) cout << " -- MeshS : "<<filename << ", space dimension "<< 3 << ", num Triangle elts " << nt << ", num Vertice "
261 << nv << " num Bordary elts " << nbe << endl;
262
263
264 ffassert(mes>=0); // add F. Hecht sep 2009.
265
266 }
267
268
269
270
load(const string & filename)271 int MeshS::load(const string & filename)
272 {
273 int bin;
274 int ver,inm,dim;
275 int lf=filename.size()+20;
276 KN<char> fileb(lf),filef(lf);
277 char *data = new char[filename.size()+1];
278 size_t ssize = filename.size()+1;
279 char *ptr;
280 char *pfile=data;
281 strncpy( data, filename.c_str(),ssize);
282 ptr = strstr(data,".mesh");
283 if( !ptr ){
284 strcpy(filef,filename.c_str());
285 strcpy(fileb,filef);
286 strcat(filef,".mesh");
287 strcat(fileb,".meshb");
288 if( (inm=GmfOpenMesh(pfile=fileb, GmfRead,&ver,&dim)) )
289 bin=true;
290 else if( (inm=GmfOpenMesh(pfile=filef, GmfRead,&ver,&dim)) )
291 bin=false;
292 else
293 if(verbosity>5){
294 cerr << " Erreur ouverture file " << (char *) fileb << " " << (char *) filef <<endl;
295 return 1;
296 }
297 }
298 else{
299 if( !(inm=GmfOpenMesh(data, GmfRead,&ver,&dim)) ){
300 if(verbosity>5)
301 cerr << " Erreur ouverture file " << (char *) data << endl;
302 return 1;
303 }
304 }
305 // data file is readed and the meshes are initilized
306 int nv=-1,nTet=-1,nTri=-1,nSeg=-1, nPt=-1;
307 nv = GmfStatKwd(inm,GmfVertices); // vertice
308 nTri= GmfStatKwd(inm,GmfTriangles); // triangles in case of volume mesh -> boundary element / in case of surface mesh -> element
309 nSeg=GmfStatKwd(inm,GmfEdges); // segment elements only present in surface mesh
310 nPt=0; //GmfStatKwd(inm,GmfEdges); // points border on border mesh, not available at this time
311
312 if (nTri>0 && nSeg>0 && nPt==0)
313 if(verbosity>1)
314 cout << "data file "<< pfile << " contains only a MeshS, the MeshL associated is created (whitout border points)." << endl;
315 if (nTri>0 && nSeg>0 && nPt>0)
316 if(verbosity>1) cout << "data file "<< pfile << " contains a MeshS and MeshL" << endl;
317 if(verbosity && !nTri && !nTet)
318 cerr << " WARNING!!! The mesh file just contains a set of vertices" << endl;
319
320
321
322 this->set(nv,nTri,nSeg);
323 nPoints=nPt;
324
325 if(nTri == 0) {
326 cerr << " A meshS type must have elements " << endl;
327 ffassert(0);exit(1);}
328
329 if(verbosity>1)
330 cout << " -- MeshS(load): "<< (char *) data << ", MeshVersionFormatted:= " << ver << ", space dimension:= "<< dim
331 << ", num Triangles elts:= " << nTri << ", num vertice:= " << nv << " num Edges boundaries:= " << nSeg << endl;
332
333 if(dim != 3) {
334 cerr << "Err dim == " << dim << " !=3 " <<endl;
335 return 2; }
336 if( nv<=0 && (nTri <=0 || nSeg <0) ) {
337 cerr << " missing data "<< endl;
338 return 3;
339 }
340 int iv[3],lab;
341 float cr[3];
342 int mxlab=0, mnlab=0;
343 // read vertices
344 GmfGotoKwd(inm,GmfVertices);
345 for(int i=0;i<this->nv;++i) {
346 if(ver<2) {
347 GmfGetLin(inm,GmfVertices,&cr[0],&cr[1],&cr[2],&lab);
348 vertices[i].x=cr[0];
349 vertices[i].y=cr[1];
350 vertices[i].z=cr[2];}
351 else
352 GmfGetLin(inm,GmfVertices,&vertices[i].x,&vertices[i].y,&vertices[i].z,&lab);
353 vertices[i].lab=lab;
354 mxlab= max(mxlab,lab);
355 mnlab= min(mnlab,lab);
356 }
357 // read triangles (element meshS)
358 if(mnlab==0 && mxlab==0 ) {
359 int kmv=0;
360 mes=0;
361 GmfGotoKwd(inm,GmfTriangles);
362 for(int i=0;i<nTri;++i) {
363 GmfGetLin(inm,GmfTriangles,&iv[0],&iv[1],&iv[2],&lab);
364 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv && iv[2]>0 && iv[2]<=nv);
365 for(int j=0;j<3;++j)
366 if(!vertices[iv[j]-1].lab) {
367 vertices[iv[j]-1].lab=1;
368 kmv++;
369 }
370 for (int j=0;j<3;++j) iv[j]--;
371 elements[i].set(vertices,iv,lab);
372 mes += elements[i].mesure();
373 }
374 if(kmv&& verbosity>1) cout << " Aucun label Hack (FH) ??? => 1 sur les triangle frontiere "<<endl;
375 }
376 else {
377 mes=0;
378 GmfGotoKwd(inm,GmfTriangles);
379 for(int i=0;i<nTri;++i) {
380 GmfGetLin(inm,GmfTriangles,&iv[0],&iv[1],&iv[2],&lab);
381 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv && iv[2]>0 && iv[2]<=nv);
382 for (int j=0;j<3;++j) iv[j]--;
383 elements[i].set(this->vertices,iv,lab);
384 mes += elements[i].mesure();
385 }
386 }
387 // read edges (boundary elements meshS)
388 mesb=0;
389 GmfGotoKwd(inm,GmfEdges);
390 for(int i=0;i<nSeg;++i) {
391 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
392 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv);
393 for (int j=0;j<2;j++) iv[j]--;
394 borderelements[i].set(this->vertices,iv,lab); //element
395 mesb += this->borderelements[i].mesure();
396 }
397
398 // the .mesh contains edges, Building the meshS
399 // for this, surface vertices must be extract of the original vertice list and a mapping must be created between surface and volume vertices
400 if (nPoints>0) {
401 meshL = new MeshL();
402 // Number of Vertex in the surface
403 meshL->mapSurf2Curv=new int[nv];
404 meshL->mapCurv2Surf=new int[nv]; // mapping to curve/surface vertices
405 for (int k=0; k<nv; k++) {
406 meshL->mapSurf2Curv[k]=-1;
407 meshL->mapCurv2Surf[k]=0;
408 }
409 // search Vertex on the surface
410 int nbv_curv=0;
411 for (int k=0; k<nbe; k++) {
412 const BorderElement & K(this->borderelements[k]);
413 for(int jj=0; jj<2; jj++) {
414 int i0=this->operator()(K[jj]);
415 if( meshL->mapSurf2Curv[i0] == -1 ) {
416 // the mapping v_num_surf -> new numbering / liste_v_num_surf[nbv_surf] -> the global num
417 meshL->mapSurf2Curv[i0] = nbv_curv;
418 meshL->mapCurv2Surf[nbv_curv]= i0;
419 nbv_curv++;
420 }
421 }
422 }
423 this->meshL->set(nbv_curv,nSeg,nPoints);
424 // save the surface vertices
425 for (int k=0; k<nbv_curv; k++) {
426 int k0 = meshL->mapCurv2Surf[k];
427 const Vertex & P = this->vertices[k0];
428 meshL->vertices[k].lab=P.lab;
429 meshL->vertices[k].x=P.x;
430 meshL->vertices[k].y=P.y;
431 meshL->vertices[k].z=P.z;
432 }
433 // read triangles and change with the surface numbering
434 int kmv=0;
435 meshL->mes=0;
436 GmfGotoKwd(inm,GmfEdges);
437 for(int i=0;i<nSeg;++i) {
438 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
439 for (int j=0;j<2;++j)
440 iv[j]=meshL->mapSurf2Curv[iv[j]-1];
441 for(int j=0;j<2;++j)
442 if(!meshL->vertices[iv[j]].lab) {
443 meshL->vertices[iv[j]].lab=1;
444 kmv++;
445 }
446 meshL->elements[i].set(meshL->vertices,iv,lab);
447 meshL->mes += meshL->elements[i].mesure();
448 }
449 // reading border points with the curv numbering not available at this moment
450 meshL->mesb=0;
451 /*GmfGotoKwd(inm,GmfEdges);
452 for(int i=0;i<nSeg;++i) {
453 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
454 for (int j=0;j<2;++j) iv[j]=meshS->mapVol2Surf[iv[j]-1];
455 // assert( iv[0]>0 && iv[0]<=nbv_surf && iv[1]>0 && iv[1]<=nbv_surf);
456 meshS->borderelements[i].set(meshS->vertices,iv,lab);
457 meshS->mesb += meshS->borderelements[i].mesure();
458 }*/
459 }
460 else
461 if(verbosity>5) cout << " use Th = buildBdMesh(Th) to build the curve mesh associated" << endl;
462
463
464 if(verbosity>1)
465 cout << " -- MeshS(load): "<< (char *) data << ", MeshVersionFormatted:= " << ver << ", space dimension:= "<< dim
466 << ", Triangle elts:= " << nt << ", num vertice:= " << nv << ", num edges boundaries:= " << nbe << endl;
467 if(verbosity>1 && meshL)
468 cout << " -- MeshS::MeshL(load): "<< (char *) data << ", MeshVersionFormatted:= " << ver << ", space dimension:= "<< dim
469 << ", Edges elts:= " << meshL->nt << ", num vertice:= " << meshL->nv << ", num border points:= " << meshL->nbe << endl;
470
471 GmfCloseMesh(inm);
472 delete[] data;
473 return 0; // OK
474 }
475
476
MeshS(const string filename,bool cleanmesh,bool removeduplicate,bool rebuildboundary,int orientation,double precis_mesh,double ridgeangledetection)477 MeshS::MeshS(const string filename, bool cleanmesh, bool removeduplicate, bool rebuildboundary, int orientation, double precis_mesh, double ridgeangledetection)
478 :mapSurf2Vol(0),mapVol2Surf(0),meshL(0) {
479
480
481 int ok=load(filename);
482 if(verbosity) {
483 cout << "read mesh ok " << ok << endl;
484 cout << ", nt " << nt << ", nv " << nv << " nbe: = " << nbe << endl;
485 }
486 if(ok)
487 {
488 ifstream f(filename.c_str());
489 if(!f) {
490 cerr << " -- MeshS: Erreur opening " << filename<<endl;ffassert(0);exit(1);}
491 if(verbosity>2)
492 cout << " -- MeshS: Read On file \"" <<filename<<"\""<< endl;
493 if(filename.rfind(".msh")==filename.length()-4)
494 readmsh(f,-1);
495 else
496 read(f);
497 }
498 if (cleanmesh) {
499 if(verbosity>3)
500 cout << "before clean meshS, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
501 clean_mesh(precis_mesh, nv, nt, nbe, vertices, elements, borderelements, removeduplicate, rebuildboundary, orientation);
502 if(verbosity>3)
503 cout << "after clean meshS, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
504 }
505
506 BuildBound();
507 BuildAdj();
508 Buildbnormalv();
509 BuildjElementConteningVertex();
510
511 // if not edges then build the edges - need access to the old adjacensce to build eges and rebuild the new adj
512 if (nbe==0) {
513 BuildBdElem();
514 delete [] TheAdjacencesLink;
515 delete [] BoundaryElementHeadLink;
516 TheAdjacencesLink=0;
517 BoundaryElementHeadLink=0;
518 BuildBound();
519 BuildAdj();
520 Buildbnormalv();
521 BuildjElementConteningVertex();
522 }
523 // else BuildMeshL(ridgeangledetection);
524
525 if(verbosity>2)
526 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
527 if(verbosity)
528 cout << " -- MeshS : "<<filename << ", d "<< 3 << ", n Tri " << nt << ", n Vtx "
529 << nv << " n Border Edges " << nbe << endl;
530 ffassert(mes>=0); // add F. Hecht sep 2009.
531 }
532
533
534
MeshS(FILE * f,int offset)535 MeshS::MeshS(FILE *f,int offset)
536 :mapSurf2Vol(0),mapVol2Surf(0),meshL(0)
537 {
538 GRead(f,offset);// remove 1
539 assert( (nt >= 0 || nbe>=0) && nv>0) ;
540 BuildBound();
541 if(verbosity>2)
542 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
543
544 BuildAdj();
545 Buildbnormalv();
546 BuildjElementConteningVertex();
547
548 // if not edges then build the edges - need access to the old adjacensce to build eges and rebuild the new adj
549 if (nbe==0) {
550 BuildBdElem();
551 delete [] TheAdjacencesLink;
552 delete [] BoundaryElementHeadLink;
553 TheAdjacencesLink=0;
554 BoundaryElementHeadLink=0;
555 BuildBound();
556 BuildAdj();
557 Buildbnormalv();
558 BuildjElementConteningVertex();
559 }
560
561 if(verbosity>2)
562 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
563
564 if(verbosity>1)
565 cout << " -- MeshS (File *), d "<< 3 << ", n Tri " << nt << ", n Vtx "
566 << nv << " n Bord " << nbe << endl;
567 ffassert(mes>=0); // add F. Hecht sep 2009.
568 }
569
570
hmin() const571 double MeshS::hmin() const
572 {
573 R3 Pinf(1e100,1e100,1e100),Psup(-1e100,-1e100,-1e100); // Extremite de la boite englobante
574 double hmin=1e10;
575
576 for (int ii=0;ii< this->nv;ii++) {
577 R3 P( vertices[ii].x, vertices[ii].y, vertices[ii].z);
578 Pinf=Minc(P,Pinf);
579 Psup=Maxc(P,Psup);
580 }
581
582 for (int k=0;k<this->nt;k++) {
583 for (int e=0;e<3;e++){
584 if( this->elements[k].lenEdge(e) < Norme2(Psup-Pinf)/1e9 ){
585 const TriangleS & K(this->elements[k]);
586 int iv[3];
587 for(int jj=0; jj <3; jj++)
588 iv[jj] = this->operator()(K[jj]);
589 if(verbosity>2) for (int eh=0;eh<3;eh++)
590 cout << "TriangleS: " << k << " edge : " << eh << " lenght "<< this->elements[k].lenEdge(eh) << endl;
591 if(verbosity>2) cout << " A triangleS with a very small edge was created " << endl;
592 return 1;
593 }
594 hmin=min(hmin,this->elements[k].lenEdge(e)); // calcul de .lenEdge pour un Mesh3
595 }
596 }
597 ffassert(hmin>Norme2(Psup-Pinf)/1e9);
598 return hmin;
599 }
600
601
602 // brute force method
GRead(FILE * ff,int offset)603 void MeshS::GRead(FILE * ff,int offset)
604 {
605 PlotStream f(ff);
606 string s;
607 f >> s;
608 ffassert( s== GsbeginS);
609 f >> nv >> nt >> nbe;
610 if(verbosity>2)
611 cout << " GRead : nv " << nv << " " << nt << " " << nbe << endl;
612 this->vertices = new Vertex[nv];
613 this->elements = new Element [nt];
614 this->borderelements = new BorderElement[nbe];
615 for (int k=0; k<nv; k++) {
616 Vertex & P = this->vertices[k];
617 f >> P.x >>P.y >> P.z >> P.lab ;
618 }
619 mes=0.;
620 mesb=0.;
621
622 if(nt == 0) {
623 cerr << " A meshS type must have elements " << endl;
624 ffassert(0);exit(1);}
625
626
627 for (int k=0; k<nt; k++) {
628 int i[3],lab;
629 Element & K(this->elements[k]);
630 f >> i[0] >> i[1] >> i[2] >> lab;
631 Add(i,3,offset);
632 K.set(this->vertices,i,lab);
633 mes += K.mesure();
634
635 }
636 for (int k=0; k<nbe; k++) {
637 int i[2],lab;
638 BorderElement & K(this->borderelements[k]);
639 f >> i[0] >> i[1] >> lab;
640 Add(i,2,offset);
641 K.set(this->vertices,i,lab);
642 mesb += K.mesure();
643
644 }
645 f >> s;
646 ffassert( s== GsendS);
647 }
648
649
Find(Rd P,R2 & Phat,bool & outside,const Element * tstart) const650 const MeshS::Element * MeshS::Find( Rd P, R2 & Phat,bool & outside,const Element * tstart) const //;
651
652 {
653 if(searchMethod != 2)
654 {
655 GenericDataFindBoundary<GMesh> * gdfb=Buildgdfb();
656 if(gdfb )
657 {
658 double l[3];
659 int loutside;// 0 inside, 1 out close, 2, out fare, , -1 inside
660 int itt =gdfb->Find(P,l,loutside);
661 outside=loutside;
662 Phat=R2(l+1);
663 Element &K=(this->elements)[itt];
664 if( verbosity > 9)
665 cout << " - Find "<< P << " -> " << K(Phat) << " outside " << loutside << " k= " << itt
666 << " dist =" << (P-K(Phat)).norme() << " :: " << Phat << endl;
667 return itt<0 ? 0: this->elements + itt; // outside
668
669 }
670 }
671 static int count =0;
672 if( verbosity && count++< 5 )
673 cerr << " MeshS::Find warning brute force to day " << endl;
674 double dmin2 = 1e200;
675 bool out = true;
676 int nopt =-1;// best triangle number
677 R2 Pha;
678 for (int i=0;i<nt;i++)
679 {
680 kthrough++;
681 const TriangleS & K(this->elements[i]);
682 R3 A(K[0]),B(K[1]),C(K[2]);
683
684 R3 n = Area3(A,B,C); // the normal n toA,B,C
685 // to build the vectorial area
686 R a[]={(Area3(P,B,C),n),(Area3(A,P,C),n),(Area3(A,B,P),n)};
687 // PB de projection on
688 R s=(a[0]+a[1]+a[2]);
689 ffassert( s>0);
690 R eps=(s)*1e-8;
691 int na[3],nn=0;
692 if( a[0] < -eps) na[nn++]=0;
693 if( a[1] < -eps) na[nn++]=1;
694 if( a[2] < -eps) na[nn++]=2;
695 double l[3]={a[0]/s,a[1]/s,a[2]/s};
696
697 if (nn==0)
698 {
699 // ok inside
700
701 R2 Ph(l+1);
702 R3 Q = K(Ph);
703 double d2 =R3(P,Q).norme2() ;
704 if( dmin2 > d2)
705 {
706 dmin2 = d2;
707 out = d2 < s*1e-2;
708 Pha = Ph;
709 nopt= i;
710 }
711
712 }
713 else // on border of K
714 {
715 // Project on edge i=i0,i1
716 int i = na[0];
717 int i0=( i+1)%3;
718 int i1=( i+2)%3;
719 R3 E[]={K[i0],K[i1]};
720 R3 AB(E[0],E[1]);
721 double lab2 = AB.norme2();
722 double eps= lab2*1e-6;
723 double le= min(1.,max(0.,(AB,R3(E[0],P))/lab2));
724 l[i]=0;
725 l[i1]=le;
726 l[i0]=1-le;
727 R2 Ph(l+1);
728 R3 Q = K(Ph);
729 double d2 =R3(P,Q).norme2() ;
730 if( dmin2 > d2)
731 {
732 nopt= i;
733 dmin2 = d2;
734 out = d2 < s*1e-6;
735 Pha = Ph;
736 }
737 }
738
739 }
740 outside=out;
741 Phat=Pha;
742 ffassert(nopt>=0);
743 if( verbosity > 9 )
744 {
745 Element &K =(this->elements)[nopt];
746 cout << " - Find B "<< P << " -> " << K(Phat) << " , " << outside << " k= " << nopt
747 << " dist =" << (P-K(Phat)).norme() << " :: " << Phat << " " << sqrt(dmin2) <<endl;
748 }
749
750 return this->elements + nopt; // outside
751
752 }
753
754
755
MeshS(int nnv,int nnt,int nnbe,Vertex3 * vv,TriangleS * tt,BoundaryEdgeS * bb,bool cleanmesh,bool removeduplicate,bool rebuildboundary,int orientation,double precis_mesh)756 MeshS::MeshS(int nnv, int nnt, int nnbe, Vertex3 *vv, TriangleS *tt, BoundaryEdgeS *bb, bool cleanmesh, bool removeduplicate, bool rebuildboundary, int orientation, double precis_mesh)
757 :mapVol2Surf(0),mapSurf2Vol(0),meshL(0)
758 {
759 nv = nnv;
760 nt = nnt;
761 nbe =nnbe;
762 vertices = vv;
763 elements = tt;
764 borderelements = bb;
765 mes=0.;
766 mesb=0.;
767
768 for (int i=0;i<nt;i++)
769 mes += this->elements[i].mesure();
770 for (int i=0;i<nbe;i++)
771 mesb += this->be(i).mesure();
772
773 if (cleanmesh) {
774 if(verbosity>5)
775 cout << "before clean meshS, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
776 clean_mesh(precis_mesh, nv, nt, nbe, vertices, elements, borderelements, removeduplicate, rebuildboundary, orientation);
777 if(verbosity>5)
778 cout << "after clean meshS, nv: " <<nv << " nt:" << nt << " nbe:" << nbe << endl;
779 }
780 BuildBound();
781 BuildAdj();
782 Buildbnormalv();
783 BuildjElementConteningVertex();
784 // if not edges then build the edges - need access to the old adjacensce to build eges and rebuild the new adj
785 if (nbe==0) {
786 if(verbosity>3)
787 cout << " building of boundary " << endl;
788 BuildBdElem();
789 delete [] TheAdjacencesLink;
790 delete [] BoundaryElementHeadLink;
791 TheAdjacencesLink=0;
792 BoundaryElementHeadLink=0;
793 BuildBound();
794 BuildAdj();
795 Buildbnormalv();
796 BuildjElementConteningVertex();
797 }
798
799
800 if(verbosity>1)
801 cout << " -- End of read meshS: mesure = " << mes << " border mesure " << mesb << endl;
802
803 assert(mes>=0.);
804 }
805
806
MeshS(const Serialize & serialized)807 MeshS::MeshS(const Serialize &serialized)
808 :GenericMesh<TriangleS,BoundaryEdgeS,Vertex3> (serialized),
809 mapVol2Surf(0), mapSurf2Vol(0), meshL(0)
810 {
811 BuildBound();
812 if(verbosity>1)
813 cout << " -- End of serialized: mesure = " << mes << " border mesure " << mesb << endl;
814
815
816 BuildAdj();
817 Buildbnormalv();
818 BuildjElementConteningVertex();
819
820
821 if(verbosity>1)
822 cout << " -- MeshS (serialized), d "<< 3 << ", n Vtx " << nv <<", n Tri " << nt << " n Bord " << nbe << endl;
823 ffassert(mes>=0); // add F. Hecht sep 2009.
824
825 }
826
827
Save(const string & filename) const828 int MeshS::Save(const string & filename) const
829 {
830 int ver = GmfDouble, outm;
831 if ( !(outm = GmfOpenMesh(filename.c_str(),GmfWrite,ver,3)) ) {
832 cerr <<" -- MeshS**::Save UNABLE TO OPEN :"<< filename << endl;
833 return(1);
834 }
835 float fx,fy,fz;
836 // write vertice (meshS)
837 GmfSetKwd(outm,GmfVertices,nv);
838 for (int k=0; k<nv; k++) {
839 const Vertex & P = vertices[k];
840 GmfSetLin(outm,GmfVertices,fx=P.x,fy=P.y,fz=P.z,P.lab);
841 }
842 // write triangles (meshS)
843 GmfSetKwd(outm,GmfTriangles,nt);
844 for (int k=0; k<nt; k++) {
845 const MeshS::Element & K(elements[k]);
846 int i0=this->operator()(K[0])+1;
847 int i1=this->operator()(K[1])+1;
848 int i2=this->operator()(K[2])+1;
849 int lab=K.lab;
850 GmfSetLin(outm,GmfTriangles,i0,i1,i2,lab);
851 }
852 // write edges (meshS)
853 GmfSetKwd(outm,GmfEdges,nbe);
854 for (int k=0; k<nbe; k++) {
855 const BorderElement & K(borderelements[k]);
856 int i0=this->operator()(K[0])+1;
857 int i1=this->operator()(K[1])+1;
858 int lab=K.lab;
859 GmfSetLin(outm,GmfEdges,i0,i1,lab);
860 }
861 GmfCloseMesh(outm);
862 return (0);
863 }
864
865
866 /* Serialize MeshS::serialize_withBorderMesh() const {
867
868 // structure here for MeshL...but not use now
869 ffassert(0);
870 }
871 */
872
873
874
875
flipSurfaceMeshS(int surface_orientation)876 void MeshS::flipSurfaceMeshS(int surface_orientation)
877 {
878 /* inverse the orientation of the surface if necessary*/
879 /* and control that all surfaces are oriented in the same way*/
880 int nbflip=0;
881 for (int i=0;i<this->nt;i++) {
882 double mes_triangleS= this->elements[i].mesure();
883
884 if( surface_orientation*mes_triangleS < 0){
885 const TriangleS &K(elements[i] );
886 int iv[3];
887
888 for (int j=0 ; j<3 ; j++)
889 iv[j] = this->operator()(K[j]);
890
891 int iv_temp=iv[1];
892 iv[1]=iv[2];
893 iv[2]=iv_temp;
894 this->elements[i].set( this->vertices, iv, K.lab ) ;
895 nbflip++;
896 }
897 }
898 assert(nbflip==0 || nbflip== this->nt);
899 }
900
901
902 // determine the boundary edge list for meshS
BuildBdElem(const double angle)903 void MeshS::BuildBdElem(const double angle) {
904
905 delete [] borderelements; // to remove the previous pointers
906 borderelements = new BoundaryEdgeS[3 * nt]; // 3 * nt upper evaluated
907
908 HashTable<SortArray<int, 2>, int> edgesI(3 * nt, nt);
909 int* AdjLink = new int[3 * nt];
910
911 int nbeS=0,nbiS=0,nk=0;
912 // Build edges from the triangle list
913 for (int i = 0; i < nt; i++)
914 for (int j = 0; j < 3; j++) {
915
916 int jt = j, it = ElementAdj(i, jt);
917 TriangleS &K(elements[i]); // current element
918
919 // True boundary edge -> no adjacence / on domain border
920 if ((it == i || it < 0)) {
921 int iv[2];
922 for(int ip=0;ip<2;ip++)
923 iv[ip] = this->operator () (K [TriangleS::nvedge[j][ip]]);
924 if(verbosity>15)
925 cout << " the edge " << iv[0] << " " << iv[1] << " is a boundary " << endl;
926 be(nbeS).set(vertices,iv,K.lab);
927 mesb += be(nbeS).mesure();
928 nbeS++;
929
930 }
931 // internal edge -- check angular and no manifold
932 else {
933 TriangleS &K_adj(elements[it]); //adjacence element
934 int iv[2];
935 for(int ip=0;ip<2;ip++)
936 iv[ip] = this->operator () (K [TriangleS::nvedge[j][ip]]);
937
938 SortArray<int, 2> key(iv[0], iv[1]);
939 typename HashTable<SortArray<int,2>,int>::iterator p= edgesI.find(key);
940 if (!p) { // 1st time the edge is seen
941 // unit normal
942 R3 Normal = K.NFrenetUnitaire();
943 R3 Normal_adj = K_adj.NFrenet();
944 Normal_adj /= Normal_adj.norme();
945 R pdt = (Normal,Normal_adj); // scalar product
946 pdt = acos(pdt); // radian angle (Normal,Normal_adj)
947 if(verbosity>15)
948 cout << "Element num: " << i << " N " << Normal << " Element adjacent num: " << it << " N_adj " << Normal_adj << " angle between N N_adj = " << pdt <<endl;
949
950 if(pdt >= angle) {
951 if(verbosity>15)
952 cout << " the edge " <<nbeS <<": [" << iv[0] << " " << iv[1] << "] is a boundary with the angular criteria" << endl;
953 int lab = min(K.lab, K_adj.lab);
954 //cout << "oldnum" <<liste_v_num_surf[iv[0]] <<" /" << liste_v_num_surf[iv[1]] << endl;
955
956 be(nbeS).set(vertices,iv,lab);
957 mesb += be(nbeS).mesure();
958 edgesI.add(key, nbeS++);
959 }
960 }
961 // the edge is internal --- manifold or no?
962 else {
963 // the edge is adjacent to more 2 triangles -> no manifold -> it's a boundary edge
964 if (p->v < 0) {
965 int nk1=-1-p->v;
966 int nk2= AdjLink[nk1];
967 int lab = min(K.lab, K_adj.lab);
968 be(nbeS).set(vertices,iv,lab);
969 mesb += be(nbeS).mesure();
970 nbeS++;
971 if(nk2>=0) { // firt time remove existing link ...
972 AdjLink[nk1]=-2;
973 AdjLink[nk2]=-2;// on no manifold border .
974 }
975 }
976 else {
977 AdjLink[nk]=p->v;
978 AdjLink[p->v]=nk;
979 p->v=-1-nk;
980 nbiS++;
981 }
982 }
983 }
984 nk++; // increment the total edge jump --- nt * 3
985
986 }
987 assert(nt*3==nk);
988 delete [] AdjLink;
989 // update the number of edges
990 nbe = nbeS;
991 if (verbosity>5) cout << " Building edges from mesh3 nbe: "<< nbeS << " nbi: " << nbiS << endl;
992
993 BuildBound();
994 delete []TheAdjacencesLink;
995 delete [] BoundaryElementHeadLink;
996 TheAdjacencesLink=0;
997 BoundaryElementHeadLink=0;
998 BuildAdj();
999 Buildbnormalv();
1000 BuildjElementConteningVertex();
1001 }
1002
1003
BuildMeshL(double angle)1004 void MeshS::BuildMeshL(double angle)
1005 {
1006
1007 if (meshL) {
1008 cout << "error, MeshS::meshL previously created " << endl;
1009 return;
1010 }
1011 if (verbosity) cout << "Build meshL from meshS.... " << endl;
1012
1013
1014 int mes = 0, mesb = 0;
1015
1016 int *v_num_curve, *map_v_num_curve;
1017 // Extraction of Vertex belongs to the surface
1018 v_num_curve=new int[nv];
1019 map_v_num_curve=new int[nv];
1020 for (int k=0; k<nv; k++){
1021 v_num_curve[k]=-1;
1022 map_v_num_curve[k]=0;
1023 }
1024 // Search Vertex on the surface
1025 int nbv_curve=0;
1026 for (int k=0; k<nbe; k++) {
1027 const BoundaryEdgeS & K(borderelements[k]);
1028 for(int jj=0; jj<2; jj++){
1029 int i0=this->operator()(K[jj]);
1030 if( v_num_curve[i0] == -1 ){
1031 v_num_curve[i0] = nbv_curve;
1032 map_v_num_curve[nbv_curve]= i0;
1033 nbv_curve++;
1034 }
1035 }
1036 }
1037
1038 // set the curve vertex in meshL
1039 ffassert(nbv_curve);
1040
1041 Vertex3 *vL = new Vertex3[nbv_curve];
1042 EdgeL *tL = new EdgeL[nbe];
1043 EdgeL *ttL = tL;
1044
1045
1046 for (int iv=0; iv<nbv_curve; iv++) {
1047 int k0 = map_v_num_curve[iv];
1048 const Vertex3 & P = vertices[k0];
1049 vL[iv].x=P.x;
1050 vL[iv].y=P.y;
1051 vL[iv].z=P.z;
1052 vL[iv].lab=P.lab;
1053 }
1054
1055 ffassert(nbe);
1056
1057
1058 for (int it=0; it<nbe; it++) {
1059 int iv[2];
1060 const BoundaryEdgeS & K(borderelements[it]);
1061 for (int j=0;j<2;++j)
1062 iv[j]=v_num_curve[this->operator()(K[j])];
1063 int lab=K.lab;
1064 (ttL)->set(vL,iv,lab);
1065 mes += ttL++->mesure();
1066 }
1067
1068 // first building without list edges
1069
1070 meshL = new MeshL(nbv_curve,nbe,0,vL,tL,0);
1071 meshL->mapSurf2Curv = new int[nv];
1072 meshL->mapCurv2Surf= new int[nv];
1073 for(int i=0 ; i<nv ; i++) {
1074 meshL->mapSurf2Curv[i] = v_num_curve[i];
1075 meshL->mapCurv2Surf[i] = map_v_num_curve[i];
1076 }
1077 meshL->BuildBdElem(angle);
1078 meshL->BuildGTree();
1079 delete [] v_num_curve;
1080 delete [] map_v_num_curve;
1081
1082 }
1083
1084
1085
MeshS(const Serialize & serialized,int withSurface)1086 MeshS::MeshS(const Serialize &serialized, int withSurface)
1087 {
1088
1089 nt=0;nv=0;nbe=0;
1090 mes=0.;mesb=0.;
1091 vertices=0;elements=0;
1092 borderelements=0;bnormalv=0;
1093 TheAdjacencesLink=0;BoundaryElementHeadLink=0;
1094 ElementConteningVertex=0; gtree=0;
1095 gdfb=0;
1096 ffassert(withSurface==1);
1097 //const int nvTet = Tet::nv;
1098 const int nvTriangle = TriangleS::nv;
1099 const int nvEdge = BoundaryEdgeS::nv;
1100 const int nvPoint = BoundaryPointL::nv;
1101
1102 const int d = Rd::d;
1103 int havebordermesh(0);
1104 int dd,nnvTriangle,nnvEdge,nnvPoint,nnvL,nnbeL,nnv,nnt,nnbe;
1105 long long l=0;
1106 size_t pp=0;
1107 serialized.get(pp, l);
1108 serialized.get(pp,dd);
1109 serialized.get(pp,havebordermesh);
1110 serialized.get(pp,nnvTriangle);
1111 serialized.get(pp,nnvEdge);
1112 serialized.get(pp,nnvPoint);
1113 serialized.get(pp,nnt);
1114 serialized.get(pp,nnv);
1115 serialized.get(pp,nnbe);
1116 // miss serialized.get(pp,nadjnomanifold);
1117 serialized.get(pp,nnvL);
1118 serialized.get(pp,nnbeL);
1119
1120 ffassert(d==dd && nvTriangle == nnvTriangle && nvEdge == nnvEdge && nvPoint == nnvPoint && havebordermesh==1);
1121
1122 set(nnv,nnt,nnbe);
1123 for (int i=0;i<nv;++i)
1124 {
1125 for(int j=0;j<d;++j)
1126 serialized.get(pp,vertices[i][j]);
1127 serialized.get(pp,vertices[i].lab);
1128 }
1129 mes=0.;
1130 for (int i=0;i<nt;++i)
1131 {
1132 int ii[nvTriangle],lab;
1133 for(int j=0;j<nvTriangle;++j)
1134 serialized.get(pp,ii[j]);
1135 serialized.get(pp,lab);
1136 mes += elements[i].set(vertices,ii,lab).mesure();
1137
1138 }
1139 mesb=0;
1140 for (int i=0;i<nbe;++i)
1141 {
1142 int ii[nvEdge],lab;
1143 for(int j=0;j<nvEdge;++j)
1144 serialized.get(pp,ii[j]);
1145 serialized.get(pp, lab);
1146 mesb += borderelements[i].set(vertices,ii,lab).mesure();
1147 }
1148 // build the meshS
1149 meshL = new MeshL();
1150 meshL->nt=0;meshL->nv=0;meshL->nbe=0;
1151 meshL->mes=0.;meshL->mesb=0.;
1152 meshL->vertices=0;meshL->elements=0;
1153 meshL->borderelements=0;meshL->bnormalv=0;
1154 meshL->TheAdjacencesLink=0;meshL->BoundaryElementHeadLink=0;
1155 meshL->ElementConteningVertex=0; meshL->gtree=0;
1156
1157 meshL->set(nnvL,nbe,nnbeL);
1158
1159 // Number of Vertex in the surface
1160 meshL->mapSurf2Curv=new int[nv];
1161 meshL->mapCurv2Surf=new int[nv];
1162 for (int k=0; k<nv; ++k) {
1163 serialized.get(pp,meshL->mapCurv2Surf);
1164 serialized.get(pp,meshL->mapSurf2Curv[k]);
1165 }
1166
1167 for (int k=0; k<nnvL; ++k) {
1168 int k0 = meshL->mapSurf2Curv[k];
1169 const Vertex & P = this->vertices[k0];
1170 meshL->vertices[k].lab=P.lab;
1171 meshL->vertices[k].x=P.x;
1172 meshL->vertices[k].y=P.y;
1173 meshL->vertices[k].z=P.z;
1174 }
1175 int iv[nvEdge];
1176 for(int i=0;i<nbe;++i) {
1177 const BorderElement & K(borderelements[i]);
1178 for (int j=0;j<nvEdge;++j)
1179 iv[j]=meshL->mapSurf2Curv[this->operator()(K[j])];
1180
1181 meshL->elements[i].set(meshL->vertices,iv,K.lab);
1182 meshL->mes += meshL->elements[i].mesure();
1183 }
1184
1185 for (int i=0;i<nnbeL;++i) {
1186 int iv[nvPoint],lab;
1187 for(int j=0;j<nvPoint;++j)
1188 serialized.get(pp,iv[j]);
1189 serialized.get(pp, lab);
1190 meshL->borderelements[i].set(meshL->vertices,iv,lab);
1191 meshL->mesb += meshL->borderelements[i].mesure();
1192 }
1193 // check pp reading
1194 assert(pp==serialized.size());
1195
1196 // end building for mesh3
1197 BuildBound();
1198 if(verbosity>1)
1199 cout << " -- End of serialized: mesure = " << mes << " border mesure " << mesb << endl;
1200
1201 if(nt > 0){
1202 BuildAdj();
1203 Buildbnormalv();
1204 BuildjElementConteningVertex();
1205 }
1206
1207 if(verbosity>1)
1208 cout << " -- MeshS (serialized), d "<< 3 << ", n Tri " << nt << ", n Vtx "
1209 << nv << " n Bord " << nbe << endl;
1210 ffassert(mes>=0);
1211 // end building meshS
1212 meshL->BuildBound();
1213 if(verbosity>1)
1214 cout << " -- End of serialized MeshS:MeshL: mesure = " << meshL->mes << " border mesure " << meshL->mesb << endl;
1215
1216 meshL->BuildAdj();
1217 meshL->Buildbnormalv();
1218 meshL->BuildjElementConteningVertex();
1219
1220 if(verbosity>1)
1221 cout << " -- MeshS:MeshL (serialized), d "<< 3 << ", n Edge " << nt << ", n Vtx "
1222 << nv << " n Bord " << nbe << endl;
1223 ffassert(meshL->mes>=0);
1224
1225 }
1226
1227 // this function merges serialize2 in serialize1
serialize_withBorderMesh() const1228 Serialize MeshS::serialize_withBorderMesh() const {
1229
1230 ffassert(meshL);
1231
1232 const int nvTriangle = TriangleS::nv;
1233 const int nvEdge = BoundaryEdgeS::nv;
1234 const int nvPoint = BoundaryPointL::nv;
1235 const int d = Rd::d;
1236 long long l=0;
1237 int nvL=meshL->nv;
1238 int nbeL=meshL->nbe;
1239 int havebordermesh(1);
1240
1241 l += sizeof(long long);
1242 l += 10*sizeof(int); //( +nv nt nbe bordermesh(=1))
1243 l += nt*(sizeof(int)*(nvTriangle + 1));
1244 l += nv*( sizeof(int) + sizeof(double)*d);
1245 l += nbe*(sizeof(int)*(nvEdge+1));
1246 // add new2old mapping
1247 l += nv*(sizeof(int)); // mapSurf2Vol
1248 l += nv*(sizeof(int)); // mapVol2Surf
1249 l += nbeL*(sizeof(int)*(nvPoint+1));
1250
1251 if(verbosity>3)
1252 cout << "Serialize gmesh " << l << " " << nvTriangle << " " << nvEdge << " " << nvPoint << endl;
1253 Serialize serialized(l,GenericMesh_magicmesh);
1254 // cout << l << magicmesh << endl;
1255 size_t pp=0;
1256 serialized.put(pp, l);
1257 serialized.put(pp,d);
1258 serialized.put(pp,havebordermesh);
1259 serialized.put(pp,nvTriangle);
1260 serialized.put(pp,nvEdge);
1261 serialized.put(pp,nvPoint);
1262 serialized.put(pp,nt);
1263 serialized.put(pp,nv);
1264 serialized.put(pp,nbe);
1265 // serialized.put(pp,nadjnomanifold);
1266 serialized.put(pp,nvL);
1267 serialized.put(pp,nbeL);
1268
1269
1270 if (verbosity>9)
1271 cout << " GenericMesh Serialized : " << l << " " << nt << " " << nv << " " << nbe << endl;
1272 for (int i=0;i<nv;i++)
1273 {
1274 for(int j=0;j<d;++j)
1275 serialized.put(pp,vertices[i][j]);
1276 serialized.put(pp,vertices[i].lab);
1277 }
1278 for (int i=0;i<nt;i++)
1279 {
1280
1281 const Element & K(elements[i]);
1282 for(int j=0;j<nvTriangle;++j)
1283 serialized.put(pp,(int) operator()(K[j]));
1284 serialized.put(pp, K.lab);
1285 }
1286 for (int i=0;i<nbe;i++)
1287 {
1288 const BorderElement & K(borderelements[i]);
1289 for(int j=0;j<nvEdge;++j)
1290 serialized.put(pp,(int) operator()(K[j]));
1291 serialized.put(pp, K.lab);
1292 }
1293 // copy mappings to the vertice of the surface
1294 int vp;
1295 for (int i=0;i<nv;i++) {
1296 vp=meshL->mapCurv2Surf[i];
1297 serialized.put(pp,vp);
1298 vp=meshL->mapSurf2Curv[i];
1299 serialized.put(pp,vp);
1300 }
1301 for (int i=0;i<nbeL;i++)
1302 {
1303 const BoundaryPointL & K(meshL->borderelements[i]);
1304 for(int j=0;j<nvPoint;++j)
1305 serialized.put(pp,(int) meshL->operator()(K[j]));
1306 serialized.put(pp, K.lab);
1307 }
1308 assert(pp==serialized.size());
1309 return serialized;
1310 }
1311
1312
BuildCurvBasis()1313 KNM<R3> MeshS::BuildCurvBasis(){
1314 KNM<R3> g(3,nv); // gx(0,:) gy(1,:) gz(2,:);
1315 for (int i=0 ; i<nv; i++) {
1316 g(0,i)=R3(0.,0.,0.);
1317 g(1,i)=R3(0.,0.,0.);
1318 g(2,i)=R3(0.,0.,0.);
1319 }
1320
1321 for (int it=0 ; it<nt; it++) {
1322 const TriangleS &K = elements[it];
1323 R3 ie1, ie2, ie3;
1324 ie1 = vertices[this->operator()(K[0])];
1325 ie2 = vertices[this->operator()(K[1])];
1326 ie3 = vertices[this->operator()(K[2])];
1327
1328 for (int i = 0; i < TriangleS::nv; i++) {
1329 int iiv = this->operator()(K[i]);
1330 g(0,iiv)+= K.mesure()*K.Edge(2);
1331 g(2,iiv)+=K.mesure()*K.NFrenet();
1332 g(1,iiv)+=K.mesure()*K.Edge(1); // sens ???
1333 }
1334 }
1335 for (int i=0 ; i<nv; i++) {
1336 g(0,i)/=g(0,i).norme();
1337 g(1,i)/=g(1,i).norme();
1338 g(2,i)/=g(2,i).norme();
1339 if (verbosity>5)
1340 cout << "NORMALIZE test covariant basis i: "<<i << " gx= " << g(0,i) << " gy= " << g(1,i) << " gz= " << g(2,i) << endl;}
1341 return g;
1342 }
1343
1344
1345
1346
1347
1348 }
1349