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, 6 months ago) by sattarpa
File size: 108391 byte(s)
Log Message:
modifier le lecture de *.txt fichier.

File Contents

# Content
1 #include <stdio.h>
2 #include "mg_file.h"
3 #include "criaqoperators.h"
4 #include "tpl_grille.h"
5 #include "ot_mathematique.h"
6 #include <math.h>
7 #include "tpl_octree.h"
8 #include "fem_maillage_outils.h"
9 #include "mg_maillage_outils.h"
10 #include <sstream>
11 #include <string.h>
12
13 CRIAQOPERATORS::CRIAQOPERATORS():MAILLEUR(false)
14 {
15 }
16 CRIAQOPERATORS::~CRIAQOPERATORS()
17 {
18 }
19 void CRIAQOPERATORS::meshdistance_compare(char* referencemagicfilename,int refmeshno,char* comparemagicfilename,int compmeshno,char* outputcomparefile)
20 {
21 MG_FILE refgest(referencemagicfilename);
22 MG_MAILLAGE* refmgmai;
23 if (refmeshno==0) refmgmai=refgest.get_mg_maillage(refmeshno); else refmgmai=refgest.get_mg_maillageid(refmeshno);
24 MG_FILE compgest(comparemagicfilename);
25 MG_MAILLAGE* compmgmai;
26 if (compmeshno==0) compmgmai=compgest.get_mg_maillage(compmeshno); else compmgmai=compgest.get_mg_maillageid(compmeshno);
27 /*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 MG_SOLUTION* mgsol=new MG_SOLUTION(refmgmai,1,(char*)"distance.sol",1,(char*)"Normaldistance",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
70 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 }
129 refgest.enregistrer(outputcomparefile);
130 }
131 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 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 {
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 if(feof(in1))
245 cout<<"correspondfile has not been well written"<<endl;
246 res=fgets(mess,500,in1);
247 if(feof(in1))
248 cout<<"correspondfile has not been well written"<<endl;
249 int cadmeshno;
250 int defcadmeshno;
251 int nb1=sscanf(mess,"%ld %ld",&cadmeshno,&defcadmeshno);
252 while (!feof(in1))
253 {
254 char chaine[500];
255 res=fgets(chaine,500,in1);
256 if(res!=NULL)
257 {
258 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 }
267 }
268 fclose(in1);
269
270 cout<<"cadmeshno= "<<cadmeshno<<" defcadmeshno= "<<defcadmeshno<<endl;
271 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
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 multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_cadmgmai;
328 multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_cadmgmai;
329 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_cadmgmai;
330 multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_defcadmgmai;
331 multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_defcadmgmai;
332 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_defcadmgmai;
333 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 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
355 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
369 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 {
377 lstkmin_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,cadmgnd) );
378 lstkmin_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmin).second));
379 }
380 }
381 for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
382 {
383 if(mess==(*itkmax).first)
384 {
385 lstkmax_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,cadmgnd) );
386 lstkmax_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmax).second));
387 }
388 }
389 for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
390 {
391 if(mess==(*itkgau).first)
392 {
393 lstkgau_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,cadmgnd) );
394 lstkgau_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkgau).second));
395 }
396 }
397 }
398 //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 // 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 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
568 cout<<"calcul courbure defcadmgmai"<<endl;
569 // calcul courbure defcadmgmai ********************************
570
571 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 {
590 lstkmin_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,defcadmgnd) );
591 lstkmin_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmin).second));
592 }
593 }
594 for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
595 {
596 if(mess==(*itkmax).first)
597 {
598 lstkmax_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,defcadmgnd) );
599 lstkmax_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmax).second));
600 }
601 }
602 for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
603 {
604 if(mess==(*itkgau).first)
605 {
606 lstkgau_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,defcadmgnd) );
607 lstkgau_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkgau).second));
608 }
609 }
610 }
611 //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 // 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 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 multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_diff;
801 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_cadmgmai;
802 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_cadmgmai;
803 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_cadmgmai;
804 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_defcadmgmai;
805 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_defcadmgmai;
806 multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_defcadmgmai;
807 cout<<"lstkmin_cadmgmai.size()= "<<lstkmin_cadmgmai.size()<<endl;
808 cout<<"lstkmin_defcadmgmai.size()= "<<lstkmin_defcadmgmai.size()<<endl;
809 cout<<"kmin diff. calc. started"<<endl;
810
811 for(itlstkmin_cadmgmai=lstkmin_cadmgmai.begin();itlstkmin_cadmgmai!=lstkmin_cadmgmai.end();itlstkmin_cadmgmai++)
812 {
813 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 cout<<"kmax diff. calc. started"<<endl;
828 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 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 // 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 // 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
1054 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 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 char* res=fgets(chaine,500,in);
1081 if(res!=NULL)
1082 {
1083 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 }
1097 }
1098 fclose(in);
1099 if (relet_curvature<1)
1100 {
1101 //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
1125 cout<<"for Kmin"<<endl;
1126 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 octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1139 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 if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
1150 lstpinsert.erase(itlstpinsert);
1151 }
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 octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
1170 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 }
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 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 //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 FILE* in2=fopen(out_gnifinspointfile_removondefect,"wt");
1280 int countfuck=0;
1281 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 fclose(in2);
1288 }
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 MG_FILE gest(magicfilename);
1292 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 nod->get_mg_element_maillage()->get_id();
1302 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 double crit=vm_coef*(maxi+mini)/2.;
1313 //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 char* res=fgets(chaine,500,in);
1322 if(res!=NULL)
1323 {
1324 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 }
1338 }
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 //double search_radius=2.; //sasan this should be as a parameter
1361 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 double search_radius1=search_radius;
1379 while(neindnb==0)
1380 {
1381 octree.rechercher(pdefect->get_x(),pdefect->get_y(),pdefect->get_z(),search_radius1,lstneipdefect);
1382 neindnb=lstneipdefect.get_nb();
1383 search_radius1=search_radius1+0.1*search_radius;
1384 }
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 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
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 || 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 || sometver->get_id()==352|| sometver->get_id()==354|| sometver->get_id()==10|| sometver->get_id()==26*/
1455
1456 )
1457 {
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 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 {
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 if (facetri->get_id()==76)
1502 {
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 void CRIAQOPERATORS::surfmaker_e2(char* outputfilename,int meshno,char* magicfilename)
1532 {
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 void CRIAQOPERATORS::surfmaker_e2t(char* outputfilename,int meshno,char* magicfilename)
1606 {
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 void CRIAQOPERATORS::surfmaker_e3(char* outputfilename,int meshno,char* magicfilename)
1697 {
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 void CRIAQOPERATORS::gnifformatmaker(char* elementxtoutfile,char* nodtxtoutfile,int meshno,char* magicfilename)
1775 {
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 void CRIAQOPERATORS::rmovscnprtfromdefrmdcad(char* outputfilename,int meshno,char* magicfilename)
1818 {
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 MG_GESTIONNAIRE* newgest=new MG_GESTIONNAIRE();
1823 MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
1824 newgest->ajouter_mg_geometrie(geonew);
1825 MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
1826 newgest->ajouter_mg_maillage(mainew);
1827
1828 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 //if(facevert->get_id()!=76 && topnd->get_dimension()!=3) //e1
1834 //if(facevert->get_id()!=5 && topnd->get_dimension()!=3) //e2
1835 {
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
1845 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 //if (faceseg->get_id()!=76 && topseg->get_dimension()!=3)
1851 if(faceseg->get_id()!=5 && topseg->get_dimension()!=3) //e2
1852 {
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
1860 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 // if (facetri->get_id()!=76 && toptri->get_dimension()!=3)
1866 if (facetri->get_id()!=5 && toptri->get_dimension()!=3) //e2
1867 {
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 newgest->enregistrer(outputfilename);
1875 //gest.enregistrer("deformedcad_sansscannedpart.magic");
1876
1877 }
1878 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1879 void CRIAQOPERATORS::msh2dmakerfrommsh3d(char* outputfilename,int meshno,char* magicfilename)
1880 {
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
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 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 newgest->enregistrer(outputfilename);
1988 }
1989 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1990 void CRIAQOPERATORS::bumpmaker(char* outputfilename,int bumpndid, double bumptip,int meshno,char* magicfilename)
1991 {
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 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2099 int CRIAQOPERATORS::import_triangulation_gnif(char* outputfilename,char* triangulationfile)
2100 {
2101 FILE* in=fopen(triangulationfile,"rt");
2102 if (in==NULL) return 0;
2103
2104
2105 MG_GESTIONNAIRE gesttmp;
2106 MG_MAILLAGE* mai=new MG_MAILLAGE(NULL);
2107 gesttmp.ajouter_mg_maillage(mai);
2108 char mess[500];
2109
2110 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 if(res!=NULL)
2123 {
2124 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 }
2160 }
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
2186 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 }