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 "ufunction.hpp"
36 #include "error.hpp"
37 #include "RNM.hpp"
38 #include "libmesh5.h"
39
40
41 #include "Mesh2dn.hpp"
42 // for plotStream (a change)
43 #include "Mesh3dn.hpp"
44 #include "MeshSn.hpp"
45 #include "MeshLn.hpp"
46 #include "rgraph.hpp"
47 #include "fem.hpp"
48 #include "PlotStream.hpp"
49
50 namespace Fem2D {
51
52 // static const int nvfaceTet[4][3] = {{3,2,1}, {0,2,3},{ 3,1,0},{ 0,1,2}} ;;
53 // static const int nvedgeTet[6][2] = { {0,1},{0,2},{0,3},{1,2},{1,3},{2,3} };;
54
55 static const int nvfaceTria[1][3] = { {0,1,2} };
56 static const int nvedgeTria[3][2] = { {1,2},{2,0},{0,1}}; // tourne de le sens trigo donc Normal ext vect(1,0) ^ perp
57
58 //static const int nvfaceSeg[1][3] = {{-1,-1,1}};
59 static const int nvedgeSeg[1][2] = { {0,1} };
60
61 static const int nvadjSeg[2][1] = { {0},{1} };
62
63
64 template<> const int (* const GenericElement<DataSeg2>::nvface)[3] = 0 ;
65 template<> const int (* const GenericElement<DataSeg2>::nvedge)[2] = nvedgeSeg; //nvedgeTria ;
66 template<> const int (* const GenericElement<DataSeg2>::nvadj)[1] = nvadjSeg ;
67
68
69 template<> const int (* const GenericElement<DataTriangle2>::nvface)[3] = nvfaceTria ;
70 template<> const int (* const GenericElement<DataTriangle2>::nvedge)[2] = nvedgeTria ;
71 template<> const int (* const GenericElement<DataTriangle2>::nvadj)[2] = nvedgeTria ;
72 template<> const int GenericElement<DataTriangle2>::nitemdim[4] = {3,3,1,0 } ;
73
74
75 static const int onWhatIsEdge2d[3][7] = { {0,1,3, 2,0,0, 0}, // edge 0
76 {3,0,1, 0,2,0, 0},
77 {1,3,0, 0,0,2, 0}};
78
79 template<>
80 const int (* const GenericElement<DataTriangle2>::onWhatBorder)[7] = onWhatIsEdge2d ;
81
82 template<> int GenericMesh<Triangle2,BoundaryEdge2,Vertex2>::kfind=0;
83 template<> int GenericMesh<Triangle2,BoundaryEdge2,Vertex2>::kthrough=0;
84
85
Mesh2(const char * filename)86 Mesh2::Mesh2(const char * filename)
87 { // read the mesh
88
89 int nt,nv,nbe;
90 int ok=load(filename);
91 if(ok)
92 {
93 ifstream f(filename);
94 if(!f) {cerr << "Mesh2::Mesh2 Erreur opening " << filename<<endl;exit(1);}
95 if(verbosity)
96 cout << " Read On file \"" <<filename<<"\""<< endl;
97 f >> nv >> nt >> nbe ;
98 this->set(nv,nt,nbe);
99 if(verbosity)
100 cout << " -- Nb of Vertex " << nv << " " << " Nb of Triangles " << nt
101 << " , Nb of border edges " << nbe << endl;
102 assert(f.good() && nt && nv) ;
103 for (int i=0;i<nv;i++)
104 {
105 f >> this->vertices[i];
106 assert(f.good());
107 }
108 mes=0;
109 for (int i=0;i<nt;i++)
110 {
111 this->t(i).Read1(f,this->vertices,nv);
112 mes += t(i).mesure();
113 }
114 mesb=0.;
115 for (int i=0;i<nbe;i++)
116 {
117 this->be(i).Read1(f,this->vertices,nv);
118 mesb += be(i).mesure();
119 }
120 }
121 BuildBound();
122 BuildAdj();
123 Buildbnormalv();
124 BuildjElementConteningVertex();
125
126 if(verbosity)
127 cout << " - mesh mesure = " << mes << " border mesure: " << mesb << endl;
128 }
129
load(const string & filename)130 int Mesh2::load(const string & filename)
131 {
132
133 int bin;
134 int ver,inm,dim;
135 int lf=filename.size()+20;
136 KN<char> fileb(lf),filef(lf);
137 char * pfile;
138 strcpy(filef,filename.c_str());
139 strcpy(fileb,filef);
140 strcat(filef,".mesh");
141 strcat(fileb,".meshb");
142 if( (inm=GmfOpenMesh(pfile=fileb, GmfRead,&ver,&dim)) )
143 bin=true;
144 else if( (inm=GmfOpenMesh(pfile=filef, GmfRead,&ver,&dim)) )
145 bin=false;
146 else
147 { cerr << " Erreur ouverture file " << (char *) fileb << " " << (char *) filef << endl;
148 return 1;
149 }
150 int nv,nt,neb;
151 nv = GmfStatKwd(inm,GmfVertices);
152 nt = GmfStatKwd(inm,GmfTriangles);
153 neb=GmfStatKwd(inm,GmfEdges);
154 this->set(nv,nt,neb);
155
156 if(verbosity)
157 cout << pfile <<": ver " << ver << ", d "<< dim << ", nt " << nt << ", nv " << nv << " nbe: = " << nbe << endl;
158 if(dim != Rd::d) {
159 cerr << "Err dim == " << dim << " != " << Rd::d << endl;
160 return 2; }
161 if( nv<=0 && nt <=0 ) {
162 cerr << " missing data "<< endl;
163 return 3;
164 }
165 int iv[4],lab;
166 float cr[3]={0,0,0};
167 // read vertices
168 GmfGotoKwd(inm,GmfVertices);
169 int mxlab=0;
170 int mnlab=0;
171 for(int i=0;i<nv;++i)
172 {
173 if(ver<2) {
174 GmfGetLin(inm,GmfVertices,&cr[0],&cr[1],&lab);
175 vertices[i].x=cr[0];
176 vertices[i].y=cr[1];
177 }
178 else
179 GmfGetLin(inm,GmfVertices,&vertices[i].x,&vertices[i].y,&lab);
180 vertices[i].lab=lab;
181
182 mxlab= max(mxlab,lab);
183 mnlab= min(mnlab,lab);
184 }
185
186
187 // /* read mesh triangles */
188 {
189 mes=0;
190 GmfGotoKwd(inm,GmfTriangles);
191 for(int i=0;i<nt;++i)
192 {
193 GmfGetLin(inm,GmfTriangles,&iv[0],&iv[1],&iv[2],&lab);
194 for (int j=0;j<3;++j)
195 iv[j]--;
196 this->elements[i].set(this->vertices,iv,lab);
197 mes += this->elements[i].mesure();
198 }
199
200 }
201
202
203 /* read mesh segement */
204 mesb=0;
205 GmfGotoKwd(inm,GmfEdges);
206 for(int i=0;i<nbe;++i)
207 {
208 GmfGetLin(inm,GmfEdges,&iv[0],&iv[1],&lab);
209 assert( iv[0]>0 && iv[0]<=nv && iv[1]>0 && iv[1]<=nv);
210 for (int j=0;j<2;j++) iv[j]--;
211 this->borderelements[i].set(vertices,iv,lab);
212 mesb += this->borderelements[i].mesure();
213 }
214
215 GmfCloseMesh(inm);
216 if(verbosity>1) cout << " mesure : "<< mes << " border mesure : " << mesb<< endl;
217 return(0); // OK
218
219 }
220
Save(const string & filename)221 int Mesh2::Save(const string & filename)
222 {
223 int ver = GmfDouble, outm;
224 if ( !(outm = GmfOpenMesh(filename.c_str(),GmfWrite,ver,2)) ) {
225 cerr <<" -- Mesh2::Save UNABLE TO OPEN :"<< filename << endl;
226 return(1);
227 }
228 double fx,fy;
229 GmfSetKwd(outm,GmfVertices,this->nv);
230 for (int k=0; k<nv; k++) {
231 const Vertex & P = this->vertices[k];
232 GmfSetLin(outm,GmfVertices,fx=P.x,fy=P.y,P.lab);
233 }
234
235 GmfSetKwd(outm,GmfTriangles,nt);
236 for (int k=0; k<nt; k++) {
237 const Element & K(this->elements[k]);
238 int i0=this->operator()(K[0])+1;
239 int i1=this->operator()(K[1])+1;
240 int i2=this->operator()(K[2])+1;
241 int lab=K.lab;
242 GmfSetLin(outm,GmfTriangles,i0,i1,i2,lab);
243 }
244
245 GmfSetKwd(outm,GmfEdges,nbe);
246 for (int k=0; k<nbe; k++) {
247 const BorderElement & K(this->borderelements[k]);
248 int i0=this->operator()(K[0])+1;
249 int i1=this->operator()(K[1])+1;
250 int lab=K.lab;
251 GmfSetLin(outm,GmfEdges,i0,i1,lab);
252 }
253
254 GmfCloseMesh(outm);
255 return (0);
256
257 }
258
259 const string Gsbegin="Mesh2::GSave v0",Gsend="end";
260
261 template<class Mesh>
GSave2(FILE * ff,const Mesh & Th)262 void GSave2(FILE * ff,const Mesh & Th)
263 {
264 PlotStream f(ff);
265
266 f << Gsbegin ;
267 int nbe=Th.nbBrdElmts();
268 f << Th.nv << Th.nt << nbe;
269 for (int k=0; k<Th.nv; k++) {
270 const typename Mesh::Vertex & P = Th(k);
271 f << P.x <<P.y << P.lab ;
272 }
273
274 for (int k=0; k<Th.nt; k++) {
275 const typename Mesh::Element & K(Th[k]);
276 int i0=Th(K[0]);
277 int i1=Th(K[1]);
278 int i2=Th(K[2]);
279
280 int lab=K.lab;
281 f << i0 << i1 << i2 << lab;
282 }
283
284
285
286 for (int k=0; k<nbe; k++) {
287 const typename Mesh::BorderElement & K(Th.be(k));
288 int i0=Th(K[0]);
289 int i1=Th(K[1]);
290 int lab=K.lab;
291 f << i0 << i1 << lab;
292 }
293 f << Gsend;
294 }
295
296 template void GSave2<Mesh>(FILE * ff,const Mesh & Th) ;
297
298
Mesh2(FILE * f)299 Mesh2::Mesh2(FILE *f)
300 {
301
302 GRead(f);
303 assert( (nt >= 0 || nbe>=0) && nv>0) ;
304 BuildBound();
305 if(verbosity>1)
306 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
307
308 if(nt > 0){
309 BuildAdj();
310 Buildbnormalv();
311 BuildjElementConteningVertex();
312 }
313
314 if(verbosity>1)
315 cout << " -- Mesh2 (File *), d "<< 2 << ", n Tet " << nt << ", n Vtx "
316 << nv << " n Bord " << nbe << endl;
317
318 }
319
GRead(FILE * ff)320 void Mesh2::GRead(FILE * ff)
321 {
322 PlotStream f(ff);
323 string s;
324 f >> s;
325 ffassert( s== Gsbegin);
326 f >> nv >> nt >> nbe;
327 if(verbosity>1)
328 cout << " GRead : nv " << nv << " " << nt << " " << nbe << endl;
329 this->vertices = new Vertex[nv];
330 this->elements = new Element [nt];
331 this->borderelements = new BorderElement[nbe];
332 for (int k=0; k<nv; k++) {
333 Vertex & P = this->vertices[k];
334 f >> P.x >>P.y >> P.lab ;
335 }
336 mes=0.;
337 mesb=0.;
338
339 if(nt != 0)
340 {
341
342 for (int k=0; k<nt; k++) {
343 int i[4],lab;
344 Element & K(this->elements[k]);
345 f >> i[0] >> i[1] >> i[2] >> lab;
346 K.set(this->vertices,i,lab);
347 mes += K.mesure();
348
349 }
350 }
351
352
353 for (int k=0; k<nbe; k++) {
354 int i[4],lab;
355 BorderElement & K(this->borderelements[k]);
356 f >> i[0] >> i[1] >> lab;
357 K.set(this->vertices,i,lab);
358 mesb += K.mesure();
359
360 }
361 f >> s;
362 ffassert( s== Gsend);
363 }
364
Mesh2(int nnv,int nnt,int nnbe,Vertex2 * vv,Triangle2 * tt,BoundaryEdge2 * bb)365 Mesh2::Mesh2(int nnv, int nnt, int nnbe, Vertex2 *vv, Triangle2 *tt, BoundaryEdge2 *bb)
366 {
367
368 nv = nnv;
369 nt = nnt;
370 nbe =nnbe;
371
372 vertices = vv;
373 elements = tt;
374 borderelements = bb;
375
376 mes=0.;
377 mesb=0.;
378
379 for (int i=0;i<nt;i++)
380 mes += this->elements[i].mesure();
381
382 for (int i=0;i<nbe;i++)
383 mesb += this->be(i).mesure();
384
385
386 if(verbosity>1)
387 cout << " -- End of read: mesure = " << mes << " border mesure " << mesb << endl;
388
389 }
390
391 }
392