ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/app/criaqoperation/src/criaqoperators.cpp
Revision: 650
Committed: Wed Feb 11 00:03:21 2015 UTC (10 years, 6 months ago) by sattarpa
File size: 108124 byte(s)
Log Message:
proximity parameters added to "mailleur2d_ins_noeud"
relative curvature calculation added to "criaqoperators"

File Contents

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