ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/app/criaqoperation/src/criaqoperators.cpp
Revision: 655
Committed: Fri Feb 20 20:09:52 2015 UTC (10 years, 5 months ago) by sattarpa
File size: 108391 byte(s)
Log Message:
modifier le lecture de *.txt fichier.

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