ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_structure.cpp
Revision: 591
Committed: Tue Nov 4 20:10:38 2014 UTC (10 years, 6 months ago) by francois
File size: 14312 byte(s)
Log Message:
gestion des unites pour le mailleur structure

File Contents

# User Rev Content
1 chebbi 463 #include <stdio.h>
2     #include <math.h>
3     #include "mg_gestionnaire.h"
4     #include "mailleur3d_structure.h"
5     #include "mg_maillage.h"
6     #include "fem_maillage.h"
7     #include "fem_hexa20.h"
8     #include "mg_geometrie.h"
9 francois 576 #include "mg_sommet_noeud.h"
10     #include "mg_arete_element.h"
11     #include "mg_face_element.h"
12     #include "mg_volume_element.h"
13 chebbi 463 #include <vector>
14    
15    
16 francois 585 struct integrale
17     {
18     double ti;
19     double li;
20     };
21    
22 francois 591 MAILLEUR3D_STRUCTURE::MAILLEUR3D_STRUCTURE(double xmintmp, double xmaxtmp, double ymintmp, double ymaxtmp, double zmintmp, double zmaxtmp,int ntmp, int mtmp, int ktmp,MG_GESTIONNAIRE*gesttmp,int avecgeo,double unit):MAILLEUR(false),xmin(xmintmp),xmax(xmaxtmp),ymin(ymintmp),ymax(ymaxtmp),zmin(zmintmp),zmax(zmaxtmp),n(ntmp),m(mtmp),k(ktmp),gest(gesttmp),geovirtuel(avecgeo),densitevariable(false),unite(unit)
23 chebbi 463 {
24     }
25    
26 francois 591 MAILLEUR3D_STRUCTURE::MAILLEUR3D_STRUCTURE( double xmintmp, double xmaxtmp, double ymintmp, double ymaxtmp, double zmintmp, double zmaxtmp,double (*fx)(double),double (*fy)(double),double (*fz)(double),MG_GESTIONNAIRE* gesttmp,int avecgeo,double unit):MAILLEUR(false),xmin(xmintmp),xmax(xmaxtmp),ymin(ymintmp),ymax(ymaxtmp),zmin(zmintmp),zmax(zmaxtmp),gest(gesttmp),geovirtuel(avecgeo),enx(fx),eny(fy),enz(fz),densitevariable(true),unite(unit)
27 francois 585 {
28     }
29    
30 chebbi 463 MAILLEUR3D_STRUCTURE::MAILLEUR3D_STRUCTURE( MAILLEUR3D_STRUCTURE& mdd):MAILLEUR(mdd),
31     xmin(mdd.xmin),xmax(mdd.xmax),ymin(mdd.ymin),ymax(mdd.ymax),zmin(mdd.zmin),zmax(mdd.zmax),n(mdd.n),m(mdd.m),k(mdd.k),gest(mdd.gest)
32     {
33     }
34    
35     MAILLEUR3D_STRUCTURE::~MAILLEUR3D_STRUCTURE()
36     {
37     }
38    
39 francois 585
40     void MAILLEUR3D_STRUCTURE::discretise(std::vector<double> &tab,double min,double max,int &nb,double (*en)(double))
41     {
42     double t1=0;
43     double t2=1;
44     std::vector<integrale> tab_integrale;
45     double ti,tii,t;
46     double longueur_metrique=0.;
47     tii=t1;
48     double dt=1./pas;
49     do
50     {
51     ti=tii;
52     tii=ti+(t2-t1)*dt;
53     if (tii>t2) tii=t2;
54     t=0.7886751345*ti+0.2113248654*tii;
55     double facteur1=(max-min)/(*en)(t);
56     t=0.7886751345*tii+0.2113248654*ti;
57     double facteur2=(max-min)/(*en)(t);
58     longueur_metrique=longueur_metrique+0.5*(tii-ti)*(facteur1+facteur2);
59     integrale pas;
60     pas.ti=ti;
61     pas.li=longueur_metrique;
62     tab_integrale.insert(tab_integrale.end(),pas);
63     }
64     while (tii<t2);
65     nb=(int)floor(longueur_metrique);
66     if (longueur_metrique-floor(longueur_metrique)>0.5) nb++;
67     if (nb<1) nb=1;
68     double valeur_cible=longueur_metrique/nb;
69     tab.push_back(min);
70     int nbcree=0;
71     int nb_pas_integrale=tab_integrale.size();
72     for (int i=0;i<nb_pas_integrale;i++)
73     {
74     while ((tab_integrale[i].li>(nbcree+1)*valeur_cible) && (nbcree<nb-1))
75     {
76     double ti=tab_integrale[i].ti;
77     double tii=t2;
78     double li=0.;
79     if (i!=0) li=tab_integrale[i-1].li;
80     double lii=tab_integrale[i].li;
81     if (i!=nb_pas_integrale-1) tii=tab_integrale[i+1].ti;
82     t=ti+(tii-ti)/(lii-li)*((nbcree+1)*valeur_cible-li);
83     double x=min+t*(max-min);
84     tab.push_back(x);
85     nbcree++;
86     }
87     }
88     tab.push_back(max);
89    
90     }
91    
92    
93 francois 551 void MAILLEUR3D_STRUCTURE::maille(MG_GROUPE_TOPOLOGIQUE* mggt)
94 chebbi 463 {
95     affiche((char*)"Création d'un maillage MG structuré");
96     double intx, inty, intz;
97 francois 585 std::vector<double> tabx,taby,tabz;
98     if (densitevariable==false)
99     {
100     intx=fabs(xmax-xmin)/n;
101     inty=fabs(ymax-ymin)/m;
102     intz=fabs(zmax-zmin)/k;
103     }
104     else
105     {
106     discretise(tabx,xmin,xmax,n,enx);
107     discretise(taby,ymin,ymax,m,eny);
108     discretise(tabz,zmin,zmax,k,enz);
109     }
110 chebbi 463 MG_GEOMETRIE* geo= NULL;
111 francois 576 MG_SOMMET_NOEUD *som1=NULL;
112     MG_SOMMET_NOEUD *som2=NULL;
113     MG_SOMMET_NOEUD *som3=NULL;
114     MG_SOMMET_NOEUD *som4=NULL;
115     MG_SOMMET_NOEUD *som5=NULL;
116     MG_SOMMET_NOEUD *som6=NULL;
117     MG_SOMMET_NOEUD *som7=NULL;
118     MG_SOMMET_NOEUD *som8=NULL;
119     MG_ARETE_ELEMENT* are1=NULL;
120     MG_ARETE_ELEMENT* are2=NULL;
121     MG_ARETE_ELEMENT* are3=NULL;
122     MG_ARETE_ELEMENT* are4=NULL;
123     MG_ARETE_ELEMENT* are5=NULL;
124     MG_ARETE_ELEMENT* are6=NULL;
125     MG_ARETE_ELEMENT* are7=NULL;
126     MG_ARETE_ELEMENT* are8=NULL;
127     MG_ARETE_ELEMENT* are9=NULL;
128     MG_ARETE_ELEMENT* are10=NULL;
129     MG_ARETE_ELEMENT* are11=NULL;
130     MG_ARETE_ELEMENT* are12=NULL;
131     MG_FACE_ELEMENT* face1=NULL;
132     MG_FACE_ELEMENT* face2=NULL;
133     MG_FACE_ELEMENT* face3=NULL;
134     MG_FACE_ELEMENT* face4=NULL;
135     MG_FACE_ELEMENT* face5=NULL;
136     MG_FACE_ELEMENT* face6=NULL;
137     MG_VOLUME_ELEMENT* vol=NULL;
138     if (geovirtuel==true)
139     {
140     geo=new MG_GEOMETRIE((char*)"VIRTUEL",(char*)"VIRTUEL");
141 francois 591 geo->change_valeur_unite(unite);
142 francois 576 gest->ajouter_mg_geometrie(geo);
143     som1=new MG_SOMMET_NOEUD(NULL);
144     geo->ajouter_mg_sommet(som1);
145     som2=new MG_SOMMET_NOEUD(NULL);
146     geo->ajouter_mg_sommet(som2);
147     som3=new MG_SOMMET_NOEUD(NULL);
148     geo->ajouter_mg_sommet(som3);
149     som4=new MG_SOMMET_NOEUD(NULL);
150     geo->ajouter_mg_sommet(som4);
151     som5=new MG_SOMMET_NOEUD(NULL);
152     geo->ajouter_mg_sommet(som5);
153     som6=new MG_SOMMET_NOEUD(NULL);
154     geo->ajouter_mg_sommet(som6);
155     som7=new MG_SOMMET_NOEUD(NULL);
156     geo->ajouter_mg_sommet(som7);
157     som8=new MG_SOMMET_NOEUD(NULL);
158     geo->ajouter_mg_sommet(som8);
159     are1=new MG_ARETE_ELEMENT();
160     geo->ajouter_mg_arete(are1);
161     are2=new MG_ARETE_ELEMENT();
162     geo->ajouter_mg_arete(are2);
163     are3=new MG_ARETE_ELEMENT();
164     geo->ajouter_mg_arete(are3);
165     are4=new MG_ARETE_ELEMENT();
166     geo->ajouter_mg_arete(are4);
167     are5=new MG_ARETE_ELEMENT();
168     geo->ajouter_mg_arete(are5);
169     are6=new MG_ARETE_ELEMENT();
170     geo->ajouter_mg_arete(are6);
171     are7=new MG_ARETE_ELEMENT();
172     geo->ajouter_mg_arete(are7);
173     are8=new MG_ARETE_ELEMENT();
174     geo->ajouter_mg_arete(are8);
175     are9=new MG_ARETE_ELEMENT();
176     geo->ajouter_mg_arete(are9);
177     are10=new MG_ARETE_ELEMENT();
178     geo->ajouter_mg_arete(are10);
179     are11=new MG_ARETE_ELEMENT();
180     geo->ajouter_mg_arete(are11);
181     are12=new MG_ARETE_ELEMENT();
182     geo->ajouter_mg_arete(are12);
183     face1=new MG_FACE_ELEMENT();
184     geo->ajouter_mg_face(face1);
185     face2=new MG_FACE_ELEMENT();
186     geo->ajouter_mg_face(face2);
187     face3=new MG_FACE_ELEMENT();
188     geo->ajouter_mg_face(face3);
189     face4=new MG_FACE_ELEMENT();
190     geo->ajouter_mg_face(face4);
191     face5=new MG_FACE_ELEMENT();
192     geo->ajouter_mg_face(face5);
193     face6=new MG_FACE_ELEMENT();
194     geo->ajouter_mg_face(face6);
195     vol=new MG_VOLUME_ELEMENT();
196     geo->ajouter_mg_volume(vol);
197     }
198     MG_MAILLAGE* mai=new MG_MAILLAGE(geo);
199 chebbi 463 gest->ajouter_mg_maillage(mai);
200     std::vector<MG_NOEUD*> vecnod;
201     for (int l=0; l<=k ;l++)
202     {
203     for (int j=0; j<=m;j++)
204     {
205     for (int i=0; i<=n;i++)
206     {
207 francois 576 MG_ELEMENT_TOPOLOGIQUE *topo=NULL;
208     if (geovirtuel)
209     {
210     if ((i==0) && (j==0) && (l==0)) topo=som1;
211     else if ((i==0) && (j==m) && (l==0)) topo=som2;
212     else if ((i==0) && (j==m) && (l==k)) topo=som3;
213     else if ((i==0) && (j==0) && (l==k)) topo=som4;
214     else if ((i==n) && (j==0) && (l==0)) topo=som5;
215     else if ((i==n) && (j==m) && (l==0)) topo=som6;
216     else if ((i==n) && (j==m) && (l==k)) topo=som7;
217     else if ((i==n) && (j==0) && (l==k)) topo=som8;
218     else if ((j==0) && (l==0)) topo=are1;
219     else if ((j==m) && (l==0)) topo=are2;
220     else if ((j==m) && (l==k)) topo=are3;
221     else if ((j==0) && (l==k)) topo=are4;
222     else if ((i==0) && (l==0)) topo=are5;
223     else if ((i==n) && (l==0)) topo=are6;
224     else if ((i==n) && (l==k)) topo=are7;
225     else if ((i==0) && (l==k)) topo=are8;
226     else if ((i==0) && (j==0)) topo=are9;
227     else if ((i==n) && (j==0)) topo=are10;
228     else if ((i==n) && (j==m)) topo=are11;
229     else if ((i==0) && (j==m)) topo=are12;
230     else if (i==0) topo=face1;
231     else if (i==n) topo=face2;
232     else if (j==0) topo=face3;
233     else if (j==m) topo=face4;
234     else if (l==0) topo=face5;
235     else if (l==k) topo=face6;
236     else topo=vol;
237     }
238 francois 585 double x,y,z;
239     if (densitevariable==false)
240     {
241     z=zmin+(l*intz);
242     y=ymin+(j*inty);
243     x=xmin+(i*intx);
244     }
245     else
246     {
247     x=tabx[i];
248     y=taby[j];
249     z=tabz[l];
250     }
251 francois 576 MG_NOEUD* nou=mai->ajouter_mg_noeud(topo,x,y,z,MAILLEUR_AUTO);
252 chebbi 463 vecnod.push_back(nou);
253 francois 576 if (geovirtuel)
254     {
255     if ((i==0) && (j==0) && (l==0)) som1->change_mg_noeud(nou);
256     else if ((i==0) && (j==m) && (l==0)) som2->change_mg_noeud(nou);
257     else if ((i==0) && (j==m) && (l==k)) som3->change_mg_noeud(nou);
258     else if ((i==0) && (j==0) && (l==k)) som4->change_mg_noeud(nou);
259     else if ((i==n) && (j==0) && (l==0)) som5->change_mg_noeud(nou);
260     else if ((i==n) && (j==m) && (l==0)) som6->change_mg_noeud(nou);
261     else if ((i==n) && (j==m) && (l==k)) som7->change_mg_noeud(nou);
262     else if ((i==n) && (j==0) && (l==k)) som8->change_mg_noeud(nou);
263     }
264     }
265 chebbi 463 }
266     }
267 francois 576 if (geovirtuel)
268     {
269     for(int i=0; i<n;i++)
270     {
271     int j=0,l=0;
272     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
273     MG_NOEUD* nod2=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
274     mai->ajouter_mg_segment(are1,nod1,nod2,MAILLEUR_AUTO);
275     j=m,l=0;
276     MG_NOEUD* nod3=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
277     MG_NOEUD* nod4=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
278     mai->ajouter_mg_segment(are2,nod3,nod4,MAILLEUR_AUTO);
279     j=m,l=k;
280     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
281     MG_NOEUD* nod6=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
282     mai->ajouter_mg_segment(are3,nod5,nod6,MAILLEUR_AUTO);
283     j=0,l=k;
284     MG_NOEUD* nod7=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
285     MG_NOEUD* nod8=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
286     mai->ajouter_mg_segment(are4,nod7,nod8,MAILLEUR_AUTO);
287     }
288     for(int j=0; j<m;j++)
289     {
290     int i=0,l=0;
291     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
292     MG_NOEUD* nod2=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
293     mai->ajouter_mg_segment(are5,nod1,nod2,MAILLEUR_AUTO);
294     i=n,l=0;
295     MG_NOEUD* nod3=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
296     MG_NOEUD* nod4=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
297     mai->ajouter_mg_segment(are6,nod3,nod4,MAILLEUR_AUTO);
298     i=n,l=k;
299     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
300     MG_NOEUD* nod6=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
301     mai->ajouter_mg_segment(are7,nod5,nod6,MAILLEUR_AUTO);
302     i=0,l=k;
303     MG_NOEUD* nod7=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
304     MG_NOEUD* nod8=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
305     mai->ajouter_mg_segment(are8,nod7,nod8,MAILLEUR_AUTO);
306     }
307     for(int l=0; l<k;l++)
308     {
309     int i=0,j=0;
310     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
311     MG_NOEUD* nod2=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
312     mai->ajouter_mg_segment(are9,nod1,nod2,MAILLEUR_AUTO);
313     i=n,j=0;
314     MG_NOEUD* nod3=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
315     MG_NOEUD* nod4=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
316     mai->ajouter_mg_segment(are10,nod3,nod4,MAILLEUR_AUTO);
317     i=n,j=m;
318     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
319     MG_NOEUD* nod6=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
320     mai->ajouter_mg_segment(are11,nod5,nod6,MAILLEUR_AUTO);
321     i=0,j=m;
322     MG_NOEUD* nod7=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
323     MG_NOEUD* nod8=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
324     mai->ajouter_mg_segment(are12,nod7,nod8,MAILLEUR_AUTO);
325     }
326     for (int j=0;j<m;j++)
327     for (int l=0;l<k;l++)
328     {
329     int i=0;
330     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
331     MG_NOEUD* nod2=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
332     MG_NOEUD* nod3=vecnod[i+((n+1)*(j+1))+((l+1)*(m+1)*(n+1))];
333     MG_NOEUD* nod4=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
334     i=n;
335     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
336     MG_NOEUD* nod6=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
337     MG_NOEUD* nod7=vecnod[i+((n+1)*(j+1))+((l+1)*(m+1)*(n+1))];
338     MG_NOEUD* nod8=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
339     mai->ajouter_mg_quadrangle(face1,nod1,nod4,nod3,nod2,MAILLEUR_AUTO);
340     mai->ajouter_mg_quadrangle(face2,nod5,nod6,nod7,nod8,MAILLEUR_AUTO);
341     }
342     for (int i=0;i<n;i++)
343     for (int l=0;l<k;l++)
344     {
345     int j=0;
346     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
347     MG_NOEUD* nod2=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
348     MG_NOEUD* nod3=vecnod[(i+1)+((n+1)*j)+((l+1)*(m+1)*(n+1))];
349     MG_NOEUD* nod4=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
350     j=m;
351     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
352     MG_NOEUD* nod6=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
353     MG_NOEUD* nod7=vecnod[(i+1)+((n+1)*j)+((l+1)*(m+1)*(n+1))];
354     MG_NOEUD* nod8=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
355     mai->ajouter_mg_quadrangle(face3,nod1,nod2,nod3,nod4,MAILLEUR_AUTO);
356     mai->ajouter_mg_quadrangle(face4,nod7,nod6,nod5,nod8,MAILLEUR_AUTO);
357     }
358     for (int i=0;i<n;i++)
359     for (int j=0;j<m;j++)
360     {
361     int l=0;
362     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
363     MG_NOEUD* nod2=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
364     MG_NOEUD* nod3=vecnod[(i+1)+((n+1)*(j+1))+(l*(m+1)*(n+1))];
365     MG_NOEUD* nod4=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
366     l=k;
367     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
368     MG_NOEUD* nod6=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
369     MG_NOEUD* nod7=vecnod[(i+1)+((n+1)*(j+1))+(l*(m+1)*(n+1))];
370     MG_NOEUD* nod8=vecnod[i+((n+1)*(j+1))+(l*(m+1)*(n+1))];
371     mai->ajouter_mg_quadrangle(face5,nod1,nod4,nod3,nod2,MAILLEUR_AUTO);
372     mai->ajouter_mg_quadrangle(face6,nod5,nod6,nod7,nod8,MAILLEUR_AUTO);
373     }
374     }
375     for(int l=0; l<k;l++)
376 chebbi 463 {
377     for(int j=0; j<m;j++)
378     {
379     for(int i=0; i<n;i++)
380     {
381     MG_NOEUD* nod1=vecnod[i+((n+1)*j)+(l*(m+1)*(n+1))];
382     MG_NOEUD* nod2=vecnod[(i+1)+((n+1)*j)+(l*(m+1)*(n+1))];
383     MG_NOEUD* nod3=vecnod[(i+1)+(n+1)*(j+1)+(l*(m+1)*(n+1))];
384     MG_NOEUD* nod4=vecnod[i+(n+1)*(j+1)+(l*(m+1)*(n+1))];
385     MG_NOEUD* nod5=vecnod[i+((n+1)*j)+((l+1)*(m+1)*(n+1))];
386     MG_NOEUD* nod6=vecnod[(i+1)+((n+1)*j)+((l+1)*(m+1)*(n+1))];
387     MG_NOEUD* nod7=vecnod[(i+1)+(n+1)*(j+1)+((l+1)*(m+1)*(n+1))];
388     MG_NOEUD* nod8=vecnod[i+(n+1)*(j+1)+((l+1)*(m+1)*(n+1))];
389 francois 576 MG_HEXA *hexa=mai->ajouter_mg_hexa(vol,nod1,nod2,nod3,nod4,nod5,nod6,nod7,nod8,MAILLEUR_AUTO);
390     }
391 chebbi 463 }
392     }
393 francois 465 BOITE_3D b(xmin,ymin,zmin,xmax,ymax,zmax);
394     mai->change_param_structure(b,n,m,k);
395 chebbi 463 affiche((char*)"Création d'un maillage FEM structuré");
396     FEM_MAILLAGE* fem = new FEM_MAILLAGE(geo, mai,1);
397     gest->ajouter_fem_maillage(fem);
398     fem->construire();
399     }