ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/app/criaqoperation/src/criaqoperators.cpp
Revision: 658
Committed: Wed Mar 4 20:35:01 2015 UTC (10 years, 5 months ago) by sattarpa
File size: 105696 byte(s)
Log Message:
debuger criaqoperation et ajouter point_insert avec types.

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