ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/app/criaqoperation/src/criaqoperators.cpp
Revision: 662
Committed: Tue Mar 10 02:21:38 2015 UTC (10 years, 5 months ago) by sattarpa
File size: 108867 byte(s)
Log Message:
compare_sample_points added to CRIAQOPERATION, samplepoint inserted by type (removed-remained) added to mailleur2d_ins_noeud.

File Contents

# User Rev Content
1 sattarpa 554 #include <stdio.h>
2     #include "mg_file.h"
3     #include "criaqoperators.h"
4 sattarpa 596 #include "tpl_grille.h"
5     #include "ot_mathematique.h"
6     #include <math.h>
7     #include "tpl_octree.h"
8 sattarpa 609 #include "fem_maillage_outils.h"
9     #include "mg_maillage_outils.h"
10 sattarpa 646 #include <sstream>
11     #include <string.h>
12 sattarpa 554
13 sattarpa 596 CRIAQOPERATORS::CRIAQOPERATORS():MAILLEUR(false)
14 sattarpa 554 {
15     }
16 sattarpa 596 CRIAQOPERATORS::~CRIAQOPERATORS()
17 sattarpa 554 {
18     }
19 sattarpa 662 void CRIAQOPERATORS::samplepoints_compare(char* samplepoints1,char* samplepoints2,char* sortedsmplpnts)
20     {
21     std::vector<double*> lstsamplepoints1;
22     FILE *in1=fopen(samplepoints1,"rt");
23     if (in1==NULL) affiche((char*) "file is not available");
24     while(!feof(in1))
25     {
26     char chaine[500];
27     char* res=fgets(chaine,500,in1);
28     if(res!=NULL)
29     {
30     double x,y,z;
31     double q1,q2,q3;
32     double* pcoordbc=new double[6];
33     int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
34     if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
35     else if(nb==6)
36     {
37     pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
38     lstsamplepoints1.push_back(pcoordbc);
39     }
40     }
41     }
42     fclose(in1);
43     std::vector<double*> lstsamplepoints2;
44     FILE *in2=fopen(samplepoints2,"rt");
45     if (in2==NULL) affiche((char*) "file is not available");
46     while(!feof(in2))
47     {
48     char chaine[500];
49     char* res=fgets(chaine,500,in2);
50     if(res!=NULL)
51     {
52     double x,y,z;
53     double q1,q2,q3;
54     double* pcoordbc=new double[6];
55     int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
56     if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
57     else if(nb==6)
58     {
59     pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
60     lstsamplepoints2.push_back(pcoordbc);
61     }
62     }
63     }
64     fclose(in2);
65     cout<<lstsamplepoints1.size()<<endl;
66     cout<<lstsamplepoints2.size()<<endl;
67    
68     std::vector<double*> lstsortedsmplpnts;
69     for(std::vector<double*>::iterator itlstsamplepoints1=lstsamplepoints1.begin();itlstsamplepoints1!=lstsamplepoints1.end();itlstsamplepoints1++)
70     {
71     double* pcoordbc=new double[7];
72     int countsmpcomp=0;
73     double* xyzsmp1=(*itlstsamplepoints1);
74     for(std::vector<double*>::iterator itlstsamplepoints2=lstsamplepoints2.begin();itlstsamplepoints2!=lstsamplepoints2.end();itlstsamplepoints2++)
75     {
76     double* xyzsmp2=(*itlstsamplepoints2);
77     if(xyzsmp1[0]==xyzsmp2[0] && xyzsmp1[1]==xyzsmp2[1] && xyzsmp1[2]==xyzsmp2[2] && xyzsmp1[3]==xyzsmp2[3] && xyzsmp1[4]==xyzsmp2[4] && xyzsmp1[5]==xyzsmp2[5])
78     countsmpcomp++;
79     }
80     if(countsmpcomp==0)
81     {
82     pcoordbc[0]=0; pcoordbc[1]=xyzsmp1[0]; pcoordbc[2]=xyzsmp1[1]; pcoordbc[3]=xyzsmp1[2]; pcoordbc[4]=xyzsmp1[3]; pcoordbc[5]=xyzsmp1[4]; pcoordbc[6]=xyzsmp1[5];
83     lstsortedsmplpnts.push_back(pcoordbc);
84     }
85     else if (countsmpcomp>0)
86     {
87     pcoordbc[0]=1; pcoordbc[1]=xyzsmp1[0]; pcoordbc[2]=xyzsmp1[1]; pcoordbc[3]=xyzsmp1[2]; pcoordbc[4]=xyzsmp1[3]; pcoordbc[5]=xyzsmp1[4]; pcoordbc[6]=xyzsmp1[5];
88     lstsortedsmplpnts.push_back(pcoordbc);
89     }
90    
91    
92     }
93     cout<<lstsortedsmplpnts.size()<<endl;
94    
95    
96     FILE* insmp=fopen(sortedsmplpnts,"wt");
97     for(std::vector<double*>::iterator itlstsortedsmplpnts=lstsortedsmplpnts.begin();itlstsortedsmplpnts!=lstsortedsmplpnts.end();itlstsortedsmplpnts++)
98     {
99     double* xyzsmpsort=(*itlstsortedsmplpnts);
100     //cout<<xyzsmpsort[0]<<" "<<xyzsmpsort[1]<<endl;
101     fprintf(insmp,"%le %le %le %le %le %le %le \n",xyzsmpsort[0],xyzsmpsort[1],xyzsmpsort[2],xyzsmpsort[3],xyzsmpsort[4],xyzsmpsort[5],xyzsmpsort[6]);
102     }
103     fclose(insmp);
104    
105     }
106    
107 sattarpa 596 void CRIAQOPERATORS::meshdistance_compare(char* referencemagicfilename,int refmeshno,char* comparemagicfilename,int compmeshno,char* outputcomparefile)
108 sattarpa 554 {
109 sattarpa 596 MG_FILE refgest(referencemagicfilename);
110     MG_MAILLAGE* refmgmai;
111 sattarpa 600 if (refmeshno==0) refmgmai=refgest.get_mg_maillage(refmeshno); else refmgmai=refgest.get_mg_maillageid(refmeshno);
112 sattarpa 596 MG_FILE compgest(comparemagicfilename);
113     MG_MAILLAGE* compmgmai;
114 sattarpa 600 if (compmeshno==0) compmgmai=compgest.get_mg_maillage(compmeshno); else compmgmai=compgest.get_mg_maillageid(compmeshno);
115 sattarpa 658
116     //octree initialization
117 sattarpa 596 double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
118     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
119     LISTE_MG_NOEUD::iterator it;
120     for (MG_NOEUD* no=refmgmai->get_premier_noeud(it);no!=NULL;no=refmgmai->get_suivant_noeud(it))
121     {
122     if (no->get_x()<xmin) xmin=no->get_x();
123     if (no->get_y()<ymin) ymin=no->get_y();
124     if (no->get_z()<zmin) zmin=no->get_z();
125     if (no->get_x()>xmax) xmax=no->get_x();
126     if (no->get_y()>ymax) ymax=no->get_y();
127     if (no->get_z()>zmax) zmax=no->get_z();
128     }
129     for (MG_NOEUD* no=compmgmai->get_premier_noeud(it);no!=NULL;no=compmgmai->get_suivant_noeud(it))
130     {
131     if (no->get_x()<xmin) xmin=no->get_x();
132     if (no->get_y()<ymin) ymin=no->get_y();
133     if (no->get_z()<zmin) zmin=no->get_z();
134     if (no->get_x()>xmax) xmax=no->get_x();
135     if (no->get_y()>ymax) ymax=no->get_y();
136     if (no->get_z()>zmax) zmax=no->get_z();
137     lstnoeud.ajouter(no);
138     }
139     OT_VECTEUR_3D vecmin(xmin,ymin,zmin);OT_VECTEUR_3D vecmmax(xmax,ymax,zmax);
140     OT_VECTEUR_3D vec(vecmmax,vecmin);
141    
142     //double search_radius=0.001*(vec.get_longueur());
143     double search_radius=3.;
144     OT_VECTEUR_3D min(xmin,ymin,zmin); OT_VECTEUR_3D max(xmax,ymax,zmax); OT_VECTEUR_3D average=(min+max)/2.; OT_VECTEUR_3D lengthvec(min,max);
145     double length=sqrt(lengthvec*lengthvec); double bxr=1.1;
146     xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
147     xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);
148    
149     TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octree;
150     octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
151     for (MG_NOEUD* no=compmgmai->get_premier_noeud(it);no!=NULL;no=compmgmai->get_suivant_noeud(it))
152     octree.inserer(no);
153    
154 sattarpa 646 MG_SOLUTION* mgsol=new MG_SOLUTION(refmgmai,1,(char*)"distance.sol",1,(char*)"Normaldistance",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
155 sattarpa 596 refgest.ajouter_mg_solution(mgsol);
156     mgsol->change_legende(0,"Normaldistance");
157     int i=0;
158     for (MG_NOEUD* nop=refmgmai->get_premier_noeud(it);nop!=NULL;nop=refmgmai->get_suivant_noeud(it))
159     {
160     TPL_MAP_ENTITE<MG_NOEUD*> lstneiclose;
161     int neindnb=lstneiclose.get_nb();
162     while(neindnb==0)
163     {
164     octree.rechercher(nop->get_x(),nop->get_y(),nop->get_z(),search_radius,lstneiclose);
165     neindnb=lstneiclose.get_nb();
166     }
167     map<double,MG_NOEUD*,std::less<double> > lstdist;
168     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itnop;
169     for(MG_NOEUD* no=lstneiclose.get_premier(itnop);no!=NULL;no=lstneiclose.get_suivant(itnop))
170     {
171     double distns=sqrt(pow(nop->get_x()-no->get_x(),2)+pow(nop->get_y()-no->get_y(),2)+pow(nop->get_z()-no->get_z(),2));
172    
173     OT_VECTEUR_3D w_n_tot(0.,0.,0.);
174     double norme_w_n_tot=0.;
175     OT_VECTEUR_3D w2_n_tot(0.,0.,0.);
176     double norme_w2_n_tot=0.;
177     int nbtri=nop->get_lien_triangle()->get_nb();
178     for (int i=0;i<nbtri;i++)
179     {
180     MG_TRIANGLE* tri=(MG_TRIANGLE*)nop->get_lien_triangle()->get(i);
181    
182     double *xyzn1=tri->get_noeud1()->get_coord();
183     double *xyzn2=tri->get_noeud2()->get_coord();
184     double *xyzn3=tri->get_noeud3()->get_coord();
185    
186     OT_VECTEUR_3D vec1(xyzn1,xyzn3);
187     OT_VECTEUR_3D vec2(xyzn1,xyzn2);
188     OT_VECTEUR_3D n_tri=vec1&vec2; // Normale au triangle
189    
190     double w=n_tri.get_longueur(); // 2 fois l'aire du triangle (pour éviter une division par 2)
191     n_tri.norme();
192     OT_VECTEUR_3D w_n=w*n_tri;
193     double norme_w_n=w_n.get_longueur();
194     w_n_tot=w_n_tot+w_n;
195     norme_w_n_tot=norme_w_n_tot+norme_w_n;
196     }
197    
198     OT_VECTEUR_3D n=w_n_tot/norme_w_n_tot; // Normale à la surface au noeud
199     n.norme();
200     OT_VECTEUR_3D distvec(nop->get_coord(),no->get_coord());
201     double normaldistance=distvec*n;
202     double normaldistanceabs=fabs(normaldistance);
203     //cout<<"normaldistanceabs= "<<normaldistanceabs<<endl;
204     lstdist.insert(pair<double,MG_NOEUD*>(normaldistanceabs,nop));
205     }
206     map<double,MG_NOEUD*,std::less<double> >::iterator itlstdist=lstdist.begin();
207     double distnsclos=(*itlstdist).first;
208     MG_NOEUD* nodcls=(*itlstdist).second;
209     //cout<<"distnsclos= "<<distnsclos<<endl;
210     mgsol->ecrire(distnsclos,i,0,0);
211     //mgsol->ecrire(distnsclos,1,0,0);
212     i++;
213 sattarpa 554 }
214 sattarpa 596 refgest.enregistrer(outputcomparefile);
215     }
216 sattarpa 609 void CRIAQOPERATORS::deformed_correspond_mgmaiadd(char* magicfilename,int meshno,int numsol1,int numsol2,int numsol3,int numchamp1,int numchamp2,int numchamp3,double coef)
217     {
218     MG_FILE gest(magicfilename);
219     FEM_MAILLAGE* femmai;
220     if (meshno==0) femmai=gest.get_fem_maillage(meshno); else femmai=gest.get_fem_maillageid(meshno);
221     FEM_SOLUTION* sol1=gest.get_fem_solutionid(numsol1);
222     FEM_SOLUTION* sol2=gest.get_fem_solutionid(numsol2);
223     FEM_SOLUTION* sol3=gest.get_fem_solutionid(numsol3);
224     femmai->calcul_deforme(sol1,numchamp1,sol2,numchamp2,sol3,numchamp3);
225     if (coef>1e-16)
226     {
227     cout<<" Création du maillage deforme avec correspondence text file"<<endl;
228     MG_MAILLAGE* defmgmai=new MG_MAILLAGE(gest.get_mg_geometrie(0));
229     gest.ajouter_mg_maillage(defmgmai);
230     {
231     if (femmai->get_degre()==1)
232     {
233     //gest=femmai->get_mg_maillage()->get_gestionnaire();
234     //gest->ajouter_mg_maillage(this);
235     std::vector<unsigned long*> lstcores;//sasan
236     LISTE_FEM_NOEUD::iterator it;
237     for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
238     {
239     double x=no->get_x(coef);
240     double y=no->get_y(coef);
241     double z=no->get_z(coef);
242     MG_NOEUD *newno=new MG_NOEUD(no->get_lien_topologie(),x,y,z,DEFORME);
243     defmgmai->ajouter_mg_noeud(newno);
244     no->change_numero(newno->get_id());
245     unsigned long* correscad_defcad=new unsigned long[2];//sasan
246     correscad_defcad[0]=no->get_mg_element_maillage()->get_id();//sasan
247     correscad_defcad[1]=newno->get_id();//sasan
248     lstcores.push_back(correscad_defcad);//sasan
249     }
250     //sasan
251     FILE* in1=fopen("correspond_deform_mgmaill.txt","wt");
252 sattarpa 659 fprintf(in1,"MG_MAILLAGE No. : deformed_MAILLAGE No.\n");
253 sattarpa 609 fprintf(in1,"%lu %lu\n",femmai->get_mg_maillage()->get_id(),defmgmai->get_id());
254     for(std::vector<unsigned long*>::iterator itlstpinsert=lstcores.begin();itlstpinsert!=lstcores.end();itlstpinsert++)
255     {
256     unsigned long* xyzpins1;
257     xyzpins1=(*itlstpinsert);
258     //cout<<" xyzpins1[0]= "<<xyzpins1[0]<<" xyzpins1[1]= "<<xyzpins1[1]<<endl;
259     fprintf(in1,"%lu %lu\n",xyzpins1[0],xyzpins1[1]);
260     }
261     fclose(in1);
262     //sasan till
263     LISTE_FEM_ELEMENT1::iterator it1;
264     for (FEM_ELEMENT1 *ele=femmai->get_premier_element1(it1);ele!=NULL;ele=femmai->get_suivant_element1(it1))
265     {
266     MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
267     MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
268     MG_SEGMENT *seg=new MG_SEGMENT(ele->get_lien_topologie(),no1,no2,DEFORME);
269     defmgmai->ajouter_mg_segment(seg);
270     }
271     LISTE_FEM_ELEMENT2::iterator it2;
272     for (FEM_ELEMENT2 *ele=femmai->get_premier_element2(it2);ele!=NULL;ele=femmai->get_suivant_element2(it2))
273     {
274     if (ele->get_nb_fem_noeud()==3)
275     {
276     MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
277     MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
278     MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
279     defmgmai->ajouter_mg_triangle(ele->get_lien_topologie(),no1,no2,no3,DEFORME);
280     }
281     if (ele->get_nb_fem_noeud()==4)
282     {
283     MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
284     MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
285     MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
286     MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
287     defmgmai->ajouter_mg_quadrangle(ele->get_lien_topologie(),no1,no2,no3,no4,DEFORME);
288     }
289     }
290     LISTE_FEM_ELEMENT3::iterator it3;
291     for (FEM_ELEMENT3 *ele=femmai->get_premier_element3(it3);ele!=NULL;ele=femmai->get_suivant_element3(it3))
292     {
293     if (ele->get_nb_fem_noeud()==4)
294     {
295     MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
296     MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
297     MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
298     MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
299     defmgmai->ajouter_mg_tetra(ele->get_lien_topologie(),no1,no2,no3,no4,DEFORME);
300     }
301     if (ele->get_nb_fem_noeud()==8)
302     {
303     MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
304     MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
305     MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
306     MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
307     MG_NOEUD* no5=defmgmai->get_mg_noeudid(ele->get_fem_noeud(4)->get_numero());
308     MG_NOEUD* no6=defmgmai->get_mg_noeudid(ele->get_fem_noeud(5)->get_numero());
309     MG_NOEUD* no7=defmgmai->get_mg_noeudid(ele->get_fem_noeud(6)->get_numero());
310     MG_NOEUD* no8=defmgmai->get_mg_noeudid(ele->get_fem_noeud(7)->get_numero());
311     defmgmai->ajouter_mg_hexa(ele->get_lien_topologie(),no1,no2,no3,no4,no5,no6,no7,no8,DEFORME);
312     }
313     }
314    
315     }
316     }
317     }
318     cout<<" Enregistrement"<<endl;
319     gest.enregistrer(magicfilename);
320     }
321 sattarpa 650 void CRIAQOPERATORS::curvtur_crit_remvsmp(char* magicfilename,char* correspondfilename,char* gnifinspointfile,char* out_gnifinspointfile_removondefect,double search_radius,double curvdif_coef,int gmsh_affiche,double relet_search_rad,int relet_curvature)
322 sattarpa 609 {
323     MG_FILE gest(magicfilename);
324     map<unsigned long, unsigned long> coresspondndidlst;
325     FILE* in1=fopen(correspondfilename,"rt");
326     if (in1==NULL) cout<<"correspondfile is not available"<<endl;
327     char mess[500];
328     char *res=fgets(mess,500,in1);
329 sattarpa 655 if(feof(in1))
330     cout<<"correspondfile has not been well written"<<endl;
331 sattarpa 609 res=fgets(mess,500,in1);
332 sattarpa 655 if(feof(in1))
333     cout<<"correspondfile has not been well written"<<endl;
334 sattarpa 609 int cadmeshno;
335     int defcadmeshno;
336 sattarpa 659 int nb1=sscanf(mess,"%d %d",&cadmeshno,&defcadmeshno);
337 sattarpa 609 while (!feof(in1))
338     {
339     char chaine[500];
340 sattarpa 655 res=fgets(chaine,500,in1);
341     if(res!=NULL)
342     {
343 sattarpa 609 unsigned long cadndid;
344     unsigned long defcadndid;
345     int nb=sscanf(chaine,"%lu %lu",&cadndid,&defcadndid);
346     if (nb!=-1 && nb!=2) cout<<"Wrong file format"<<endl;
347     else if (nb==2)
348     {
349     coresspondndidlst.insert(pair<unsigned long, unsigned long> (defcadndid,cadndid) );
350     }
351 sattarpa 655 }
352 sattarpa 609 }
353     fclose(in1);
354 sattarpa 650
355     cout<<"cadmeshno= "<<cadmeshno<<" defcadmeshno= "<<defcadmeshno<<endl;
356 sattarpa 609 MG_MAILLAGE* cadmgmai;
357     if (cadmeshno==0) cadmgmai=gest.get_mg_maillage(cadmeshno); else cadmgmai=gest.get_mg_maillageid(cadmeshno);
358     MG_MAILLAGE* defcadmgmai;
359     if (defcadmeshno==0) defcadmgmai=gest.get_mg_maillage(defcadmeshno); else defcadmgmai=gest.get_mg_maillageid(defcadmeshno);
360 sattarpa 650
361     cout<<"octree initialization "<<endl;
362     //CAD octree initialization
363     LISTE_MG_NOEUD::iterator itcadndcurv;
364     double xmincad=1e308,ymincad=1e308,zmincad=1e308,xmaxcad=-1e308,ymaxcad=-1e308,zmaxcad=-1e308;
365     TPL_MAP_ENTITE<MG_NOEUD*> lstcadnoeud;
366     for(MG_NOEUD* cadndcurv=cadmgmai->get_premier_noeud(itcadndcurv);cadndcurv!=NULL;cadndcurv=cadmgmai->get_suivant_noeud(itcadndcurv))
367     {
368     if (cadndcurv->get_x()<xmincad) xmincad=cadndcurv->get_x();
369     if (cadndcurv->get_y()<ymincad) ymincad=cadndcurv->get_y();
370     if (cadndcurv->get_z()<zmincad) zmincad=cadndcurv->get_z();
371     if (cadndcurv->get_x()>xmaxcad) xmaxcad=cadndcurv->get_x();
372     if (cadndcurv->get_y()>ymaxcad) ymaxcad=cadndcurv->get_y();
373     if (cadndcurv->get_z()>zmaxcad) zmaxcad=cadndcurv->get_z();
374     lstcadnoeud.ajouter(cadndcurv);
375     }
376     OT_VECTEUR_3D vecmincad(xmincad,ymincad,zmincad);
377     OT_VECTEUR_3D vecmmaxcad(xmaxcad,ymaxcad,zmaxcad);
378     OT_VECTEUR_3D vecad(vecmmaxcad,vecmincad);
379     OT_VECTEUR_3D mincad(xmincad,ymincad,zmincad); OT_VECTEUR_3D maxcad(xmaxcad,ymaxcad,zmaxcad); OT_VECTEUR_3D averagecad=(mincad+maxcad)/2.; OT_VECTEUR_3D lengthvecad(mincad,maxcad);
380     double lengthcad=sqrt(lengthvecad*lengthvecad); double bxrcad=1.1;
381     xmincad=averagecad.get_x()-(lengthcad*bxrcad);ymincad=averagecad.get_y()-(lengthcad*bxrcad);zmincad=averagecad.get_z()-(lengthcad*bxrcad);
382     xmaxcad=averagecad.get_x()+(lengthcad*bxrcad);ymaxcad=averagecad.get_y()+(lengthcad*bxrcad);zmaxcad=averagecad.get_z()+(lengthcad*bxrcad);
383     TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreecad;
384     octreecad.initialiser(&lstcadnoeud,1,xmincad,ymincad,zmincad,xmaxcad,ymaxcad,zmaxcad);
385     for(MG_NOEUD* cadndcurv=cadmgmai->get_premier_noeud(itcadndcurv);cadndcurv!=NULL;cadndcurv=cadmgmai->get_suivant_noeud(itcadndcurv))
386     octreecad.inserer(cadndcurv);
387     //deformed_CAD octree initialization
388     LISTE_MG_NOEUD::iterator itdefcadndcurv;
389     double xmindefcad=1e308,ymindefcad=1e308,zmindefcad=1e308,xmaxdefcad=-1e308,ymaxdefcad=-1e308,zmaxdefcad=-1e308;
390     TPL_MAP_ENTITE<MG_NOEUD*> lstdefcadnoeud;
391     for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
392     {
393     if (defcadndcurv->get_x()<xmindefcad) xmindefcad=defcadndcurv->get_x();
394     if (defcadndcurv->get_y()<ymindefcad) ymindefcad=defcadndcurv->get_y();
395     if (defcadndcurv->get_z()<zmindefcad) zmindefcad=defcadndcurv->get_z();
396     if (defcadndcurv->get_x()>xmaxdefcad) xmaxdefcad=defcadndcurv->get_x();
397     if (defcadndcurv->get_y()>ymaxdefcad) ymaxdefcad=defcadndcurv->get_y();
398     if (defcadndcurv->get_z()>zmaxdefcad) zmaxdefcad=defcadndcurv->get_z();
399     lstdefcadnoeud.ajouter(defcadndcurv);
400     }
401     OT_VECTEUR_3D vecmindefcad(xmindefcad,ymindefcad,zmindefcad);OT_VECTEUR_3D vecmmaxdefcad(xmaxdefcad,ymaxdefcad,zmaxdefcad);
402     OT_VECTEUR_3D vecdefcad(vecmmaxdefcad,vecmindefcad);
403     OT_VECTEUR_3D mindefcad(xmindefcad,ymindefcad,zmindefcad); OT_VECTEUR_3D maxdefcad(xmaxdefcad,ymaxdefcad,zmaxdefcad); OT_VECTEUR_3D averagedefcad=(mindefcad+maxdefcad)/2.; OT_VECTEUR_3D lengthvecdefcad(mindefcad,maxdefcad);
404     double lengthdefcad=sqrt(lengthvecdefcad*lengthvecdefcad); double bxrdefcad=1.1;
405     xmindefcad=averagedefcad.get_x()-(lengthdefcad*bxrdefcad);ymindefcad=averagedefcad.get_y()-(lengthdefcad*bxrdefcad);zmindefcad=averagedefcad.get_z()-(lengthdefcad*bxrdefcad);
406     xmaxdefcad=averagedefcad.get_x()+(lengthdefcad*bxrdefcad);ymaxdefcad=averagedefcad.get_y()+(lengthdefcad*bxrdefcad);zmaxdefcad=averagedefcad.get_z()+(lengthdefcad*bxrdefcad);
407     TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreedefcad;
408     octreedefcad.initialiser(&lstdefcadnoeud,1,xmindefcad,ymindefcad,zmindefcad,xmaxdefcad,ymaxdefcad,zmaxdefcad);
409     for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
410     octreedefcad.inserer(defcadndcurv);
411    
412 sattarpa 609 multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_cadmgmai;
413     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_cadmgmai;
414 sattarpa 646 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_cadmgmai;
415 sattarpa 609 multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_defcadmgmai;
416     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_defcadmgmai;
417 sattarpa 646 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_defcadmgmai;
418 sattarpa 650 map<MG_NOEUD*,double > lstkmin_cad;
419     map<MG_NOEUD*,double > lstkmax_cad;
420     map<MG_NOEUD*,double > lstkgau_cad;
421     map<MG_NOEUD*,double > lstkmin_defcad;
422     map<MG_NOEUD*,double > lstkmax_defcad;
423     map<MG_NOEUD*,double > lstkgau_defcad;
424     multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_cadmgmai;
425     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_cadmgmai;
426     multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_cadmgmai;
427     multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_defcadmgmai;
428     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_defcadmgmai;
429     multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_defcadmgmai;
430 sattarpa 609 std::map<std::string,double> liste_kmin;
431     std::map<std::string,double> liste_kmax;
432     std::map<std::string,double> liste_gau;
433     std::map<std::string,double> liste_courburemax1;
434     std::map<std::string,double> liste_courburemax2;
435     std::map<std::string,double> liste_courburemax3;
436     std::map<std::string,double> liste_courburemax4;
437     std::map<std::string,double> liste_courburemax5;
438     int opensurafece=1;
439 sattarpa 650
440 sattarpa 609 cout<<"calcul courbure cadmgmai"<<endl;
441     // calcul courbure cadmgmai ********************************
442     for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv))
443     {
444     double kmin1=0.;
445     double kmax1=0.;
446     double kmin2=0.;
447     double kmax2=0.;
448     MG_MAILLAGE_OUTILS *cou;
449     cou->calcul_courbure_face_arete_sommet(cadmgnd,kmin1,kmax1,kmin2,kmax2,liste_kmin,liste_kmax,liste_gau,liste_courburemax1,liste_courburemax2,liste_courburemax3,liste_courburemax4,liste_courburemax5,opensurafece);
450     }
451     for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv))
452     {
453 sattarpa 650
454 sattarpa 609 unsigned long idface=cadmgnd->get_lien_topologie()->get_id();
455     unsigned long idnoeud=cadmgnd->get_id();
456     char mess[500];
457     sprintf(mess,"%lu_%lu",idnoeud,idface);
458     for(std::map<std::string,double>::iterator itkmin=liste_kmin.begin();itkmin!=liste_kmin.end();itkmin++)
459     {
460     if(mess==(*itkmin).first)
461 sattarpa 650 {
462 sattarpa 609 lstkmin_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,cadmgnd) );
463 sattarpa 650 lstkmin_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmin).second));
464     }
465 sattarpa 609 }
466     for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
467     {
468     if(mess==(*itkmax).first)
469 sattarpa 650 {
470 sattarpa 609 lstkmax_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,cadmgnd) );
471 sattarpa 650 lstkmax_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmax).second));
472     }
473 sattarpa 609 }
474 sattarpa 646 for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
475     {
476     if(mess==(*itkgau).first)
477 sattarpa 650 {
478 sattarpa 646 lstkgau_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,cadmgnd) );
479 sattarpa 650 lstkgau_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkgau).second));
480     }
481 sattarpa 646 }
482 sattarpa 609 }
483 sattarpa 650 //calc reletive curvature CAD
484     map<MG_NOEUD*,double >::iterator itlstkmin_cad;
485     map<MG_NOEUD*,double >::iterator itlstkmax_cad;
486     map<MG_NOEUD*,double >::iterator itlstkgau_cad;
487     for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv))
488     {
489     itlstkmin_cad=lstkmin_cad.find(cadmgnd);
490     double relet_kmincad=(*itlstkmin_cad).second;
491     itlstkmax_cad=lstkmax_cad.find(cadmgnd);
492     double relet_kmaxcad=(*itlstkmax_cad).second;
493     itlstkgau_cad=lstkgau_cad.find(cadmgnd);
494     double relet_kgaucad=(*itlstkgau_cad).second;
495     TPL_MAP_ENTITE<MG_NOEUD*> lstnei_relcad;
496     int relcad_neindnb=lstnei_relcad.get_nb();
497     double relet_search_rad1=relet_search_rad;
498     while(relcad_neindnb<1)
499     {
500     octreecad.rechercher(cadmgnd->get_x(),cadmgnd->get_y(),cadmgnd->get_z(),relet_search_rad1,lstnei_relcad);
501     relcad_neindnb=lstnei_relcad.get_nb();
502     relet_search_rad1=relet_search_rad1+0.1*relet_search_rad;
503     }
504     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
505     for(MG_NOEUD* cadndnei=lstnei_relcad.get_premier(itndnei);cadndnei!=NULL;cadndnei=lstnei_relcad.get_suivant(itndnei))
506     {
507     itlstkmin_cad=lstkmin_cad.find(cadndnei);
508     relet_kmincad=relet_kmincad+(*itlstkmin_cad).second;
509     itlstkmax_cad=lstkmax_cad.find(cadndnei);
510     relet_kmaxcad=relet_kmaxcad+(*itlstkmax_cad).second;
511     itlstkgau_cad=lstkgau_cad.find(cadndnei);
512     relet_kgaucad=relet_kgaucad+(*itlstkgau_cad).second;
513     }
514     relet_kmincad=relet_kmincad/lstnei_relcad.get_nb();
515     relet_kmaxcad=relet_kmaxcad/lstnei_relcad.get_nb();
516     relet_kgaucad=relet_kgaucad/lstnei_relcad.get_nb();
517     lstkmin_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmincad,cadmgnd));
518     lstkmax_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmaxcad,cadmgnd));
519     lstkgau_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kgaucad,cadmgnd));
520     }
521     // affichage reletive sur GMSH ******************************
522     if (gmsh_affiche>1)
523     {
524     cout<<"affiche: calcul_courbure_cadmgmai_relet "<<endl;
525     std::ostringstream oss;
526     oss << "calcul_courbure_cadmgmai_curvdif_relet_coef:" <<curvdif_coef ;
527     std::string namefic = oss.str();
528     char *cstr = new char[namefic.length() + 1];
529     strcpy(cstr, namefic.c_str());
530     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
531     gest.ajouter_mg_solution(sol1);
532     sol1->change_legende(0,"Kmin_relet");
533     sol1->change_legende(1,"Kmax_relet");
534     sol1->change_legende(2,"K_gau_relet");
535     LISTE_MG_TRIANGLE::iterator ittri;
536     int compte=0;
537     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
538     {
539     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
540     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
541     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
542     MG_NOEUD *no1=tri->get_noeud1();
543     MG_NOEUD *no2=tri->get_noeud2();
544     MG_NOEUD *no3=tri->get_noeud3();
545     double valeur0;
546     double valeur1;
547     double valeur2;
548     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
549     {
550     if (no1==(*itlstkmin_relt_cadmgmai).second)
551     valeur0=(*itlstkmin_relt_cadmgmai).first;
552     }
553     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
554     {
555     if (no1==(*itlstkmax_relt_cadmgmai).second)
556     valeur1=(*itlstkmax_relt_cadmgmai).first;
557     }
558     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
559     {
560     if (no1==(*itlstkgau_relt_cadmgmai).second)
561     valeur2=(*itlstkgau_relt_cadmgmai).first;
562     }
563     sol1->ecrire(valeur0,compte,0,0);
564     sol1->ecrire(valeur1,compte,1,0);
565     sol1->ecrire(valeur2,compte,2,0);
566     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
567     {
568     if (no2==(*itlstkmin_relt_cadmgmai).second)
569     valeur0=(*itlstkmin_relt_cadmgmai).first;
570     }
571     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
572     {
573     if (no2==(*itlstkmax_relt_cadmgmai).second)
574     valeur1=(*itlstkmax_relt_cadmgmai).first;
575     }
576     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
577     {
578     if (no2==(*itlstkgau_relt_cadmgmai).second)
579     valeur2=(*itlstkgau_relt_cadmgmai).first;
580     }
581     sol1->ecrire(valeur0,compte,0,1);
582     sol1->ecrire(valeur1,compte,1,1);
583     sol1->ecrire(valeur2,compte,2,1);
584     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
585     {
586     if (no3==(*itlstkmin_relt_cadmgmai).second)
587     valeur0=(*itlstkmin_relt_cadmgmai).first;
588     }
589     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
590     {
591     if (no3==(*itlstkmax_relt_cadmgmai).second)
592     valeur1=(*itlstkmax_relt_cadmgmai).first;
593     }
594     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
595     {
596     if (no3==(*itlstkgau_relt_cadmgmai).second)
597     valeur2=(*itlstkgau_relt_cadmgmai).first;
598     }
599     sol1->ecrire(valeur0,compte,0,2);
600     sol1->ecrire(valeur1,compte,1,2);
601     sol1->ecrire(valeur2,compte,2,2);
602     compte++;
603     }
604     }
605 sattarpa 646 // affichage sur GMSH ******************************
606     if (gmsh_affiche>1)
607     {
608     cout<<"affiche: calcul_courbure_cadmgmai "<<endl;
609     std::ostringstream oss;
610     oss << "calcul_courbure_cadmgmai_curvdif_coef:" <<curvdif_coef ;
611     std::string namefic = oss.str();
612     char *cstr = new char[namefic.length() + 1];
613     strcpy(cstr, namefic.c_str());
614     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
615     gest.ajouter_mg_solution(sol1);
616     sol1->change_legende(0,"Kmin");
617     sol1->change_legende(1,"Kmax");
618     sol1->change_legende(2,"K_gau");
619     LISTE_MG_TRIANGLE::iterator ittri;
620     int compte=0;
621     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
622     {
623     MG_NOEUD *no1=tri->get_noeud1();
624     MG_NOEUD *no2=tri->get_noeud2();
625     MG_NOEUD *no3=tri->get_noeud3();
626     char mess[500];
627     sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
628     double valeur0=liste_kmin[mess];
629     double valeur1=liste_kmax[mess];
630     double valeur2=liste_gau[mess];
631     sol1->ecrire(valeur0,compte,0,0);
632     sol1->ecrire(valeur1,compte,1,0);
633     sol1->ecrire(valeur2,compte,2,0);
634     sprintf(mess,"%lu_%lu",no2->get_id(),tri->get_lien_topologie()->get_id());
635     valeur0=liste_kmin[mess];
636     valeur1=liste_kmax[mess];
637     valeur2=liste_gau[mess];
638     sol1->ecrire(valeur0,compte,0,1);
639     sol1->ecrire(valeur1,compte,1,1);
640     sol1->ecrire(valeur2,compte,2,1);
641     sprintf(mess,"%lu_%lu",no3->get_id(),tri->get_lien_topologie()->get_id());
642     valeur0=liste_kmin[mess];
643     valeur1=liste_kmax[mess];
644     valeur2=liste_gau[mess];
645     sol1->ecrire(valeur0,compte,0,2);
646     sol1->ecrire(valeur1,compte,1,2);
647     sol1->ecrire(valeur2,compte,2,2);
648     compte++;
649     }
650     }
651 sattarpa 609 liste_kmin.clear();liste_kmax.clear();liste_gau.clear();liste_courburemax1.clear();liste_courburemax2.clear();liste_courburemax3.clear();liste_courburemax4.clear();liste_courburemax5.clear();
652 sattarpa 650
653     cout<<"calcul courbure defcadmgmai"<<endl;
654 sattarpa 609 // calcul courbure defcadmgmai ********************************
655 sattarpa 554
656 sattarpa 609 for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
657     {
658     double kmin1=0.;
659     double kmax1=0.;
660     double kmin2=0.;
661     double kmax2=0.;
662     MG_MAILLAGE_OUTILS *cou;
663     cou->calcul_courbure_face_arete_sommet(defcadmgnd,kmin1,kmax1,kmin2,kmax2,liste_kmin,liste_kmax,liste_gau,liste_courburemax1,liste_courburemax2,liste_courburemax3,liste_courburemax4,liste_courburemax5,opensurafece);
664     }
665     for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
666     {
667     unsigned long idface=defcadmgnd->get_lien_topologie()->get_id();
668     unsigned long idnoeud=defcadmgnd->get_id();
669     char mess[500];
670     sprintf(mess,"%lu_%lu",idnoeud,idface);
671     for(std::map<std::string,double>::iterator itkmin=liste_kmin.begin();itkmin!=liste_kmin.end();itkmin++)
672     {
673     if(mess==(*itkmin).first)
674 sattarpa 650 {
675 sattarpa 609 lstkmin_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,defcadmgnd) );
676 sattarpa 650 lstkmin_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmin).second));
677     }
678 sattarpa 609 }
679     for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
680     {
681     if(mess==(*itkmax).first)
682 sattarpa 650 {
683 sattarpa 609 lstkmax_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,defcadmgnd) );
684 sattarpa 650 lstkmax_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmax).second));
685     }
686 sattarpa 609 }
687 sattarpa 646 for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
688     {
689     if(mess==(*itkgau).first)
690 sattarpa 650 {
691 sattarpa 646 lstkgau_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,defcadmgnd) );
692 sattarpa 650 lstkgau_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkgau).second));
693     }
694 sattarpa 646 }
695 sattarpa 609 }
696 sattarpa 650 //calc reletive curvature defCAD
697     map<MG_NOEUD*,double >::iterator itlstkmin_defcad;
698     map<MG_NOEUD*,double >::iterator itlstkmax_defcad;
699     map<MG_NOEUD*,double >::iterator itlstkgau_defcad;
700     for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
701     {
702     itlstkmin_defcad=lstkmin_defcad.find(defcadmgnd);
703     double relet_kmindefcad=(*itlstkmin_defcad).second;
704     itlstkmax_defcad=lstkmax_defcad.find(defcadmgnd);
705     double relet_kmaxdefcad=(*itlstkmax_defcad).second;
706     itlstkgau_defcad=lstkgau_cad.find(defcadmgnd);
707     double relet_kgaudefcad=(*itlstkgau_defcad).second;
708     TPL_MAP_ENTITE<MG_NOEUD*> lstnei_reldefcad;
709     int reldefcad_neindnb=lstnei_reldefcad.get_nb();
710     double relet_search_rad1=relet_search_rad;
711     while(reldefcad_neindnb<1)
712     {
713     octreecad.rechercher(defcadmgnd->get_x(),defcadmgnd->get_y(),defcadmgnd->get_z(),relet_search_rad1,lstnei_reldefcad);
714     reldefcad_neindnb=lstnei_reldefcad.get_nb();
715     relet_search_rad1=relet_search_rad1+0.1*relet_search_rad;
716     }
717     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
718     for(MG_NOEUD* defcadndnei=lstnei_reldefcad.get_premier(itndnei);defcadndnei!=NULL;defcadndnei=lstnei_reldefcad.get_suivant(itndnei))
719     {
720     itlstkmin_cad=lstkmin_cad.find(defcadndnei);
721     relet_kmindefcad=relet_kmindefcad+(*itlstkmin_defcad).second;
722     itlstkmax_cad=lstkmax_cad.find(defcadndnei);
723     relet_kmaxdefcad=relet_kmaxdefcad+(*itlstkmax_defcad).second;
724     itlstkgau_cad=lstkgau_cad.find(defcadndnei);
725     relet_kgaudefcad=relet_kgaudefcad+(*itlstkgau_defcad).second;
726     }
727     relet_kmindefcad=relet_kmindefcad/lstnei_reldefcad.get_nb();
728     relet_kmaxdefcad=relet_kmaxdefcad/lstnei_reldefcad.get_nb();
729     relet_kgaudefcad=relet_kgaudefcad/lstnei_reldefcad.get_nb();
730     lstkmin_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmindefcad,defcadmgnd));
731     lstkmax_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmaxdefcad,defcadmgnd));
732     lstkgau_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kgaudefcad,defcadmgnd));
733     }
734     // affichage reletive defcad sur GMSH ******************************
735     if (gmsh_affiche>1)
736     {
737     cout<<"affiche: calcul_courbure_cadmgmai_relet "<<endl;
738     std::ostringstream oss;
739     oss << "calcul_courbure_cadmgmai_curvdif_relet_coef:" <<curvdif_coef ;
740     std::string namefic = oss.str();
741     char *cstr = new char[namefic.length() + 1];
742     strcpy(cstr, namefic.c_str());
743     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
744     gest.ajouter_mg_solution(sol1);
745     sol1->change_legende(0,"Kmin_relet");
746     sol1->change_legende(1,"Kmax_relet");
747     sol1->change_legende(2,"K_gau_relet");
748     LISTE_MG_TRIANGLE::iterator ittri;
749     int compte=0;
750     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
751     {
752     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
753     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
754     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
755     MG_NOEUD *no1=tri->get_noeud1();
756     MG_NOEUD *no2=tri->get_noeud2();
757     MG_NOEUD *no3=tri->get_noeud3();
758     double valeur0;
759     double valeur1;
760     double valeur2;
761     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
762     {
763     if (no1==(*itlstkmin_relt_cadmgmai).second)
764     valeur0=(*itlstkmin_relt_cadmgmai).first;
765     }
766     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
767     {
768     if (no1==(*itlstkmax_relt_cadmgmai).second)
769     valeur1=(*itlstkmax_relt_cadmgmai).first;
770     }
771     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
772     {
773     if (no1==(*itlstkgau_relt_cadmgmai).second)
774     valeur2=(*itlstkgau_relt_cadmgmai).first;
775     }
776     sol1->ecrire(valeur0,compte,0,0);
777     sol1->ecrire(valeur1,compte,1,0);
778     sol1->ecrire(valeur2,compte,2,0);
779     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
780     {
781     if (no2==(*itlstkmin_relt_cadmgmai).second)
782     valeur0=(*itlstkmin_relt_cadmgmai).first;
783     }
784     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
785     {
786     if (no2==(*itlstkmax_relt_cadmgmai).second)
787     valeur1=(*itlstkmax_relt_cadmgmai).first;
788     }
789     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
790     {
791     if (no2==(*itlstkgau_relt_cadmgmai).second)
792     valeur2=(*itlstkgau_relt_cadmgmai).first;
793     }
794     sol1->ecrire(valeur0,compte,0,1);
795     sol1->ecrire(valeur1,compte,1,1);
796     sol1->ecrire(valeur2,compte,2,1);
797     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
798     {
799     if (no3==(*itlstkmin_relt_cadmgmai).second)
800     valeur0=(*itlstkmin_relt_cadmgmai).first;
801     }
802     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
803     {
804     if (no3==(*itlstkmax_relt_cadmgmai).second)
805     valeur1=(*itlstkmax_relt_cadmgmai).first;
806     }
807     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
808     {
809     if (no3==(*itlstkgau_relt_cadmgmai).second)
810     valeur2=(*itlstkgau_relt_cadmgmai).first;
811     }
812     sol1->ecrire(valeur0,compte,0,2);
813     sol1->ecrire(valeur1,compte,1,2);
814     sol1->ecrire(valeur2,compte,2,2);
815     compte++;
816     }
817     }
818 sattarpa 646 // affichage sur GMSH ******************************
819     if (gmsh_affiche>1)
820     {
821     cout<<"affiche: calcul_courbure_defcadmgmai "<<endl;
822     std::ostringstream oss;
823     oss << "calcul_courbure_defcadmgmai_curvdif_coef:" <<curvdif_coef ;
824     std::string namefic = oss.str();
825     char *cstr = new char[namefic.length() + 1];
826     strcpy(cstr, namefic.c_str());
827     MG_SOLUTION* sol2=new MG_SOLUTION(defcadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
828     gest.ajouter_mg_solution(sol2);
829     sol2->change_legende(0,"Kmin");
830     sol2->change_legende(1,"Kmax");
831     sol2->change_legende(2,"K_gau");
832     LISTE_MG_TRIANGLE::iterator ittri;
833     int compte=0;
834     for (MG_TRIANGLE* tri=defcadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=defcadmgmai->get_suivant_triangle(ittri))
835     {
836     MG_NOEUD *no1=tri->get_noeud1();
837     MG_NOEUD *no2=tri->get_noeud2();
838     MG_NOEUD *no3=tri->get_noeud3();
839     char mess[500];
840     sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
841     double valeur0=liste_kmin[mess];
842     double valeur1=liste_kmax[mess];
843     double valeur2=liste_gau[mess];
844     sol2->ecrire(valeur0,compte,0,0);
845     sol2->ecrire(valeur1,compte,1,0);
846     sol2->ecrire(valeur2,compte,2,0);
847     sprintf(mess,"%lu_%lu",no2->get_id(),tri->get_lien_topologie()->get_id());
848     valeur0=liste_kmin[mess];
849     valeur1=liste_kmax[mess];
850     valeur2=liste_gau[mess];
851     sol2->ecrire(valeur0,compte,0,1);
852     sol2->ecrire(valeur1,compte,1,1);
853     sol2->ecrire(valeur2,compte,2,1);
854     sprintf(mess,"%lu_%lu",no3->get_id(),tri->get_lien_topologie()->get_id());
855     valeur0=liste_kmin[mess];
856     valeur1=liste_kmax[mess];
857     valeur2=liste_gau[mess];
858     sol2->ecrire(valeur0,compte,0,2);
859     sol2->ecrire(valeur1,compte,1,2);
860     sol2->ecrire(valeur2,compte,2,2);
861     compte++;
862     }
863     char chaine[1000];
864     char* p;
865     p=strrchr(magicfilename,'.');
866     strncpy(chaine,magicfilename,p-magicfilename);
867     std::string namfic=chaine;
868     namfic=namfic + "courbcal.magic";
869     gest.enregistrer(namfic.c_str());
870     }
871     liste_kmin.clear();liste_kmax.clear();liste_gau.clear();liste_courburemax1.clear();liste_courburemax2.clear();liste_courburemax3.clear();liste_courburemax4.clear();liste_courburemax5.clear();
872    
873    
874     ///make correspondence between cadmgmai and defcadmgmai nodes
875 sattarpa 609 map<unsigned long, unsigned long>::iterator itcoresspondndidlst;
876     for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
877     {
878     itcoresspondndidlst=coresspondndidlst.find(defcadndcurv->get_id());
879     defcadndcurv->change_nouveau_numero((*itcoresspondndidlst).second);
880     }
881    
882     cout<<"solution calculation"<<endl;
883     multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_diff;
884     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_diff;
885 sattarpa 646 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_diff;
886 sattarpa 609 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_cadmgmai;
887     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_cadmgmai;
888 sattarpa 646 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_cadmgmai;
889 sattarpa 609 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_defcadmgmai;
890     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_defcadmgmai;
891 sattarpa 646 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_defcadmgmai;
892 sattarpa 609 cout<<"lstkmin_cadmgmai.size()= "<<lstkmin_cadmgmai.size()<<endl;
893     cout<<"lstkmin_defcadmgmai.size()= "<<lstkmin_defcadmgmai.size()<<endl;
894 sattarpa 646 cout<<"kmin diff. calc. started"<<endl;
895 sattarpa 650
896 sattarpa 609 for(itlstkmin_cadmgmai=lstkmin_cadmgmai.begin();itlstkmin_cadmgmai!=lstkmin_cadmgmai.end();itlstkmin_cadmgmai++)
897 sattarpa 554 {
898 sattarpa 609 MG_NOEUD* cadndcurv=(*itlstkmin_cadmgmai).second;
899     double kmindcad;
900     double kmindcad_def;
901     kmindcad=(*itlstkmin_cadmgmai).first;
902     for(itlstkmin_defcadmgmai=lstkmin_defcadmgmai.begin();itlstkmin_defcadmgmai!=lstkmin_defcadmgmai.end();itlstkmin_defcadmgmai++)
903     {
904     MG_NOEUD* defcadndcurv=(*itlstkmin_defcadmgmai).second;
905     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
906     kmindcad_def=(*itlstkmin_defcadmgmai).first;
907     }
908    
909     double kmin_diff=kmindcad_def-kmindcad;
910     lstkmin_diff.insert(pair<double,MG_NOEUD*>(kmin_diff,cadndcurv) );
911     }
912 sattarpa 646 cout<<"kmax diff. calc. started"<<endl;
913 sattarpa 609 for(itlstkmax_cadmgmai=lstkmax_cadmgmai.begin();itlstkmax_cadmgmai!=lstkmax_cadmgmai.end();itlstkmax_cadmgmai++)
914     {
915     MG_NOEUD* cadndcurv=(*itlstkmax_cadmgmai).second;
916     double kmaxndcad;
917     double kmaxndcad_def;
918     kmaxndcad=(*itlstkmax_cadmgmai).first;
919     for(itlstkmax_defcadmgmai=lstkmax_defcadmgmai.begin();itlstkmax_defcadmgmai!=lstkmax_defcadmgmai.end();itlstkmax_defcadmgmai++)
920     {
921     MG_NOEUD* defcadndcurv=(*itlstkmax_defcadmgmai).second;
922     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
923     {
924     kmaxndcad_def=(*itlstkmax_defcadmgmai).first;
925     }
926     }
927     double kmax_diff=kmaxndcad_def-kmaxndcad;
928     lstkmax_diff.insert(pair<double,MG_NOEUD*>(kmax_diff,cadndcurv) );
929     }
930 sattarpa 646 cout<<"kgau diff. calc. started"<<endl;
931     for(itlstkgau_cadmgmai=lstkgau_cadmgmai.begin();itlstkgau_cadmgmai!=lstkgau_cadmgmai.end();itlstkgau_cadmgmai++)
932     {
933     MG_NOEUD* cadndcurv=(*itlstkgau_cadmgmai).second;
934     double kgaundcad;
935     double kgaundcad_def;
936     kgaundcad=(*itlstkgau_cadmgmai).first;
937     for(itlstkgau_defcadmgmai=lstkgau_defcadmgmai.begin();itlstkgau_defcadmgmai!=lstkgau_defcadmgmai.end();itlstkgau_defcadmgmai++)
938     {
939     MG_NOEUD* defcadndcurv=(*itlstkgau_defcadmgmai).second;
940     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
941     {
942     kgaundcad_def=(*itlstkgau_defcadmgmai).first;
943     }
944     }
945     double kgau_diff=kgaundcad_def-kgaundcad;
946     lstkgau_diff.insert(pair<double,MG_NOEUD*>(kgau_diff,cadndcurv) );
947     }
948 sattarpa 650 // relative diff
949     multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_diff;
950     multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_diff;
951     multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_diff;
952     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
953     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
954     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
955     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_defcadmgmai;
956     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_defcadmgmai;
957     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_defcadmgmai;
958 sattarpa 658 cout<<"kmin _relt diff. calc. started"<<endl;
959     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
960 sattarpa 650 {
961     MG_NOEUD* cadndcurv=(*itlstkmin_relt_cadmgmai).second;
962     double kmindcad;
963     double kmindcad_relt_def;
964     kmindcad=(*itlstkmin_relt_cadmgmai).first;
965     for(itlstkmin_relt_defcadmgmai=lstkmin_relt_defcadmgmai.begin();itlstkmin_relt_defcadmgmai!=lstkmin_relt_defcadmgmai.end();itlstkmin_relt_defcadmgmai++)
966     {
967     MG_NOEUD* defcadndcurv=(*itlstkmin_relt_defcadmgmai).second;
968     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
969     kmindcad_relt_def=(*itlstkmin_relt_defcadmgmai).first;
970     }
971    
972     double kmin_relt_diff=kmindcad_relt_def-kmindcad;
973     lstkmin_relt_diff.insert(pair<double,MG_NOEUD*>(kmin_relt_diff,cadndcurv) );
974     }
975     cout<<"kmax _relt diff. calc. started"<<endl;
976     for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
977     {
978     MG_NOEUD* cadndcurv=(*itlstkmax_relt_cadmgmai).second;
979     double kmaxndcad;
980     double kmaxndcad_relt_def;
981     kmaxndcad=(*itlstkmax_relt_cadmgmai).first;
982     for(itlstkmax_relt_defcadmgmai=lstkmax_relt_defcadmgmai.begin();itlstkmax_relt_defcadmgmai!=lstkmax_relt_defcadmgmai.end();itlstkmax_relt_defcadmgmai++)
983     {
984     MG_NOEUD* defcadndcurv=(*itlstkmax_relt_defcadmgmai).second;
985     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
986     {
987     kmaxndcad_relt_def=(*itlstkmax_relt_defcadmgmai).first;
988     }
989     }
990     double kmax_relt_diff=kmaxndcad_relt_def-kmaxndcad;
991     lstkmax_relt_diff.insert(pair<double,MG_NOEUD*>(kmax_relt_diff,cadndcurv) );
992     }
993     cout<<"kgau_relt diff. calc. started"<<endl;
994     for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
995     {
996     MG_NOEUD* cadndcurv=(*itlstkgau_relt_cadmgmai).second;
997     double kgaundcad;
998     double kgaundcad_relt_def;
999     kgaundcad=(*itlstkgau_relt_cadmgmai).first;
1000     for(itlstkgau_relt_defcadmgmai=lstkgau_relt_defcadmgmai.begin();itlstkgau_relt_defcadmgmai!=lstkgau_relt_defcadmgmai.end();itlstkgau_relt_defcadmgmai++)
1001     {
1002     MG_NOEUD* defcadndcurv=(*itlstkgau_relt_defcadmgmai).second;
1003     if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
1004     {
1005     kgaundcad_relt_def=(*itlstkgau_relt_defcadmgmai).first;
1006     }
1007     }
1008     double kgau_relt_diff=kgaundcad_relt_def-kgaundcad;
1009     lstkgau_relt_diff.insert(pair<double,MG_NOEUD*>(kgau_relt_diff,cadndcurv) );
1010     }
1011     // affichage reletive diff sur GMSH ******************************
1012     if (gmsh_affiche>0)
1013     {
1014     cout<<"affiche: calcul_relet_difference_courbure "<<endl;
1015     std::ostringstream oss;
1016     oss << "KMaxMinGaudiff_curvdifcoef_relet:" <<curvdif_coef ;
1017     std::string namefic = oss.str();
1018     char *cstr = new char[namefic.length() + 1];
1019     strcpy(cstr, namefic.c_str());
1020     MG_SOLUTION* sol4=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"Kdiff_relt_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
1021     gest.ajouter_mg_solution(sol4);
1022     sol4->change_legende(0,"Kmin_diff");
1023     sol4->change_legende(1,"Kmax_diff");
1024     sol4->change_legende(2,"Kgau_diff");
1025     LISTE_MG_TRIANGLE::iterator ittri;
1026     int compte=0;
1027     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
1028     {
1029     MG_NOEUD *no1=tri->get_noeud1();
1030     MG_NOEUD *no2=tri->get_noeud2();
1031     MG_NOEUD *no3=tri->get_noeud3();
1032     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diff;
1033     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diff;
1034     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_diff;
1035     double valeurKmin1, valeurKmin2, valeurKmin3;
1036     double valeurKmax1, valeurKmax2, valeurKmax3;
1037     double valeurKgau1, valeurKgau2, valeurKgau3;
1038     for (itlstkmin_relt_diff=lstkmin_relt_diff.begin();itlstkmin_relt_diff!=lstkmin_relt_diff.end();itlstkmin_relt_diff++)
1039     {
1040     if (no1==(*itlstkmin_relt_diff).second)
1041     valeurKmin1=(*itlstkmin_relt_diff).first;
1042     if (no2==(*itlstkmin_relt_diff).second)
1043     valeurKmin2=(*itlstkmin_relt_diff).first;
1044     if (no3==(*itlstkmin_relt_diff).second)
1045     valeurKmin3=(*itlstkmin_relt_diff).first;
1046     }
1047     for (itlstkmax_relt_diff=lstkmax_relt_diff.begin();itlstkmax_relt_diff!=lstkmax_relt_diff.end();itlstkmax_relt_diff++)
1048     {
1049     if (no1==(*itlstkmax_relt_diff).second)
1050     valeurKmax1=(*itlstkmax_relt_diff).first;
1051     if (no2==(*itlstkmax_relt_diff).second)
1052     valeurKmax2=(*itlstkmax_relt_diff).first;
1053     if (no3==(*itlstkmax_relt_diff).second)
1054     valeurKmax3=(*itlstkmax_relt_diff).first;
1055     }
1056     for (itlstkgau_relt_diff=lstkgau_relt_diff.begin();itlstkgau_relt_diff!=lstkgau_relt_diff.end();itlstkgau_relt_diff++)
1057     {
1058     if (no1==(*itlstkgau_relt_diff).second)
1059     valeurKgau1=(*itlstkgau_relt_diff).first;
1060     if (no2==(*itlstkgau_relt_diff).second)
1061     valeurKgau2=(*itlstkgau_relt_diff).first;
1062     if (no3==(*itlstkgau_relt_diff).second)
1063     valeurKgau3=(*itlstkgau_relt_diff).first;
1064     }
1065     sol4->ecrire(valeurKmin1,compte,0,0);
1066     sol4->ecrire(valeurKmax1,compte,1,0);
1067     sol4->ecrire(valeurKgau1,compte,2,0);
1068    
1069     sol4->ecrire(valeurKmin2,compte,0,1);
1070     sol4->ecrire(valeurKmax2,compte,1,1);
1071     sol4->ecrire(valeurKgau2,compte,2,1);
1072    
1073     sol4->ecrire(valeurKmin3,compte,0,2);
1074     sol4->ecrire(valeurKmax3,compte,1,2);
1075     sol4->ecrire(valeurKgau3,compte,2,2);
1076     compte++;
1077     }
1078     }
1079    
1080 sattarpa 646 // affichage sur GMSH ******************************
1081 sattarpa 658 if (gmsh_affiche>1)
1082 sattarpa 646 {
1083     cout<<"affiche: calcul_difference_courbure "<<endl;
1084     std::ostringstream oss;
1085     oss << "KMaxMinGaudiff_curvdifcoef:" <<curvdif_coef ;
1086     std::string namefic = oss.str();
1087     char *cstr = new char[namefic.length() + 1];
1088     strcpy(cstr, namefic.c_str());
1089     MG_SOLUTION* sol3=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"Kdiff_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
1090     gest.ajouter_mg_solution(sol3);
1091     sol3->change_legende(0,"Kmin_diff");
1092     sol3->change_legende(1,"Kmax_diff");
1093     sol3->change_legende(2,"Kgau_diff");
1094     LISTE_MG_TRIANGLE::iterator ittri;
1095     int compte=0;
1096     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
1097     {
1098     MG_NOEUD *no1=tri->get_noeud1();
1099     MG_NOEUD *no2=tri->get_noeud2();
1100     MG_NOEUD *no3=tri->get_noeud3();
1101     //char mess[500];
1102     //sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
1103     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_diff;
1104     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_diff;
1105     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_diff;
1106     double valeurKmin1, valeurKmin2, valeurKmin3;
1107     double valeurKmax1, valeurKmax2, valeurKmax3;
1108     double valeurKgau1, valeurKgau2, valeurKgau3;
1109     for (itlstkmin_diff=lstkmin_diff.begin();itlstkmin_diff!=lstkmin_diff.end();itlstkmin_diff++)
1110     {
1111     if (no1==(*itlstkmin_diff).second)
1112     valeurKmin1=(*itlstkmin_diff).first;
1113     if (no2==(*itlstkmin_diff).second)
1114     valeurKmin2=(*itlstkmin_diff).first;
1115     if (no3==(*itlstkmin_diff).second)
1116     valeurKmin3=(*itlstkmin_diff).first;
1117     }
1118     for (itlstkmax_diff=lstkmax_diff.begin();itlstkmax_diff!=lstkmax_diff.end();itlstkmax_diff++)
1119     {
1120     if (no1==(*itlstkmax_diff).second)
1121     valeurKmax1=(*itlstkmax_diff).first;
1122     if (no2==(*itlstkmax_diff).second)
1123     valeurKmax2=(*itlstkmax_diff).first;
1124     if (no3==(*itlstkmax_diff).second)
1125     valeurKmax3=(*itlstkmax_diff).first;
1126     }
1127     for (itlstkgau_diff=lstkgau_diff.begin();itlstkgau_diff!=lstkgau_diff.end();itlstkgau_diff++)
1128     {
1129     if (no1==(*itlstkgau_diff).second)
1130     valeurKgau1=(*itlstkgau_diff).first;
1131     if (no2==(*itlstkgau_diff).second)
1132     valeurKgau2=(*itlstkgau_diff).first;
1133     if (no3==(*itlstkgau_diff).second)
1134     valeurKgau3=(*itlstkgau_diff).first;
1135     }
1136     sol3->ecrire(valeurKmin1,compte,0,0);
1137     sol3->ecrire(valeurKmax1,compte,1,0);
1138     sol3->ecrire(valeurKgau1,compte,2,0);
1139 sattarpa 609
1140 sattarpa 646 sol3->ecrire(valeurKmin2,compte,0,1);
1141     sol3->ecrire(valeurKmax2,compte,1,1);
1142     sol3->ecrire(valeurKgau2,compte,2,1);
1143    
1144     sol3->ecrire(valeurKmin3,compte,0,2);
1145     sol3->ecrire(valeurKmax3,compte,1,2);
1146     sol3->ecrire(valeurKgau3,compte,2,2);
1147     compte++;
1148     }
1149     char chaine[1000];
1150     char* p;
1151     p=strrchr(magicfilename,'.');
1152     strncpy(chaine,magicfilename,p-magicfilename);
1153     std::string namfic=chaine;
1154     namfic=namfic + "courbcal.magic";
1155     gest.enregistrer(namfic.c_str());
1156     }
1157 sattarpa 609 cout<<"read insert points in lstpinsert"<<endl;
1158     //read insert points in lstpinsert
1159     std::vector<double*> lstpinsert;
1160     double prop=0.001;
1161     FILE *in=fopen(gnifinspointfile,"rt");
1162     if (in==NULL) cout<<"file is not available"<<endl;
1163     while (!feof(in))
1164     {
1165     char chaine[500];
1166 sattarpa 655 char* res=fgets(chaine,500,in);
1167     if(res!=NULL)
1168     {
1169 sattarpa 609 double x,y,z;
1170     double q1,q2,q3;
1171     int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
1172     q1=prop*q1;
1173     q2=prop*q2;
1174     q3=prop*q3;
1175     if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
1176     else if (nb==6)
1177     {
1178     double* pcoordbc=new double[6];
1179     pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
1180     lstpinsert.push_back(pcoordbc);
1181     }
1182 sattarpa 655 }
1183 sattarpa 609 }
1184     fclose(in);
1185 sattarpa 650 if (relet_curvature<1)
1186     {
1187 sattarpa 609 //criterion
1188 sattarpa 658 multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmin_diff=lstkmin_diff.rbegin();
1189     multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmax_diff=lstkmax_diff.rbegin();
1190     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_diff=lstkmin_diff.begin();
1191     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_diff=lstkmax_diff.begin();
1192     double kmin_diff_min=(*ritlstkmin_diff).first;
1193     double kmin_diff_max=(*itlstkmin_diff).first;
1194     double kmax_diff_min=(*ritlstkmax_diff).first;
1195     double kmax_diff_max=(*itlstkmax_diff).first;
1196     cout<<"kmin_diff_min= "<<kmin_diff_min<<" kmin_diff_max= "<<kmin_diff_max<<endl;
1197     cout<<"kmax_diff_min= "<<kmax_diff_min<<" kmax_diff_max= "<<kmax_diff_max<<endl;
1198    
1199     double crikmin_min=curvdif_coef*kmin_diff_min;
1200     double crikmin_max=curvdif_coef*kmin_diff_max;
1201     double crikmax_min=curvdif_coef*kmax_diff_min;
1202     double crikmax_max=curvdif_coef*kmax_diff_max;
1203     cout<<"crikmin_min= "<<crikmin_min<<" crikmin_max= "<<crikmin_max<<" crikmax_min= "<<crikmax_min<<" crikmax_max= "<<crikmax_max<<endl;
1204 sattarpa 609 cout<<"lstkmin_diff.size()= "<<lstkmin_diff.size()<<endl;
1205     cout<<"lstkmax_diff.size()= "<<lstkmax_diff.size()<<endl;
1206 sattarpa 650
1207     cout<<"for Kmin"<<endl;
1208 sattarpa 609 TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
1209     for (itlstkmin_diff=lstkmin_diff.begin();itlstkmin_diff!=lstkmin_diff.end();itlstkmin_diff++)
1210     {
1211     {
1212     if((crikmin_max<(*itlstkmin_diff).first) || (crikmin_min>(*itlstkmin_diff).first))
1213     {
1214     MG_NOEUD* piefect=(*itlstkmin_diff).second;
1215     int neindnb=lstneipdefect.get_nb();
1216     double search_radius1=search_radius;
1217     while(neindnb==0)
1218     {
1219 sattarpa 650 octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1220 sattarpa 609 neindnb=lstneipdefect.get_nb();
1221     search_radius1=search_radius1+0.1*search_radius;
1222     }
1223     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
1224     for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
1225     {
1226     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1227     {
1228     double* xyzpins;
1229     xyzpins=(*itlstpinsert);
1230 sattarpa 650 if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
1231     lstpinsert.erase(itlstpinsert);
1232 sattarpa 609 }
1233     }
1234     lstneipdefect.vide();
1235     }
1236     }
1237     }
1238     cout<<"for Kmax"<<endl;
1239     for (itlstkmax_diff=lstkmax_diff.begin();itlstkmax_diff!=lstkmax_diff.end();itlstkmax_diff++)
1240     {
1241     if((crikmax_max<(*itlstkmax_diff).first) || (crikmax_min>(*itlstkmax_diff).first))
1242     {
1243     MG_NOEUD* piefect=(*itlstkmax_diff).second;
1244     TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
1245     int neindnb=lstneipdefect.get_nb();
1246     double search_radius1=search_radius;
1247     while(neindnb==0)
1248     {
1249 sattarpa 650 octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1250 sattarpa 609 neindnb=lstneipdefect.get_nb();
1251     search_radius1=search_radius1+0.1*search_radius;
1252     }
1253     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
1254     for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
1255     {
1256     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1257     {
1258     double* xyzpins;
1259     xyzpins=(*itlstpinsert);
1260     if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
1261     lstpinsert.erase(itlstpinsert);
1262     }
1263     }
1264    
1265     }
1266     }
1267 sattarpa 650 }
1268     if (relet_curvature>0)
1269     {
1270     //criterion
1271 sattarpa 658 multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmin_relt_diff=lstkmin_relt_diff.rbegin();
1272     multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmax_relt_diff=lstkmax_relt_diff.rbegin();
1273     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diff=lstkmin_relt_diff.begin();
1274     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diff=lstkmax_relt_diff.begin();
1275     double kmin_relt_diff_min=(*ritlstkmin_relt_diff).first;
1276     double kmin_relt_diff_max=(*itlstkmin_relt_diff).first;
1277     double kmax_relt_diff_min=(*ritlstkmax_relt_diff).first;
1278     double kmax_relt_diff_max=(*itlstkmax_relt_diff).first;
1279     cout<<"kmin_relt_diff_min= "<<kmin_relt_diff_min<<" kmin_relt_diff_max= "<<kmin_relt_diff_max<<endl;
1280     cout<<"kmax_relt_diff_min= "<<kmax_relt_diff_min<<" kmax_relt_diff_max= "<<kmax_relt_diff_max<<endl;
1281    
1282     double crikmin_min=curvdif_coef*kmin_relt_diff_min;
1283     double crikmin_max=curvdif_coef*kmin_relt_diff_max;
1284     double crikmax_min=curvdif_coef*kmax_relt_diff_min;
1285     double crikmax_max=curvdif_coef*kmax_relt_diff_max;
1286     cout<<"crikmin_min= "<<crikmin_min<<" crikmin_max= "<<crikmin_max<<" crikmax_min= "<<crikmax_min<<" crikmax_max= "<<crikmax_max<<endl;
1287 sattarpa 650 cout<<"lstkmin_diff.size()= "<<lstkmin_diff.size()<<endl;
1288     cout<<"lstkmax_diff.size()= "<<lstkmax_diff.size()<<endl;
1289    
1290 sattarpa 658 cout<<"for Kmin_relet_diff"<<endl;
1291     cout<<"lstkmin_relt_diff.size()= "<<lstkmin_relt_diff.size()<<endl;
1292 sattarpa 650 TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
1293     for (itlstkmin_relt_diff=lstkmin_relt_diff.begin();itlstkmin_relt_diff!=lstkmin_relt_diff.end();itlstkmin_relt_diff++)
1294     {
1295     {
1296     if((crikmin_max<(*itlstkmin_relt_diff).first) || (crikmin_min>(*itlstkmin_relt_diff).first))
1297     {
1298     MG_NOEUD* piefect=(*itlstkmin_relt_diff).second;
1299     int neindnb=lstneipdefect.get_nb();
1300     double search_radius1=search_radius;
1301     while(neindnb==0)
1302     {
1303     octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1304     neindnb=lstneipdefect.get_nb();
1305     search_radius1=search_radius1+0.1*search_radius;
1306     }
1307     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
1308     for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
1309     {
1310     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1311     {
1312     double* xyzpins;
1313     xyzpins=(*itlstpinsert);
1314 sattarpa 658
1315 sattarpa 650 if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
1316 sattarpa 658 {
1317 sattarpa 650 lstpinsert.erase(itlstpinsert);
1318 sattarpa 658 itlstpinsert=lstpinsert.begin();
1319     }
1320    
1321 sattarpa 650 }
1322     }
1323     lstneipdefect.vide();
1324     }
1325     }
1326     }
1327 sattarpa 658 cout<<"for Kmax_relet_diff"<<endl;
1328    
1329 sattarpa 650 for (itlstkmax_relt_diff=lstkmax_relt_diff.begin();itlstkmax_relt_diff!=lstkmax_relt_diff.end();itlstkmax_relt_diff++)
1330     {
1331     if((crikmax_max<(*itlstkmax_relt_diff).first) || (crikmax_min>(*itlstkmax_relt_diff).first))
1332     {
1333     MG_NOEUD* piefect=(*itlstkmax_relt_diff).second;
1334     TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
1335     int neindnb=lstneipdefect.get_nb();
1336     double search_radius1=search_radius;
1337     while(neindnb==0)
1338     {
1339     octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1340     neindnb=lstneipdefect.get_nb();
1341     search_radius1=search_radius1+0.1*search_radius;
1342     }
1343     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
1344     for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
1345     {
1346     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1347     {
1348     double* xyzpins;
1349     xyzpins=(*itlstpinsert);
1350     if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
1351 sattarpa 658 {
1352     lstpinsert.erase(itlstpinsert);
1353     itlstpinsert=lstpinsert.begin();
1354     }
1355 sattarpa 650 }
1356     }
1357    
1358     }
1359     }
1360     }
1361 sattarpa 609 FILE* in2=fopen(out_gnifinspointfile_removondefect,"wt");
1362 sattarpa 655 int countfuck=0;
1363 sattarpa 609 for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1364     {
1365     double* xyzpins1;
1366     xyzpins1=(*itlstpinsert);
1367 sattarpa 658 fprintf(in2,"%le %le %le %le %le %le\n",xyzpins1[0],xyzpins1[1],xyzpins1[2],xyzpins1[3]*1000.,xyzpins1[4]*1000.,xyzpins1[5]*1000.);
1368 sattarpa 609 }
1369 sattarpa 650 fclose(in2);
1370 sattarpa 609 }
1371     void CRIAQOPERATORS::vm_crit_remvsmp(int femmeshno,char* magicfilename,int femsolid,int femsubsolno,char* gnifinspointfile,char* out_gnifinspointfile_removondefect,double search_radius,double vm_coef)
1372     {
1373 sattarpa 554 MG_FILE gest(magicfilename);
1374 sattarpa 596 FEM_MAILLAGE* femmai;
1375     if (femmeshno==0) femmai=gest.get_fem_maillage(femmeshno); else femmai=gest.get_fem_maillageid(femmeshno);
1376     FEM_SOLUTION* femsol=gest.get_fem_solutionid(femsolid);
1377     femsol->active_solution(femsubsolno);
1378    
1379     std::map<double,FEM_NOEUD*,std::greater<double> > mpsol;
1380     LISTE_FEM_NOEUD::iterator itfemnds;
1381     for(FEM_NOEUD* nod=femmai->get_premier_noeud(itfemnds);nod!=NULL;nod=femmai->get_suivant_noeud(itfemnds))
1382     {
1383 sattarpa 609 nod->get_mg_element_maillage()->get_id();
1384 sattarpa 596 mpsol.insert(pair<double,FEM_NOEUD*> (nod->get_solution(),nod));
1385     //cout<<"nod->get_solution()= "<<nod->get_solution()<<endl;
1386     //if(nod->get_solution()<=1.e-20 && nod->get_solution()!=0.)
1387     //cout<<nod->get_id()<<" : "<<nod->get_solution()<<endl;
1388     }
1389     std::map<double,FEM_NOEUD*,std::greater<double> >::iterator itmpsol=mpsol.begin();
1390    
1391     double maxi=(*itmpsol).first;
1392     itmpsol=mpsol.end();
1393     double mini=(*itmpsol).first;
1394 sattarpa 609 double crit=vm_coef*(maxi+mini)/2.;
1395 sattarpa 596 //cout<<"critreion value= "<<crit<<endl;
1396     std::vector<double*> lstpinsert;
1397     double prop=0.001;
1398     FILE *in=fopen(gnifinspointfile,"rt");
1399     if (in==NULL) cout<<"file is not available"<<endl;
1400     while (!feof(in))
1401     {
1402     char chaine[500];
1403 sattarpa 655 char* res=fgets(chaine,500,in);
1404     if(res!=NULL)
1405     {
1406 sattarpa 596 double x,y,z;
1407     double q1,q2,q3;
1408     int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
1409     q1=prop*q1;
1410     q2=prop*q2;
1411     q3=prop*q3;
1412     if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
1413     else if (nb==6)
1414     {
1415     double* pcoordbc=new double[6];
1416     pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
1417     lstpinsert.push_back(pcoordbc);
1418     }
1419 sattarpa 655 }
1420 sattarpa 596 }
1421     fclose(in);
1422    
1423    
1424     //octree initialization
1425     double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
1426     TPL_MAP_ENTITE<FEM_NOEUD*> lstnoeud;
1427     LISTE_FEM_NOEUD::iterator it;
1428     for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
1429     {
1430     if (no->get_x()<xmin) xmin=no->get_x();
1431     if (no->get_y()<ymin) ymin=no->get_y();
1432     if (no->get_z()<zmin) zmin=no->get_z();
1433     if (no->get_x()>xmax) xmax=no->get_x();
1434     if (no->get_y()>ymax) ymax=no->get_y();
1435     if (no->get_z()>zmax) zmax=no->get_z();
1436     lstnoeud.ajouter(no);
1437     }
1438     OT_VECTEUR_3D vecmin(xmin,ymin,zmin);OT_VECTEUR_3D vecmmax(xmax,ymax,zmax);
1439     OT_VECTEUR_3D vec(vecmmax,vecmin);
1440    
1441     //double search_radius=0.001*(vec.get_longueur());
1442 sattarpa 609 //double search_radius=2.; //sasan this should be as a parameter
1443 sattarpa 596 OT_VECTEUR_3D min(xmin,ymin,zmin); OT_VECTEUR_3D max(xmax,ymax,zmax); OT_VECTEUR_3D average=(min+max)/2.; OT_VECTEUR_3D lengthvec(min,max);
1444     double length=sqrt(lengthvec*lengthvec); double bxr=1.1;
1445     xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
1446     xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);
1447    
1448     TPL_OCTREE<FEM_NOEUD*,FEM_NOEUD*> octree;
1449     octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
1450     for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
1451     octree.inserer(no);
1452    
1453     for (itmpsol=mpsol.begin();itmpsol!=mpsol.end();itmpsol++)
1454     {
1455     if(crit<(*itmpsol).first)
1456     {
1457     FEM_NOEUD* pdefect=(*itmpsol).second;
1458     TPL_MAP_ENTITE<FEM_NOEUD*> lstneipdefect;
1459     int neindnb=lstneipdefect.get_nb();
1460 sattarpa 609 double search_radius1=search_radius;
1461 sattarpa 596 while(neindnb==0)
1462     {
1463 sattarpa 609 octree.rechercher(pdefect->get_x(),pdefect->get_y(),pdefect->get_z(),search_radius1,lstneipdefect);
1464 sattarpa 596 neindnb=lstneipdefect.get_nb();
1465 sattarpa 609 search_radius1=search_radius1+0.1*search_radius;
1466 sattarpa 596 }
1467     TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR itfendnei;
1468     for(FEM_NOEUD* fendnei=lstneipdefect.get_premier(itfendnei);fendnei!=NULL;fendnei=lstneipdefect.get_suivant(itfendnei))
1469     {
1470     // std::map<double*,double*,std::less<double*> >::iterator itlstpinsert=lstpinsert.begin();//itlstpinsert!=lstpinsert.end();itlstpinsert++
1471     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1472     {
1473     double* xyzpins;
1474     xyzpins=(*itlstpinsert);
1475     if(fendnei->get_x()==xyzpins[0] && fendnei->get_y()==xyzpins[1] && fendnei->get_z()==xyzpins[2])
1476     {
1477     lstpinsert.erase(itlstpinsert);
1478 sattarpa 658 itlstpinsert=lstpinsert.begin();
1479 sattarpa 596 }
1480     }
1481     }
1482    
1483     }
1484    
1485     }
1486    
1487     FILE* in1=fopen(out_gnifinspointfile_removondefect,"wt");
1488     for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
1489     {
1490     double* xyzpins1;
1491     xyzpins1=(*itlstpinsert);
1492 sattarpa 658 fprintf(in1,"%le %le %le %le %le %le\n",xyzpins1[0],xyzpins1[1],xyzpins1[2],xyzpins1[3]*1000.,xyzpins1[4]*1000.,xyzpins1[5]*1000.);
1493 sattarpa 596 }
1494     fclose(in1);
1495     }
1496     void CRIAQOPERATORS::surfmaker_e1(char* outputfilename,int meshno,char* magicfilename)
1497     {
1498     MG_FILE gest(magicfilename);
1499 sattarpa 554 MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
1500     MG_MAILLAGE* mai;
1501     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1502     //MG_FILE gest((char*)"longflxprtwithpockets_3d_3mm.magic");
1503     //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
1504     //MG_MAILLAGE* mai=gest.get_mg_maillageid(5674824);
1505     MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
1506     MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1507     gestnew->ajouter_mg_geometrie(geonew);
1508     MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1509     gestnew->ajouter_mg_maillage(mainew);
1510    
1511     LISTE_MG_NOEUD::iterator itnod;
1512     for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
1513     {
1514     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
1515     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
1516     MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
1517     MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
1518 sattarpa 596
1519     if(facevert->get_id()==76
1520     || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==13 || aretver->get_id()==110||aretver->get_id()==98
1521     || aretver->get_id()==81 || aretver->get_id()==83 || aretver->get_id()==10|| aretver->get_id()==12|| aretver->get_id()==90|| aretver->get_id()==97
1522     /* facevert->get_id()==133
1523 sattarpa 554 || aretver->get_id()==355 || aretver->get_id()==362|| aretver->get_id()==374|| aretver->get_id()==298 || aretver->get_id()==291||aretver->get_id()==338
1524     || aretver->get_id()==326 || aretver->get_id()==333|| aretver->get_id()==319|| aretver->get_id()==312 || aretver->get_id()==305||aretver->get_id()==169
1525     || aretver->get_id()==162 || aretver->get_id()==141|| aretver->get_id()==174|| aretver->get_id()==148 || aretver->get_id()==155||aretver->get_id()==217
1526     || aretver->get_id()==212 || aretver->get_id()==191|| aretver->get_id()==184|| aretver->get_id()==64 || aretver->get_id()==280||aretver->get_id()==32
1527     || aretver->get_id()==387 || aretver->get_id()==198|| aretver->get_id()==205|| aretver->get_id()==269 || aretver->get_id()==262||aretver->get_id()==234
1528     || aretver->get_id()==227 || aretver->get_id()==255|| aretver->get_id()==248|| aretver->get_id()==241 || aretver->get_id()==380||aretver->get_id()==96
1529     || aretver->get_id()==369 || aretver->get_id()==344|| aretver->get_id()==128|| aretver->get_id()==274
1530     || sometver->get_id()==368|| sometver->get_id()==361|| sometver->get_id()==106|| sometver->get_id()==122|| sometver->get_id()==74|| sometver->get_id()==90
1531     || sometver->get_id()==233|| sometver->get_id()==240|| sometver->get_id()==247|| sometver->get_id()==254|| sometver->get_id()==261|| sometver->get_id()==268
1532     || sometver->get_id()==226|| sometver->get_id()==224|| sometver->get_id()==190|| sometver->get_id()==197|| sometver->get_id()==204|| sometver->get_id()==211
1533     || sometver->get_id()==183|| sometver->get_id()==181|| sometver->get_id()==154|| sometver->get_id()==147|| sometver->get_id()==140|| sometver->get_id()==161
1534     || sometver->get_id()==168|| sometver->get_id()==138|| sometver->get_id()==318|| sometver->get_id()==311|| sometver->get_id()==304|| sometver->get_id()==297
1535     || sometver->get_id()==325|| sometver->get_id()==332|| sometver->get_id()==290|| sometver->get_id()==288|| sometver->get_id()==42|| sometver->get_id()==58
1536 sattarpa 596 || sometver->get_id()==352|| sometver->get_id()==354|| sometver->get_id()==10|| sometver->get_id()==26*/
1537    
1538     )
1539 sattarpa 554 {
1540     double x=vertices->get_x();
1541     double y=vertices->get_y();
1542     double z=vertices->get_z();
1543     //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
1544     MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
1545     mainew->ajouter_mg_noeud(newno);
1546     vertices->change_nouveau_numero(newno->get_id()); // what does this do?
1547     }
1548     }
1549    
1550     LISTE_MG_SEGMENT::iterator itSeg;
1551     for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
1552     {
1553     //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
1554     //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
1555     //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
1556     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1557     MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
1558 sattarpa 596 if (faceseg->get_id()==76
1559     || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==13 || aretver->get_id()==110||aretver->get_id()==98
1560     //faceseg->get_id()==133
1561     //|| aretver->get_id()==355 || aretver->get_id()==362|| aretver->get_id()==374|| aretver->get_id()==298 || aretver->get_id()==291||aretver->get_id()==338
1562     //|| aretver->get_id()==326 || aretver->get_id()==333|| aretver->get_id()==319|| aretver->get_id()==312 || aretver->get_id()==305||aretver->get_id()==169
1563     //|| aretver->get_id()==162 || aretver->get_id()==141|| aretver->get_id()==174|| aretver->get_id()==148 || aretver->get_id()==155||aretver->get_id()==217
1564     //|| aretver->get_id()==212 || aretver->get_id()==191|| aretver->get_id()==184|| aretver->get_id()==64 || aretver->get_id()==280||aretver->get_id()==32
1565     //|| aretver->get_id()==387 || aretver->get_id()==198|| aretver->get_id()==205|| aretver->get_id()==269 || aretver->get_id()==262||aretver->get_id()==234
1566     // || aretver->get_id()==227 || aretver->get_id()==255|| aretver->get_id()==248|| aretver->get_id()==241 || aretver->get_id()==380||aretver->get_id()==96
1567     // || aretver->get_id()==369 || aretver->get_id()==344|| aretver->get_id()==128|| aretver->get_id()==274
1568    
1569     )
1570 sattarpa 554 {
1571     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1572     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1573     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
1574     //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
1575     mainew->ajouter_mg_segment(seg);
1576     }
1577     }
1578    
1579     LISTE_MG_TRIANGLE::iterator itTri;
1580     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
1581     {
1582     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
1583 sattarpa 596 if (facetri->get_id()==76)
1584 sattarpa 554 {
1585     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
1586     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
1587     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
1588     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
1589     }
1590     }
1591    
1592     //gestnew->enregistrer("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump.magic");
1593     gestnew->enregistrer(outputfilename);
1594     //////////////////
1595    
1596     ////////////////////////
1597    
1598    
1599     /*
1600     mainew->supprimer_mg_triangleid(106625);
1601     mainew->supprimer_mg_triangleid(106662);
1602     mainew->supprimer_mg_triangleid(106626);
1603     mainew->supprimer_mg_triangleid(106701);
1604     mainew->supprimer_mg_triangleid(106702);
1605     mainew->supprimer_mg_triangleid(106639);*/
1606     //cout<<"output mesh no.= 495557 "<<endl;
1607     //gest.enregistrer("2dsurfmesh_cad.magic"); //CAD model output
1608     //scanned part output
1609    
1610    
1611     }
1612     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1613 sattarpa 596 void CRIAQOPERATORS::surfmaker_e2(char* outputfilename,int meshno,char* magicfilename)
1614 sattarpa 554 {
1615     // surfmeshgenerator a partially scanned (deformed CAD model) for GNIF examples for simple exam with hole (badisplacenmentvector + smplwithpocket)
1616     MG_FILE gest(magicfilename);
1617     MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
1618     MG_MAILLAGE* mai;
1619     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1620     //MG_FILE gest((char*)"2d_smplwithole_23mm.magic");
1621     //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
1622     //MG_MAILLAGE* mai=gest.get_mg_maillageid(1197);
1623     MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
1624     MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1625     gestnew->ajouter_mg_geometrie(geonew);
1626     MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1627     gestnew->ajouter_mg_maillage(mainew);
1628    
1629     LISTE_MG_NOEUD::iterator itnod;
1630     for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
1631     {
1632     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
1633     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
1634     MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
1635     MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
1636     if(facevert->get_id()==76
1637     || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==98 || aretver->get_id()==13||aretver->get_id()==110
1638     || sometver->get_id()==10|| sometver->get_id()==12|| sometver->get_id()==81|| sometver->get_id()==83|| sometver->get_id()==90|| sometver->get_id()==97)
1639     {
1640     double x=vertices->get_x();
1641     double y=vertices->get_y();
1642     double z=vertices->get_z();
1643     //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
1644     MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
1645     mainew->ajouter_mg_noeud(newno);
1646     vertices->change_nouveau_numero(newno->get_id()); // what does this do?
1647     }
1648     }
1649    
1650     LISTE_MG_SEGMENT::iterator itSeg;
1651     for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
1652     {
1653     //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
1654     //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
1655     //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
1656     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1657     //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
1658     MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
1659     //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
1660     if (faceseg->get_id()==76 ||
1661     aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==98 || aretver->get_id()==13|| aretver->get_id()==110)
1662     {
1663     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1664     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1665     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
1666     //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
1667     mainew->ajouter_mg_segment(seg);
1668     }
1669     }
1670    
1671     LISTE_MG_TRIANGLE::iterator itTri;
1672     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
1673     {
1674     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
1675     if (facetri->get_id()==76)
1676     {
1677     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
1678     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
1679     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
1680     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
1681     }
1682     }
1683    
1684     //gestnew->enregistrer("cadsurfmsh_smplwithole_23mm.magic");
1685     gestnew->enregistrer(outputfilename);
1686     }
1687 sattarpa 596 void CRIAQOPERATORS::surfmaker_e2t(char* outputfilename,int meshno,char* magicfilename)
1688 sattarpa 554 {
1689     // surfmeshgenerator a partially scanned (deformed CAD model) for GNIF examples for simple exam with hole (badisplacenmentvector + smplwithpocket)
1690     MG_FILE gest(magicfilename);
1691     MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
1692     MG_MAILLAGE* mai;
1693     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1694     //MG_FILE gest((char*)"2d_smplwithole_23mm.magic");
1695     //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
1696     //MG_MAILLAGE* mai=gest.get_mg_maillageid(1197);
1697     MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
1698     MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1699     gestnew->ajouter_mg_geometrie(geonew);
1700     MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1701     gestnew->ajouter_mg_maillage(mainew);
1702    
1703     LISTE_MG_NOEUD::iterator itnod;
1704     for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
1705     {
1706     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
1707     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
1708     MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
1709     MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
1710     if(facevert->get_id()==37
1711     || aretver->get_id()==20 || aretver->get_id()==108|| aretver->get_id()==113|| aretver->get_id()==101 || aretver->get_id()==59||aretver->get_id()==64
1712     || aretver->get_id()==45 || aretver->get_id()==52|| aretver->get_id()==93|| aretver->get_id()==74 || aretver->get_id()==81||aretver->get_id()==88
1713     || sometver->get_id()==10|| sometver->get_id()==19|| sometver->get_id()==51|| sometver->get_id()==58|| sometver->get_id()==42|| sometver->get_id()==44
1714     || sometver->get_id()==87|| sometver->get_id()==73|| sometver->get_id()==80|| sometver->get_id()==71|| sometver->get_id()==100|| sometver->get_id()==107
1715     //for e2t
1716     //|| aretver->get_id()==20 || aretver->get_id()==141|| aretver->get_id()==129|| aretver->get_id()==136 || aretver->get_id()==66||aretver->get_id()==59
1717     //|| aretver->get_id()==73 || aretver->get_id()==78|| aretver->get_id()==45|| aretver->get_id()==52 || aretver->get_id()==95||aretver->get_id()==102
1718     //|| aretver->get_id()==88 || aretver->get_id()==121|| aretver->get_id()==116|| aretver->get_id()==109
1719     //|| sometver->get_id()==10|| sometver->get_id()==19|| sometver->get_id()==65|| sometver->get_id()==58|| sometver->get_id()==51|| sometver->get_id()==42
1720     //|| sometver->get_id()==44|| sometver->get_id()==72|| sometver->get_id()==115|| sometver->get_id()==87|| sometver->get_id()==108|| sometver->get_id()==101
1721     //|| sometver->get_id()==94|| sometver->get_id()==85|| sometver->get_id()==128|| sometver->get_id()==135
1722     )
1723     {
1724     double x=vertices->get_x();
1725     double y=vertices->get_y();
1726     double z=vertices->get_z();
1727     //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
1728     MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
1729     mainew->ajouter_mg_noeud(newno);
1730     vertices->change_nouveau_numero(newno->get_id()); // what does this do?
1731     }
1732     }
1733    
1734     LISTE_MG_SEGMENT::iterator itSeg;
1735     for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
1736     {
1737     //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
1738     //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
1739     //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
1740     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1741     //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
1742     MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
1743     //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
1744     if (faceseg->get_id()==37
1745     || aretver->get_id()==20 || aretver->get_id()==108|| aretver->get_id()==113|| aretver->get_id()==101 || aretver->get_id()==59||aretver->get_id()==64
1746     || aretver->get_id()==45 || aretver->get_id()==52|| aretver->get_id()==93|| aretver->get_id()==74 || aretver->get_id()==81||aretver->get_id()==88
1747     //|| aretver->get_id()==20 || aretver->get_id()==141|| aretver->get_id()==129|| aretver->get_id()==136 || aretver->get_id()==66||aretver->get_id()==59
1748     //|| aretver->get_id()==73 || aretver->get_id()==78|| aretver->get_id()==45|| aretver->get_id()==52 || aretver->get_id()==95||aretver->get_id()==102
1749     //|| aretver->get_id()==88 || aretver->get_id()==121|| aretver->get_id()==116|| aretver->get_id()==109
1750    
1751     )
1752     {
1753     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1754     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1755     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
1756     //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
1757     mainew->ajouter_mg_segment(seg);
1758     }
1759     }
1760    
1761     LISTE_MG_TRIANGLE::iterator itTri;
1762     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
1763     {
1764     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
1765     if (facetri->get_id()==37)
1766     {
1767     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
1768     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
1769     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
1770     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
1771     }
1772     }
1773    
1774     //gestnew->enregistrer("cadsurfmsh_smplwithole_23mm.magic");
1775     gestnew->enregistrer(outputfilename);
1776     }
1777     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1778 sattarpa 596 void CRIAQOPERATORS::surfmaker_e3(char* outputfilename,int meshno,char* magicfilename)
1779 sattarpa 554 {
1780     //surfmesh maker (smplwithpocket+smpl withpocket)
1781     MG_FILE gest(magicfilename);
1782     MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
1783     MG_MAILLAGE* mai;
1784     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1785     //MG_FILE gest((char*)"smplwithpocket_3d_13mm.magic");
1786     //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
1787     //MG_MAILLAGE* mai=gest.get_mg_maillageid(2982821);
1788     MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
1789     MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1790     gestnew->ajouter_mg_geometrie(geonew);
1791     MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1792     gestnew->ajouter_mg_maillage(mainew);
1793    
1794     LISTE_MG_NOEUD::iterator itnod;
1795     for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
1796     {
1797     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
1798     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
1799     MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
1800     MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
1801     if(facevert->get_id()==5
1802     || aretver->get_id()==32 || aretver->get_id()==61|| aretver->get_id()==42|| aretver->get_id()==56 || aretver->get_id()==49||aretver->get_id()==27
1803     || aretver->get_id()==85 || aretver->get_id()==90|| aretver->get_id()==78|| aretver->get_id()==71 || aretver->get_id()==20||aretver->get_id()==13
1804     || sometver->get_id()==10|| sometver->get_id()==26|| sometver->get_id()==41|| sometver->get_id()==55|| sometver->get_id()==39|| sometver->get_id()==48
1805     || sometver->get_id()==84|| sometver->get_id()==77|| sometver->get_id()==70|| sometver->get_id()==68|| sometver->get_id()==12|| sometver->get_id()==19)
1806     {
1807     double x=vertices->get_x();
1808     double y=vertices->get_y();
1809     double z=vertices->get_z();
1810     //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
1811     MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
1812     mainew->ajouter_mg_noeud(newno);
1813     vertices->change_nouveau_numero(newno->get_id()); // what does this do?
1814     }
1815     }
1816    
1817     LISTE_MG_SEGMENT::iterator itSeg;
1818     for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
1819     {
1820     //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
1821     //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
1822     //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
1823     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1824     //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
1825     MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
1826     //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
1827     if (faceseg->get_id()==5
1828     || aretver->get_id()==32 || aretver->get_id()==61|| aretver->get_id()==42|| aretver->get_id()==56 || aretver->get_id()==49||aretver->get_id()==27
1829     || aretver->get_id()==85 || aretver->get_id()==90|| aretver->get_id()==78|| aretver->get_id()==71 || aretver->get_id()==20||aretver->get_id()==13)
1830     {
1831     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1832     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1833     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
1834     //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
1835     mainew->ajouter_mg_segment(seg);
1836     }
1837     }
1838    
1839     LISTE_MG_TRIANGLE::iterator itTri;
1840     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
1841     {
1842     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
1843     if (facetri->get_id()==5)
1844     {
1845     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
1846     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
1847     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
1848     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
1849     }
1850     }
1851    
1852     //gestnew->enregistrer("scnsurf_smplwithpocket_3d_13mm_232mmdeform.magic");
1853     gestnew->enregistrer(outputfilename);
1854     }
1855     //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1856 sattarpa 596 void CRIAQOPERATORS::gnifformatmaker(char* elementxtoutfile,char* nodtxtoutfile,int meshno,char* magicfilename)
1857 sattarpa 554 {
1858     // gnifformatmaker
1859     MG_FILE gest(magicfilename);
1860     MG_MAILLAGE* mai;
1861     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1862     //MG_FILE gest((char*)"scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump.magic");
1863     // MG_MAILLAGE* mai=gest.get_mg_maillageid(2);
1864     //std::string fichiermaillage="CAD_output_elements.txt";
1865     //FILE* in=fopen(fichiermaillage.c_str(),"wt");
1866     //FILE* in1=fopen("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump_nodes.txt","wt");
1867     FILE* in1=fopen(nodtxtoutfile,"wt");
1868     //FILE* in1=fopen("scannedsurface_nodes.txt","wt");
1869     LISTE_MG_NOEUD::iterator itnod;
1870     unsigned long firstno=mai->get_premier_noeud(itnod)->get_id();
1871     firstno=firstno-1;
1872     cout<<firstno<<endl;
1873     int counter1=1;
1874     cout<<"mai->get_nb_mg_noeud()="<<mai->get_nb_mg_noeud()<<endl;
1875     for (MG_NOEUD* nod=mai->get_premier_noeud(itnod);nod!=NULL;nod=mai->get_suivant_noeud(itnod))
1876     {
1877     nod->change_nouveau_numero(counter1);
1878     //nod->change_id(counter1);
1879     //fprintf(in1,"%d node coordinate %le %le %le\n",nod->get_id(), nod->get_x(),nod->get_y(), nod->get_z() );
1880 sattarpa 658 //fprintf(in,"%d node coordinate %le %le %le\n",nod->get_id(), nod->get_x(),nod->get_y(), nod->get_z() );
1881     fprintf(in1,"%19d%32.9le%16.9le%8d\n",nod->get_nouveau_numero(), nod->get_x(),nod->get_y(),nod->get_nouveau_numero());
1882     fprintf(in1,"%7d%16.9le\n",nod->get_nouveau_numero(), nod->get_z());
1883 sattarpa 554 counter1++;
1884     }
1885     //FILE* in2=fopen("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump_elements.txt","wt"); // CAD text output elements
1886     FILE* in2=fopen(elementxtoutfile,"wt");
1887     //FILE* in2=fopen("scannedsurface_elements.txt","wt"); //scanned text output elements
1888     int counter=1;
1889     LISTE_MG_TRIANGLE::iterator ittri;
1890     for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittri);tri!=NULL;tri=mai->get_suivant_triangle(ittri))
1891     {
1892     fprintf(in2,"%10d 1%8d%8d%8d\n",counter,tri->get_noeud1()->get_nouveau_numero(),tri->get_noeud2()->get_nouveau_numero(),tri->get_noeud3()->get_nouveau_numero() );
1893     counter++;
1894     }
1895     fclose(in1);
1896     fclose(in2);
1897     }
1898     /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1899 sattarpa 596 void CRIAQOPERATORS::rmovscnprtfromdefrmdcad(char* outputfilename,int meshno,char* magicfilename)
1900 sattarpa 554 {
1901     MG_FILE gest(magicfilename);
1902     MG_MAILLAGE* mai;
1903     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
1904 sattarpa 596 MG_GESTIONNAIRE* newgest=new MG_GESTIONNAIRE();
1905 sattarpa 554 MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1906 sattarpa 596 newgest->ajouter_mg_geometrie(geonew);
1907 sattarpa 554 MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1908 sattarpa 596 newgest->ajouter_mg_maillage(mainew);
1909    
1910 sattarpa 554 LISTE_MG_NOEUD::iterator itnod;
1911     for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
1912     {
1913     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
1914     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
1915 sattarpa 596 //if(facevert->get_id()!=76 && topnd->get_dimension()!=3) //e1
1916     //if(facevert->get_id()!=5 && topnd->get_dimension()!=3) //e2
1917 sattarpa 554 {
1918     double x=vertices->get_x();
1919     double y=vertices->get_y();
1920     double z=vertices->get_z();
1921     MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
1922     mainew->ajouter_mg_noeud(newno);
1923     vertices->change_nouveau_numero(newno->get_id()); // what does this do?
1924     }
1925     }
1926 sattarpa 596
1927 sattarpa 554 LISTE_MG_SEGMENT::iterator itSeg;
1928     for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
1929     {
1930     MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
1931     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1932 sattarpa 596 //if (faceseg->get_id()!=76 && topseg->get_dimension()!=3)
1933     if(faceseg->get_id()!=5 && topseg->get_dimension()!=3) //e2
1934 sattarpa 554 {
1935     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1936     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1937     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
1938     mainew->ajouter_mg_segment(seg);
1939     }
1940     }
1941 sattarpa 596
1942 sattarpa 554 LISTE_MG_TRIANGLE::iterator itTri;
1943     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
1944     {
1945     MG_ELEMENT_TOPOLOGIQUE *toptri=triangle->get_lien_topologie();
1946     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
1947 sattarpa 596 // if (facetri->get_id()!=76 && toptri->get_dimension()!=3)
1948     if (facetri->get_id()!=5 && toptri->get_dimension()!=3) //e2
1949 sattarpa 554 {
1950     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
1951     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
1952     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
1953     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
1954     }
1955     }
1956 sattarpa 596 newgest->enregistrer(outputfilename);
1957 sattarpa 554 //gest.enregistrer("deformedcad_sansscannedpart.magic");
1958    
1959     }
1960     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1961 sattarpa 596 void CRIAQOPERATORS::msh2dmakerfrommsh3d(char* outputfilename,int meshno,char* magicfilename)
1962 sattarpa 554 {
1963     //to make a 2D mesh out of 3D mesh :)
1964     /*
1965     MG_MAILLAGE* mainew=mai->dupliquer(&gest);
1966     LISTE_MG_TETRA::iterator ittet;
1967     for(MG_TETRA* )
1968     {
1969    
1970     }
1971    
1972     cout<<"mainew->get_nb_mg_noeud()= "<<mainew->get_nb_mg_noeud()<<endl;
1973     cout<<"mainew->get_nb_mg_segment()= "<<mainew->get_nb_mg_segment()<<endl;
1974     cout<<"mainew->get_nb_mg_triangle()= "<<mainew->get_nb_mg_triangle()<<endl;
1975     cout<<"mainew->get_nb_mg_tetra()= "<<mainew->get_nb_mg_tetra()<<endl;
1976    
1977     LISTE_MG_TRIANGLE::iterator ittri;
1978     for (MG_TRIANGLE* triangle=mainew->get_premier_triangle(ittri);triangle!=NULL;triangle=mainew->get_suivant_triangle(ittri))
1979     {
1980     MG_ELEMENT_TOPOLOGIQUE *top=triangle->get_lien_topologie();
1981     int topdimension=top->get_dimension();
1982     if (topdimension==3)
1983     mainew->supprimer_mg_triangleid(triangle->get_id());
1984     }
1985    
1986     LISTE_MG_TETRA::iterator ittet;
1987     for (MG_TETRA* tetra=mainew->get_premier_tetra(ittet);tetra!=NULL;tetra=mainew->get_suivant_tetra(ittet))
1988     mainew->supprimer_mg_tetraid(tetra->get_id());
1989     //MG_FACE* face=geo->get_mg_faceid(5);
1990    
1991    
1992     LISTE_MG_SEGMENT::iterator itseg;
1993     for (MG_SEGMENT* segment=mainew->get_premier_segment(itseg);segment!=NULL;segment=mainew->get_suivant_segment(itseg))
1994     {
1995     MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
1996     if (faceseg->get_id()!=181)
1997     mainew->supprimer_mg_segmentid(segment->get_id());
1998    
1999     MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
2000     if(topseg->get_dimension()==3)
2001     cout<<"topseg->get_dimension()==3"<<endl;
2002     }
2003    
2004     LISTE_MG_NOEUD::iterator itnod;
2005     for (MG_NOEUD* vertices=mainew->get_premier_noeud(itnod);vertices!=NULL;vertices=mainew->get_suivant_noeud(itnod))
2006     {
2007     MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
2008     MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
2009     //if (facevert->get_id()!=181)
2010     if (topnd->get_dimension()==0)
2011     mainew->supprimer_mg_segmentid(vertices->get_id());
2012     }
2013    
2014     for (MG_TRIANGLE* triangle=mainew->get_premier_triangle(ittri);triangle!=NULL;triangle=mainew->get_suivant_triangle(ittri))
2015     {
2016     MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
2017     if (facetri->get_id()!=181)
2018     mainew->supprimer_mg_triangleid(triangle->get_id());
2019     }
2020     */
2021 sattarpa 596
2022     MG_FILE gest(magicfilename);
2023     MG_MAILLAGE* mai;
2024     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
2025 sattarpa 554 MG_GESTIONNAIRE* newgest=new MG_GESTIONNAIRE();
2026     MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
2027     newgest->ajouter_mg_geometrie(geonew);
2028     MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
2029     newgest->ajouter_mg_maillage(mainew);
2030    
2031     LISTE_MG_NOEUD::iterator itNod;
2032     for (MG_NOEUD* nd=mai->get_premier_noeud(itNod);nd!=NULL;nd=mai->get_suivant_noeud(itNod))
2033     {
2034     MG_ELEMENT_TOPOLOGIQUE *topnd=nd->get_lien_topologie();
2035     if(topnd->get_dimension()<3)
2036     {
2037     double x=nd->get_x();
2038     double y=nd->get_y();
2039     double z=nd->get_z();
2040     MG_NOEUD *newno=new MG_NOEUD(nd->get_lien_topologie(),x,y,z,nd->get_origine());
2041     mainew->ajouter_mg_noeud(newno);
2042     nd->change_nouveau_numero(newno->get_id());
2043     }
2044     }
2045     LISTE_MG_SEGMENT::iterator itSeg;
2046     for (MG_SEGMENT* segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
2047     {
2048     MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
2049     if(topseg->get_dimension()<3)
2050     {
2051     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
2052     MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
2053     MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
2054     mainew->ajouter_mg_segment(seg);
2055     }
2056     }
2057     LISTE_MG_TRIANGLE::iterator itTri;
2058     for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))
2059     {
2060     MG_ELEMENT_TOPOLOGIQUE *toptri=triangle->get_lien_topologie();
2061     if(toptri->get_dimension()<3)
2062     {
2063     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
2064     MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
2065     MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
2066     mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
2067     }
2068     }
2069 sattarpa 596 newgest->enregistrer(outputfilename);
2070 sattarpa 554 }
2071     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2072 sattarpa 596 void CRIAQOPERATORS::bumpmaker(char* outputfilename,int bumpndid, double bumptip,int meshno,char* magicfilename)
2073 sattarpa 554 {
2074     MG_FILE gest(magicfilename);
2075     //MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
2076     MG_MAILLAGE* mai;
2077     if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
2078     //MG_FILE gest((char*)"scnsurf_smplwithpocket_3d_13mm_64mmdeform.magic");
2079     //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
2080     //MG_MAILLAGE* mai=gest.get_mg_maillageid(2);
2081     MG_NOEUD* bptipnd;
2082     TPL_MAP_ENTITE<MG_NOEUD*> usednd;
2083     LISTE_MG_NOEUD::iterator itnod;
2084     for (MG_NOEUD* nd=mai->get_premier_noeud(itnod);nd!=NULL;nd=mai->get_suivant_noeud(itnod))
2085     {
2086     MG_ELEMENT_TOPOLOGIQUE *topnd=nd->get_lien_topologie();
2087     if (topnd->get_dimension() <3 && nd->get_id()==bumpndid)
2088     {
2089     bptipnd=nd;
2090     double new_xyz[3];
2091     new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+bumptip;
2092     nd->change_coord(new_xyz);
2093     usednd.ajouter(nd);
2094     }
2095     }
2096     TPL_MAP_ENTITE<MG_NOEUD*> usedndnd;
2097     TPL_MAP_ENTITE<MG_NOEUD*> neignd;
2098     for(int ineignd=0;ineignd<bptipnd->get_lien_segment()->get_nb();ineignd++)
2099     {
2100     MG_SEGMENT* segnei=bptipnd->get_lien_segment()->get(ineignd);
2101     if(segnei->get_noeud1()!=bptipnd) neignd.ajouter(segnei->get_noeud1());
2102     else if(segnei->get_noeud2()!=bptipnd) neignd.ajouter(segnei->get_noeud2());
2103     }
2104     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd;
2105     for (MG_NOEUD* nd=neignd.get_premier(itneignd);nd!=NULL;nd=neignd.get_suivant(itneignd))
2106     {
2107     double new_xyz[3];
2108     new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*3./4.);
2109     nd->change_coord(new_xyz);
2110     usednd.ajouter(nd);
2111     usedndnd.ajouter(nd);
2112     }
2113    
2114     TPL_MAP_ENTITE<MG_NOEUD*> neignd1;
2115     TPL_MAP_ENTITE<MG_NOEUD*> usedndnd1;
2116     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusednd;
2117     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusedndnd;
2118     for(MG_NOEUD* nd=usedndnd.get_premier(itusedndnd);nd!=NULL;nd=usedndnd.get_suivant(itusedndnd))
2119     {
2120     for(int ineignd=0;ineignd<nd->get_lien_segment()->get_nb();ineignd++)
2121     {
2122     MG_SEGMENT* segnei=nd->get_lien_segment()->get(ineignd);
2123     int sgnd1=0;
2124     int sgnd2=0;
2125     for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
2126     {
2127     if(segnei->get_noeud1()==nd1) sgnd1++;
2128     }
2129     for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
2130     {
2131     if(segnei->get_noeud2()==nd1) sgnd2++;
2132     }
2133     if(sgnd1==0) neignd1.ajouter(segnei->get_noeud1());
2134     else if(sgnd2==0) neignd1.ajouter(segnei->get_noeud2());
2135     }
2136     }
2137     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd1;
2138     for (MG_NOEUD* nd=neignd1.get_premier(itneignd1);nd!=NULL;nd=neignd1.get_suivant(itneignd1))
2139     {
2140     double new_xyz[3];
2141     new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*2./4.);
2142     nd->change_coord(new_xyz);
2143     usednd.ajouter(nd);
2144     usedndnd1.ajouter(nd);
2145     }
2146    
2147     TPL_MAP_ENTITE<MG_NOEUD*> neignd2;
2148     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusedndnd1;
2149     for(MG_NOEUD* nd=usedndnd1.get_premier(itusedndnd1);nd!=NULL;nd=usedndnd1.get_suivant(itusedndnd1))
2150     {
2151     for(int ineignd=0;ineignd<nd->get_lien_segment()->get_nb();ineignd++)
2152     {
2153     MG_SEGMENT* segnei=nd->get_lien_segment()->get(ineignd);
2154     int sgnd1=0;
2155     int sgnd2=0;
2156     for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
2157     {
2158     if(segnei->get_noeud1()==nd1) sgnd1++;
2159     }
2160     for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
2161     {
2162     if(segnei->get_noeud2()==nd1) sgnd2++;
2163     }
2164     if(sgnd1==0) neignd2.ajouter(segnei->get_noeud1());
2165     else if(sgnd2==0) neignd2.ajouter(segnei->get_noeud2());
2166     }
2167     }
2168     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd2;
2169     for (MG_NOEUD* nd=neignd2.get_premier(itneignd2);nd!=NULL;nd=neignd2.get_suivant(itneignd2))
2170     {
2171     double new_xyz[3];
2172     new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*1./4.);
2173     nd->change_coord(new_xyz);
2174     }
2175    
2176    
2177     //gest.enregistrer("scnsurf_smplwithpocket_3d_13mm_64mmdeform_withbump5mm.magic");
2178     gest.enregistrer(outputfilename);
2179     }
2180 sattarpa 596 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2181     int CRIAQOPERATORS::import_triangulation_gnif(char* outputfilename,char* triangulationfile)
2182     {
2183     FILE* in=fopen(triangulationfile,"rt");
2184     if (in==NULL) return 0;
2185 sattarpa 554
2186    
2187 sattarpa 596 MG_GESTIONNAIRE gesttmp;
2188     MG_MAILLAGE* mai=new MG_MAILLAGE(NULL);
2189     gesttmp.ajouter_mg_maillage(mai);
2190     char mess[500];
2191 sattarpa 554
2192 sattarpa 596 double xmax=-1e308;
2193     double xmin=1e308;
2194     double ymax=-1e308;
2195     double ymin=1e308;
2196     double zmax=-1e308;
2197     double zmin=1e308;
2198     while (!feof(in))
2199     {
2200     double x,y,z;
2201     char nom[20];
2202     char *res=fgets(mess,500,in);
2203     if (feof(in)) break;
2204 sattarpa 655 if(res!=NULL)
2205     {
2206 sattarpa 658 sscanf(mess,"%le %le %le",&x,&y,&z);
2207 sattarpa 596 if (x<xmin) xmin=x;
2208     if (x>xmax) xmax=x;
2209     if (y<ymin) ymin=y;
2210     if (y>ymax) ymax=y;
2211     if (z<zmin) zmin=z;
2212     if (z>zmax) zmax=z;
2213     if (feof(in)) return 1;
2214     MG_NOEUD* noeud1=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
2215     noeud1->change_nouveau_numero(-1);
2216     res=fgets(mess,500,in);
2217     if (feof(in)) break;
2218 sattarpa 658 sscanf(mess,"%le %le %le",&x,&y,&z);
2219 sattarpa 596 if (x<xmin) xmin=x;
2220     if (x>xmax) xmax=x;
2221     if (y<ymin) ymin=y;
2222     if (y>ymax) ymax=y;
2223     if (z<zmin) zmin=z;
2224     if (z>zmax) zmax=z;
2225     if (feof(in)) return 1;
2226     res=fgets(mess,500,in);
2227     if (feof(in)) break;
2228     MG_NOEUD* noeud2=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
2229     noeud2->change_nouveau_numero(-1);
2230 sattarpa 658 sscanf(mess,"%le %le %le",&x,&y,&z);
2231 sattarpa 596 if (x<xmin) xmin=x;
2232     if (x>xmax) xmax=x;
2233     if (y<ymin) ymin=y;
2234     if (y>ymax) ymax=y;
2235     if (z<zmin) zmin=z;
2236     if (z>zmax) zmax=z;
2237     if (feof(in)) return 1;
2238     MG_NOEUD* noeud3=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
2239     noeud3->change_nouveau_numero(-1);
2240     mai->ajouter_mg_triangle(NULL,noeud1,noeud2,noeud3,TRIANGULATION);
2241 sattarpa 655 }
2242 sattarpa 596 }
2243     fclose(in);
2244     //to merge the nodes of each triangulation which are close (the same)
2245     double eps=1e-4;
2246     BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
2247     double rayon=boite.get_rayon()*1.1;
2248     double x=boite.get_xcentre();
2249     double y=boite.get_ycentre();
2250     double z=boite.get_zcentre();
2251     boite.reinit(x-rayon,y-rayon,z-rayon,x+rayon,y+rayon,z+rayon);
2252     TPL_GRILLE<MG_NOEUD*> grille;
2253     int nb_noeud=mai->get_nb_mg_noeud();
2254     int pas=(int)(pow(nb_noeud,0.33333333333333333))+1;
2255     grille.initialiser(boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax(),pas,pas,pas);
2256     std::vector<MG_NOEUD*> lst_noeud;
2257     LISTE_MG_NOEUD::iterator it;
2258     int i=0;
2259     for (MG_NOEUD* noeud=mai->get_premier_noeud(it);noeud!=NULL;noeud=mai->get_suivant_noeud(it))
2260     {
2261     //MG_NOEUD* noeud=mai->get_mg_noeud(i);
2262     grille.inserer(noeud);
2263     noeud->change_nouveau_numero(i);
2264     lst_noeud.insert(lst_noeud.end(),noeud);
2265     i++;
2266     }
2267 sattarpa 554
2268 sattarpa 596 int nb_cell=grille.get_nb_cellule();
2269     for (int i=0;i<nb_cell;i++)
2270     {
2271     TPL_CELLULE_GRILLE<MG_NOEUD*> *cell=grille.get_cellule(i);
2272     int nb_noeud_cell=cell->get_nb_entite();
2273     for (int j=0;j<nb_noeud_cell;j++)
2274     {
2275     MG_NOEUD* noeud1=cell->get_entite(j);
2276     // if (noeud1->get_nouveau_numero()!=(-1)) continue;
2277     // noeud1->change_nouveau_numero(noeud1->get_id());
2278     double *xyz1=noeud1->get_coord();
2279     int nb=noeud1->get_lien_segment()->get_nb();
2280     double longref=0.;
2281     for (int jj=0;jj<nb;jj++)
2282     longref=longref+noeud1->get_lien_segment()->get(jj)->get_longueur();
2283     longref=longref/nb;
2284     for (int k=j+1;k<nb_noeud_cell;k++)
2285     {
2286     MG_NOEUD* noeud2=cell->get_entite(k);
2287     double *xyz2=noeud2->get_coord();
2288     OT_VECTEUR_3D vec(xyz1,xyz2);
2289     double longueur=vec.get_longueur();
2290     if (longueur<eps*longref)
2291     {
2292     lst_noeud[noeud2->get_nouveau_numero()]=noeud1;
2293     }
2294     }
2295     }
2296     }
2297     MG_GESTIONNAIRE gest;
2298     MG_MAILLAGE* mai2=new MG_MAILLAGE(NULL);
2299     gest.ajouter_mg_maillage(mai2);
2300     for (int i=0;i<nb_noeud;i++)
2301     {
2302     MG_NOEUD* noeud=lst_noeud[i];
2303     MG_NOEUD* nvnoeud;
2304     if (noeud->get_nouveau_numero()==i) nvnoeud=mai2->ajouter_mg_noeud(NULL,noeud->get_x(),noeud->get_y(),noeud->get_z(),TRIANGULATION);
2305     else nvnoeud=lst_noeud[noeud->get_nouveau_numero()];
2306     lst_noeud[i]=nvnoeud;
2307     }
2308     int nb_tri=mai->get_nb_mg_triangle();
2309     for (int i=0;i<nb_tri;i++)
2310     {
2311     MG_TRIANGLE* tri=mai->get_mg_triangle(i);
2312     MG_NOEUD* noeud1=tri->get_noeud1();
2313     MG_NOEUD* noeud2=tri->get_noeud2();
2314     MG_NOEUD* noeud3=tri->get_noeud3();
2315     mai2->ajouter_mg_triangle(NULL,lst_noeud[noeud1->get_nouveau_numero()],lst_noeud[noeud2->get_nouveau_numero()],lst_noeud[noeud3->get_nouveau_numero()],TRIANGULATION);
2316     }
2317     gest.enregistrer(outputfilename);
2318     return 1;
2319     }