ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/occ_surface.cpp
Revision: 1075
Committed: Tue Aug 10 17:02:54 2021 UTC (4 years ago) by francois
File size: 35250 byte(s)
Log Message:
suppression de warning avec le dernier compilateur

File Contents

# Content
1 //---------------------------------------------------------------------------
2 //------------------------------------------------------------
3 //------------------------------------------------------------
4 // MAGiC
5 // Jean Christophe Cuilli�re et Vincent FRANCOIS
6 // D�partement de G�nie M�canique - UQTR
7 //------------------------------------------------------------
8 // Le projet MAGIC est un projet de recherche du d�partement
9 // de g�nie m�canique de l'Universit� du Qu�bec �
10 // Trois Rivi�res
11 // Les librairies ne peuvent �tre utilis�es sans l'accord
12 // des auteurs (contact : francois@uqtr.ca)
13 //------------------------------------------------------------
14 //------------------------------------------------------------
15 //
16 // OCC_Surface.cpp
17 //
18 //------------------------------------------------------------
19 //------------------------------------------------------------
20 // COPYRIGHT 2000
21 // Version du 02/03/2006 � 11H22
22 //------------------------------------------------------------
23 //------------------------------------------------------------
24
25 #pragma hdrstop
26 #include "gestionversion.h"
27 #ifdef BREP_OCC
28
29 #include "occ_surface.h"
30 #include <gp_Pnt.hxx>
31 #include <gp_Vec.hxx>
32 #include "GeomAPI_ProjectPointOnSurf.hxx"
33 #include <ShapeAnalysis_Surface.hxx>
34 #include <BRep_Tool.hxx>
35 #include <BRepAdaptor_Surface.hxx>
36 //***********************************
37 #include <BRep_Tool.hxx>
38 #include <Poly_Triangulation.hxx>
39 #include <Geom_Plane.hxx>
40 #include <gp_Pln.hxx>
41 #include <Geom_CylindricalSurface.hxx>
42 #include <gp_Cylinder.hxx>
43 #include <Geom_ConicalSurface.hxx>
44 #include <gp_Cone.hxx>
45 #include <Geom_SphericalSurface.hxx>
46 #include <IGESSolid_SphericalSurface.hxx>
47 #include <gp_Sphere.hxx>
48 #include "Geom_ToroidalSurface.hxx"
49 #include <gp_Torus.hxx>
50 #include <Geom_BSplineSurface.hxx>
51 #include <Geom_BezierSurface.hxx>
52 #include <GeomConvert.hxx>
53 #include <BRepBuilderAPI_NurbsConvert.hxx>
54 #include <BRepLib_FindSurface.hxx>
55 #include "mg_gestionnaire.h"
56 #include "constantegeo.h"
57 #include "ot_mathematique.h"
58 #include <Geom_RectangularTrimmedSurface.hxx>
59 #include <BRepClass_FaceClassifier.hxx>
60 #include <gp_Pnt2d.hxx>
61 #include <Geom_SurfaceOfRevolution.hxx>
62
63 #pragma package(smart_init)
64
65
66
67 class NOEUDARETE
68 {
69 public:
70 MG_NOEUD* no;
71 MG_ARETE* are;
72 double t;
73 };
74
75
76
77 OCC_SURFACE::OCC_SURFACE(unsigned long num, TopoDS_Face srf, OCC_FONCTION* fonc):MG_SURFACE(num),face(srf), fonction1(fonc)
78 {
79 surface=BRep_Tool::Surface(face);
80 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
81 Handle(Standard_Type) type=surface->DynamicType();
82 if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
83 {
84 Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
85 surface = RTSurface->BasisSurface();
86 }
87 double u1;
88 double u2;
89 double v1;
90 double v2;
91 surface->Bounds(u1,u2,v1, v2);
92 u_min=u1;
93 u_max=u2;
94 v_min=v1;
95 v_max=v2;
96 estperiodeu=surface->IsUClosed();
97 estperiodev=surface->IsVClosed();
98 if (surface->IsUPeriodic()) periode_u=surface->UPeriod();
99 else if (surface->IsUClosed()) periode_u= u_max-u_min;
100 else periode_u=0;
101 if (surface->IsVPeriodic()) periode_v=surface->VPeriod();
102 else if (surface->IsVClosed()) periode_v= v_max-v_min;
103 else periode_v=0;
104 if (type==STANDARD_TYPE(Geom_BSplineSurface))
105 analyse_bspline();
106
107 //**************************************************************
108 // ShapeAnalysis::GetFaceUVBounds(face,u_min,u_max,v_min,v_max);
109 //**************************************************************
110
111 }
112
113 OCC_SURFACE::OCC_SURFACE(TopoDS_Face srf, OCC_FONCTION* fonc):MG_SURFACE(),face(srf), fonction1(fonc)
114 {
115 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
116 surface=BRep_Tool::Surface(face);
117 Handle(Standard_Type) type=surface->DynamicType();
118 if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
119 {
120 Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
121 surface = RTSurface->BasisSurface();
122 }
123 double u1;
124 double u2;
125 double v1;
126 double v2;
127 surface->Bounds(u1,u2,v1, v2);
128 u_min=u1;
129 u_max=u2;
130 v_min=v1;
131 v_max=v2;
132 estperiodeu=surface->IsUClosed();
133 estperiodev=surface->IsVClosed();
134 if (surface->IsUPeriodic()) periode_u=surface->UPeriod();
135 else if (surface->IsUClosed()) periode_u= u_max-u_min;
136 else periode_u=0;
137 if (surface->IsVPeriodic()) periode_v=surface->VPeriod();
138 else if (surface->IsVClosed()) periode_v= v_max-v_min;
139 else periode_v=0;
140 if (type==STANDARD_TYPE(Geom_BSplineSurface))
141 analyse_bspline();
142 //**************************************************************
143 // ShapeAnalysis::GetFaceUVBounds(face,u_min,u_max,v_min,v_max);
144 //**************************************************************
145
146
147 // std::cout << " " <<std::endl;
148 // Handle(Standard_Type) type=surface->DynamicType();
149 // if (type==STANDARD_TYPE(Geom_Plane)) std::cout << "plan" << std::endl;
150 // else if (type==STANDARD_TYPE(Geom_CylindricalSurface)) std::cout << "Cylindre" << std::endl;
151 // else if (type==STANDARD_TYPE(Geom_ConicalSurface)) std::cout << "Cone" << std::endl;
152 // else if (type==STANDARD_TYPE(Geom_SphericalSurface)) std::cout << "Sphere" << std::endl;
153 // else if (type==STANDARD_TYPE(Geom_ToroidalSurface)) std::cout << "Tore" << std::endl;
154 // else if (type==STANDARD_TYPE(Geom_BSplineSurface)) std::cout << "Bspline surf" << std::endl;
155 // else std::cout << type->Name() << std::endl;
156 //
157 // std::cout << surface->IsUClosed()<< " " <<surface->IsVClosed() <<std::endl;;
158 // std::cout << surface->IsUPeriodic() <<" " <<surface->IsVPeriodic() << std::endl;;
159 // if (surface->IsUPeriodic()) std::cout << "Pu=" << surface->UPeriod() << std::endl;;
160 // if (surface->IsVPeriodic()) std::cout << "Pv=" << surface->VPeriod() << std::endl;;
161 }
162
163 OCC_SURFACE::OCC_SURFACE(OCC_SURFACE& mdd):MG_SURFACE(mdd),face(mdd.face), fonction1(mdd.fonction1)
164 {
165 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
166 surface=BRep_Tool::Surface(face);
167 Handle(Standard_Type) type=surface->DynamicType();
168 if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
169 {
170 Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
171 surface = RTSurface->BasisSurface();
172 }
173 u_min=mdd.get_umin();
174 u_max=mdd.get_umax();
175 v_min=mdd.get_vmin();
176 v_max=mdd.get_vmax();
177 estperiodeu=mdd.estperiodeu;
178 estperiodev=mdd.estperiodev;
179 periode_u=mdd.periode_u;
180 periode_v=mdd.periode_v;
181 }
182 OCC_SURFACE::~OCC_SURFACE()
183 {
184 }
185 void OCC_SURFACE::evaluer(double *uv,double *xyz)
186 {
187
188 //const Handle(Geom_Surface) &surface=BRep_Tool::Surface(face);
189 gp_Pnt P;
190 double u=uv[0];
191 double v=uv[1];
192
193 surface->D0(u,v,P);
194
195 xyz[0]=P.X();
196 xyz[1]=P.Y();
197 xyz[2]=P.Z();
198 }
199 void OCC_SURFACE::deriver(double *uv,double *xyzdu, double *xyzdv)
200 {
201 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
202 double u=uv[0];
203 double v=uv[1];
204 gp_Vec D1U;
205 gp_Vec D1V;
206 gp_Pnt P;
207
208 surface->D1(u,v,P,D1U,D1V);
209 xyzdu[0]=D1U.X();
210 xyzdu[1]=D1U.Y();
211 xyzdu[2]=D1U.Z();
212
213 xyzdv[0]=D1V.X();
214 xyzdv[1]=D1V.Y();
215 xyzdv[2]=D1V.Z();
216
217 }
218 void OCC_SURFACE::deriver_seconde(double *uv,double* xyzduu,double* xyzduv,double* xyzdvv,double *xyz, double *xyzdu, double *xyzdv)
219 {
220 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
221 double u=uv[0];
222 double v=uv[1];
223 gp_Pnt P;
224 gp_Vec D1U;
225 gp_Vec D1V;
226 gp_Vec D2U;
227 gp_Vec D2V;
228 gp_Vec D2UV;
229
230 surface->D2(u,v,P,D1U,D1V,D2U,D2V,D2UV);
231
232 xyz[0]=P.X();
233 xyz[1]=P.Y();
234 xyz[2]=P.Z();
235
236 xyzdu[0]=D1U.X();
237 xyzdu[1]=D1U.Y();
238 xyzdu[2]=D1U.Z();
239
240 xyzdv[0]=D1V.X();
241 xyzdv[1]=D1V.Y();
242 xyzdv[2]=D1V.Z();
243
244 xyzduu[0]=D2U.X();
245 xyzduu[1]=D2U.Y();
246 xyzduu[2]=D2U.Z();
247
248 xyzdvv[0]=D2V.X();
249 xyzdvv[1]=D2V.Y();
250 xyzdvv[2]=D2V.Z();
251
252 xyzduv[0]=D2UV.X();
253 xyzduv[1]=D2UV.Y();
254 xyzduv[2]=D2UV.Z();
255
256 }
257
258 void OCC_SURFACE::inverser(double *uv,double *xyz,double precision)
259 {
260 // Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
261 double u=xyz[0];
262 double v=xyz[1];
263 double w=xyz[2];
264 gp_Pnt P(u,v,w);
265
266 //ShapeAnalysis_Surface SAS(surface);
267 //gp_Pnt2d pnt2d=SAS.ValueOfUV(P, precision);
268 GeomAPI_ProjectPointOnSurf PPS(P,surface, precision);
269 double dist=PPS.LowerDistance();
270 if (dist>precision)
271 inverser2(uv,xyz,precision);
272 else
273 {
274 double UU, VV;
275 PPS.LowerDistanceParameters(UU,VV);
276 uv[0]=UU;
277 uv[1]=VV;
278 }
279 }
280
281 bool OCC_SURFACE::est_sur_surface(double* xyz, double precision)
282 {
283 #ifdef ALL_OCC
284 GeomAPI_ProjectPointOnSurf projecteur;
285 double u0, u1;
286 double v0, v1;
287 ShapeAnalysis::GetFaceUVBounds(face,u0,u1,v0,v1);
288 gp_Pnt point(xyz[0],xyz[1],xyz[2]);
289 projecteur.Init(point,surface,u0,u1,v0,v1);
290 if(projecteur.NbPoints()==0) return false;
291 if(projecteur.LowerDistance()<precision) return true;
292 return false;
293 #else
294 std::cerr << "OCC_SURFACE::est_sur_surface(double* xyz, double precision) --> Non disponible" << std::endl;
295 #endif
296
297 }
298
299 int OCC_SURFACE::est_periodique_u(void)
300 {
301 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
302 //return surface->IsUClosed();
303 return estperiodeu;
304 }
305
306 int OCC_SURFACE::est_periodique_v(void)
307 {
308 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
309 // return surface->IsVPeriodic();
310 // return surface->IsVClosed();
311 return estperiodev;
312 }
313 double OCC_SURFACE::get_periode_u(void)
314 {
315 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
316 //if (surface->IsUPeriodic()) return surface->UPeriod();
317 //if (surface->IsUClosed()) return u_max-u_min;
318 //return 1;
319 return periode_u;
320 }
321
322 double OCC_SURFACE::get_periode_v(void)
323 {
324 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
325 //if (surface->IsVPeriodic()) return surface->VPeriod();
326 //if (surface->IsVClosed()) return v_max-v_min;
327 //return 0;
328 return periode_v;
329
330 }
331
332 void OCC_SURFACE::enregistrer(std::ostream& o,double version)
333 {
334 if(fonction1->get_version()=="OCCV2017")
335 {
336 o <<"%"<<get_id()<< "=SURFACE_OCC("<< get_idoriginal() << ");" << std::endl;
337 }
338 else
339 {
340 o <<"%"<<get_id()<< "=SURFACE_OCC("<< fonction1->GetID(face)<< ");" << std::endl;
341 }
342 }
343 int OCC_SURFACE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
344 {
345 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
346 Handle(Standard_Type) type=surface->DynamicType();
347 //******plan
348 if (type==STANDARD_TYPE(Geom_Plane))
349 {
350 Handle(Geom_Plane) Pln=Handle(Geom_Plane)::DownCast(surface);
351 gp_Pln plan=Pln->Pln();
352
353 double origine[3];
354 gp_Pnt centre=plan.Location();
355
356 origine[0]=centre.X();
357 origine[1]=centre.Y();
358 origine[2]=centre.Z();
359
360 double normal[3];
361 gp_Ax1 axe=plan.Axis();
362 gp_Dir direction=axe.Direction();
363
364 normal[0]=direction.X();
365 normal[1]=direction.Y();
366 normal[2]=direction.Z();
367
368 param.ajouter(origine[0]);
369 param.ajouter(origine[1]);
370 param.ajouter(origine[2]);
371 param.ajouter(normal[0]);
372 param.ajouter(normal[1]);
373 param.ajouter(normal[2]);
374
375 return MGCo_PLAN;
376 }
377 //******cylindre
378 if (type==STANDARD_TYPE(Geom_CylindricalSurface))
379 {
380 Handle(Geom_CylindricalSurface) cylinder=Handle(Geom_CylindricalSurface)::DownCast(surface);
381 gp_Cylinder cylin=cylinder->Cylinder();
382
383 double origine[3];
384 gp_Pnt centre=cylin.Location();
385 origine[0]=centre.X();
386 origine[1]=centre.Y();
387 origine[2]=centre.Z();
388 double directionX[3];
389 double direction[3];
390 double rayon;
391 rayon=cylin.Radius();
392
393 gp_Ax3 repere=cylin.Position();
394 gp_Dir dirX=repere.XDirection();
395 gp_Dir dirZ=repere.Direction();
396 directionX[0]=dirX.X();
397 directionX[1]=dirX.Y();
398 directionX[2]=dirX.Z();
399 direction[0]=dirZ.X();
400 direction[1]=dirZ.Y();
401 direction[2]=dirZ.Z();
402
403
404
405 param.ajouter(origine[0]);
406 param.ajouter(origine[1]);
407 param.ajouter(origine[2]);
408 param.ajouter(directionX[0]);
409 param.ajouter(directionX[1]);
410 param.ajouter(directionX[2]);
411 param.ajouter(direction[0]);
412 param.ajouter(direction[1]);
413 param.ajouter(direction[2]);
414 param.ajouter(rayon);
415 param.ajouter(rayon);
416
417 return MGCo_CYLINDRE;
418 }
419 //******Cone
420 if (type==STANDARD_TYPE(Geom_ConicalSurface))
421 {
422 Handle(Geom_ConicalSurface) cone=Handle(Geom_ConicalSurface)::DownCast(surface);
423 gp_Cone con=cone->Cone();
424
425 double origine[3];
426 gp_Pnt centre=con.Location();
427 origine[0]=centre.X();
428 origine[1]=centre.Y();
429 origine[2]=centre.Z();
430 double direction[3];
431 double directionX[3];
432 double rayon;
433 rayon=con.RefRadius();
434 double angle;
435 angle=con.SemiAngle();
436
437 gp_Ax3 repere=con.Position();
438 gp_Dir dirX=repere.XDirection();
439 gp_Dir dirZ=repere.Direction();
440 directionX[0]=dirX.X();
441 directionX[1]=dirX.Y();
442 directionX[2]=dirX.Z();
443 direction[0]=dirZ.X();
444 direction[1]=dirZ.Y();
445 direction[2]=dirZ.Z();
446 param.ajouter(origine[0]);
447 param.ajouter(origine[1]);
448 param.ajouter(origine[2]);
449 param.ajouter(directionX[0]);
450 param.ajouter(directionX[1]);
451 param.ajouter(directionX[2]);
452 param.ajouter(direction[0]);
453 param.ajouter(direction[1]);
454 param.ajouter(direction[2]);
455 param.ajouter(rayon);
456 param.ajouter(rayon);
457 param.ajouter(cos(angle));
458 param.ajouter(sin(angle));
459
460 return MGCo_CONE;
461
462 }
463 //*****Sphere
464 if (type==STANDARD_TYPE(Geom_SphericalSurface))
465 {
466 Handle(Geom_SphericalSurface) sphere=Handle(Geom_SphericalSurface)::DownCast(surface);
467 gp_Sphere sph=sphere->Sphere();
468
469 double origine[3];
470 gp_Pnt centre=sph.Location();
471 origine[0]=centre.X();
472 origine[1]=centre.Y();
473 origine[2]=centre.Z();
474 double rayon;
475 rayon=sph.Radius();
476 param.ajouter(origine[0]);
477 param.ajouter(origine[1]);
478 param.ajouter(origine[2]);
479 param.ajouter(rayon);
480
481 return MGCo_SPHERE;
482
483 }
484 //****** Tore
485 if (type==STANDARD_TYPE(Geom_ToroidalSurface))
486 {
487 Handle(Geom_ToroidalSurface) tore=Handle(Geom_ToroidalSurface)::DownCast(surface);
488 gp_Torus tro=tore->Torus();
489 double origine[3];
490 gp_Pnt centre=tro.Location();
491 origine[0]=centre.X();
492 origine[1]=centre.Y();
493 origine[2]=centre.Z();
494 double direction[3];
495 double directionX[3];
496 double Grayon;
497 Grayon=tro.MajorRadius();
498 double Prayon;
499 Prayon=tro.MinorRadius();
500 gp_Ax3 repere=tro.Position();
501 gp_Dir dirX=repere.XDirection();
502 gp_Dir dirZ=repere.Direction();
503 directionX[0]=dirX.X();
504 directionX[1]=dirX.Y();
505 directionX[2]=dirX.Z();
506 direction[0]=dirZ.X();
507 direction[1]=dirZ.Y();
508 direction[2]=dirZ.Z();
509 param.ajouter(origine[0]);
510 param.ajouter(origine[1]);
511 param.ajouter(origine[2]);
512 param.ajouter(directionX[0]);
513 param.ajouter(directionX[1]);
514 param.ajouter(directionX[2]);
515 param.ajouter(direction[0]);
516 param.ajouter(direction[1]);
517 param.ajouter(direction[2]);
518 param.ajouter(Grayon);
519 param.ajouter(Prayon);
520
521 return MGCo_TORE;
522
523 }
524 //*******BSpline
525 if (type==STANDARD_TYPE(Geom_BSplineSurface))
526 {
527 Handle(Geom_BSplineSurface) bspline=Handle(Geom_BSplineSurface)::DownCast(surface);
528
529 //nombre des noeuds suivant Udirection
530 int nb_Uknot=bspline->NbUKnots();
531 //valeur de Unoeud
532 for (int i=1; i<=nb_Uknot; i++)
533 {
534 double Uvaleur=bspline->UKnot(i);
535 param.ajouter(Uvaleur);
536 }
537 //nombre des noeuds suivants Vdirection
538 int nb_Vknot=bspline->NbVKnots();
539 //valeur de Vnoeud
540 for (int j=1; j<=nb_Vknot; j++)
541 {
542 double Vvaleur=bspline->VKnot(j);
543 param.ajouter(Vvaleur);
544 }
545 //point de controle et poids
546 gp_Pnt pctr;
547 double poids;
548 for (int u=1; u<=bspline->NbUPoles(); u++)
549 for (int v=1; v<=bspline->NbVPoles(); v++)
550 {
551 pctr=bspline->Pole(u, v);
552
553 param.ajouter(pctr.X());
554 param.ajouter(pctr.Y());
555 param.ajouter(pctr.Z());
556
557 poids=bspline->Weight(u, v);
558 param.ajouter(poids);
559 }
560
561 double uDegree=bspline->UDegree();
562 param.ajouter(uDegree);
563 double vDegree=bspline->VDegree();
564 param.ajouter(vDegree);
565 return MGCo_BSPLINES;
566 }
567 return 0;
568 }
569 void OCC_SURFACE::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
570 {
571 //Conversion of the complete geometry of a shape into
572 //NURBS geometry
573 BRepBuilderAPI_NurbsConvert NURBS(face);
574 Handle(Geom_Surface) surface=BRepLib_FindSurface(NURBS).Surface();
575 Handle(Geom_BSplineSurface) bspline=GeomConvert::SurfaceToBSplineSurface(surface) ;
576
577 // The first parameter indicate the code access
578 param.ajouter(2);
579 //The follewing two parameters of the list indicate the orders of the net points
580 param.ajouter( bspline->UDegree()+1);
581 param.ajouter(bspline->VDegree()+1);
582
583 //The follewing two parameters indicate the number of rows and colons of the control points
584 //respectively to the two parameters directions
585 param.ajouter(bspline->NbUPoles());
586 param.ajouter(bspline->NbVPoles());
587
588 // this present the knot vector in the u-direction
589 for (unsigned int i=1;i<=bspline->NbUKnots();i++)
590 {
591 param.ajouter(bspline->UKnot(i));
592 }
593 //This present the knot vector in the v-direction
594 for (unsigned int j=1;j<=bspline->NbVKnots();j++)
595 {
596 param.ajouter(bspline->VKnot(j));
597 }
598 for (int v=1;v<=bspline->NbVPoles();v++)
599 {
600 for (int u=1;u<=bspline->NbUPoles();u++)
601 {
602 double w=bspline->Weight(u,v);
603 gp_Pnt point=bspline->Pole(u, v);
604 double x=point.X();
605 double y=point.Y();
606 double z=point.Z();
607 param.ajouter(x);
608 param.ajouter(y);
609 param.ajouter(z);
610 param.ajouter(w);
611 }
612
613 }
614 indx_premier_ptctr=5+bspline->NbUKnots()+bspline->NbVKnots();
615
616
617
618
619
620 }
621
622 void OCC_SURFACE::get_triangulation(MG_MAILLAGE* mai,MG_FACE* mgface,std::multimap<double,MG_NOEUD*,std::less<double> >& tabnoeudfus,double eps,int mode)
623 {
624 TPL_MAP_ENTITE<MG_SOMMET*> listsom;
625 TPL_MAP_ENTITE<MG_ARETE*> listare;
626 //std::map<unsigned long,bool> aretemaille;
627 std::map<unsigned long,std::map<double,NOEUDARETE,std::less<double> >,std::less<unsigned long> > areteamaille;
628 int nbboucle=mgface->get_nb_mg_boucle();
629 for (int i=0;i<nbboucle;i++)
630 {
631 MG_BOUCLE* bou=mgface->get_mg_boucle(i);
632 int nbarete=bou->get_nb_mg_coarete();
633 for (int j=0;j<nbarete;j++)
634 {
635 MG_ARETE* are=bou->get_mg_coarete(j)->get_arete();
636 listare.ajouter(are);
637 //bool amailler=true;
638 //if (are->get_lien_maillage()->get_nb()>0) amailler=false;
639 //aretemaille[are->get_id()]=amailler;
640 listsom.ajouter(bou->get_mg_coarete(j)->get_arete()->get_cosommet1()->get_sommet());
641 listsom.ajouter(bou->get_mg_coarete(j)->get_arete()->get_cosommet2()->get_sommet());
642 std::map<double,NOEUDARETE,std::less<double> > tmp;
643 std::pair<unsigned long,std::map<double,NOEUDARETE,std::less<double> > > maptmp(are->get_id(),tmp);
644 areteamaille.insert(maptmp);
645 }
646 }
647
648
649 TopLoc_Location L;
650 Handle (Poly_Triangulation) pt=BRep_Tool::Triangulation(face,L);
651 int nbnoeud=pt->NbNodes();
652 int nbmaille=pt->NbTriangles();
653 const TColgp_Array1OfPnt& nodes = pt->Nodes();
654 const Poly_Array1OfTriangle& triangles = pt->Triangles();
655 const TColgp_Array1OfPnt2d& uvNodes = pt->UVNodes();
656 std::vector<MG_NOEUD*> tabnoeud;
657 for ( Standard_Integer i = 0; i < nbnoeud; i++ )
658 {
659 gp_Pnt p1=nodes(i+1);
660 double xx=p1.X();
661 double yy=p1.Y();
662 double zz=p1.Z();
663 double key=fabs(xx)+fabs(yy)+fabs(zz);
664 MG_NOEUD* nvnoeud=NULL;
665 if (mode>1)
666 {
667 std::multimap<double,MG_NOEUD*,std::less<double> >::iterator it,itbas,ithaut;
668 itbas=tabnoeudfus.lower_bound(key*0.99);
669 ithaut=tabnoeudfus.upper_bound(key*1.1010101);
670 for ( it=itbas ; it != ithaut; it++ )
671 {
672 MG_NOEUD* ntmp=(*it).second;
673 double xtmp=ntmp->get_x();
674 double ytmp=ntmp->get_y();
675 double ztmp=ntmp->get_z();
676 OT_VECTEUR_3D vec(xtmp-xx,ytmp-yy,ztmp-zz);
677 if (vec.get_longueur()<1e-6*eps) {
678 nvnoeud=ntmp;
679 break;
680 }
681 }
682 }
683 if (nvnoeud==NULL)
684 {
685 MG_ELEMENT_TOPOLOGIQUE *topo=mgface;
686 TPL_MAP_ENTITE<MG_SOMMET*>::ITERATEUR it1;
687 for (MG_SOMMET* som=listsom.get_premier(it1);som!=NULL;som=listsom.get_suivant(it1))
688 {
689 double xyz[3];
690 som->get_point()->evaluer(xyz);
691 OT_VECTEUR_3D vec(xyz[0]-xx,xyz[1]-yy,xyz[2]-zz);
692 if (vec.get_longueur()<1e-6*eps)
693 {
694 topo=som;
695 break;
696
697 }
698 }
699 double param_t;
700 if (topo==mgface)
701 {
702 TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR it2;
703 for (MG_ARETE* are=listare.get_premier(it2);are!=NULL;are=listare.get_suivant(it2))
704 {
705 double t;
706 double xyz[3]={xx,yy,zz};
707 are->inverser(t,xyz);
708 if (are->get_courbe()->est_periodique())
709 if (t< are->get_tmin()) t=t+are->get_courbe()->get_periode();
710 double xyztmp[3];
711 are->evaluer(t,xyztmp );
712 OT_VECTEUR_3D vec(xyz,xyztmp);
713 if (vec.get_longueur()<1e-6*eps)
714 if ((t>are->get_tmin()) && (t<are->get_tmax()))
715 {
716 topo=are;
717 param_t=t;
718 break;
719 }
720 }
721 }
722 nvnoeud=new MG_NOEUD(topo,xx,yy,zz,MAGIC::ORIGINE::TRIANGULATION);
723 mai->ajouter_mg_noeud(nvnoeud);
724 std::pair<double,MG_NOEUD*> tmp(key,nvnoeud);
725 tabnoeudfus.insert(tmp);
726 if (topo->get_dimension()==1)
727 {
728 NOEUDARETE na;
729 na.no=nvnoeud;
730 na.are=(MG_ARETE*)topo;
731 na.t=param_t;
732 std::pair<double,NOEUDARETE> tmp(na.t,na);
733 areteamaille[topo->get_id()].insert(tmp);
734 }
735 }
736
737
738 tabnoeud.insert(tabnoeud.end(),nvnoeud);
739 }
740 for ( Standard_Integer i = 0; i < nbmaille; i++ )
741 {
742 int n1,n2,n3;
743 Poly_Triangle triangle = triangles( i + 1 );
744 bool face_reversed = (face.Orientation() == TopAbs_REVERSED);
745 if ( face_reversed )
746 triangle.Get( n1, n3, n2 );
747 else
748 triangle.Get( n1, n2, n3 );
749 MG_NOEUD* noeud1=tabnoeud[n1-1];
750 MG_NOEUD* noeud2=tabnoeud[n2-1];
751 MG_NOEUD* noeud3=tabnoeud[n3-1];
752 if (noeud1==noeud2) continue;
753 if (noeud1==noeud3) continue;
754 if (noeud2==noeud3) continue;
755 mai->ajouter_mg_triangle(mgface,noeud1,noeud2,noeud3,MAGIC::ORIGINE::TRIANGULATION);
756 /*if (noeud1->get_lien_topologie()->get_dimension()==0)
757 {
758 MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();
759 som->get_lien_maillage()->ajouter(noeud1);
760 }
761 if (noeud2->get_lien_topologie()->get_dimension()==0)
762 {
763 MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();
764 som->get_lien_maillage()->ajouter(noeud2);
765 }
766 if (noeud3->get_lien_topologie()->get_dimension()==0)
767 {
768 MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();
769 som->get_lien_maillage()->ajouter(noeud3);
770 }*/
771 /*if (noeud1->get_lien_topologie()==noeud2->get_lien_topologie())
772 if (noeud1->get_lien_topologie()->get_dimension()==1)
773 {
774 MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
775 if (aretemaille[are->get_id()]==true)
776 {
777 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud2,TRIANGULATION);
778 are->get_lien_maillage()->ajouter(seg);
779 }
780 }
781 if (noeud1->get_lien_topologie()==noeud3->get_lien_topologie())
782 if (noeud1->get_lien_topologie()->get_dimension()==1)
783 {
784 MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
785 if (aretemaille[are->get_id()]==true)
786 {
787 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud3,TRIANGULATION);
788 are->get_lien_maillage()->ajouter(seg);
789 }
790 }#include <../gmsh/tutorial/t8.geo>
791 if (noeud3->get_lien_topologie()==noeud2->get_lien_topologie())
792 if (noeud3->get_lien_topologie()->get_dimension()==1)
793 {
794 MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
795 if (aretemaille[are->get_id()]==true)
796 {
797 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud2,TRIANGULATION);
798 are->get_lien_maillage()->ajouter(seg);
799 }
800 }
801 if (noeud1->get_lien_topologie()->get_dimension()==1)
802 if (noeud2->get_lien_topologie()->get_dimension()==0)
803 {
804 MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
805 MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();
806 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
807 {
808 if (aretemaille[are->get_id()]==true)
809 {
810 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud2,TRIANGULATION);
811 are->get_lien_maillage()->ajouter(seg);
812 }
813 }
814 }
815 if (noeud1->get_lien_topologie()->get_dimension()==1)
816 if (noeud3->get_lien_topologie()->get_dimension()==0)
817 {
818 MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
819 MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();
820 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
821 {
822 if (aretemaille[are->get_id()]==true)
823 {
824 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud3,TRIANGULATION);
825 are->get_lien_maillage()->ajouter(seg);
826 }
827 }
828 }
829 if (noeud2->get_lien_topologie()->get_dimension()==1)
830 if (noeud1->get_lien_topologie()->get_dimension()==0)
831 {
832 MG_ARETE* are=(MG_ARETE*)noeud2->get_lien_topologie();
833 MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();
834 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
835 {
836 if (aretemaille[are->get_id()]==true)
837 {
838 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud2,noeud1,TRIANGULATION);
839 are->get_lien_maillage()->ajouter(seg);
840 }
841 }
842 }
843 if (noeud2->get_lien_topologie()->get_dimension()==1)
844 if (noeud3->get_lien_topologie()->get_dimension()==0)
845 {
846 MG_ARETE* are=(MG_ARETE*)noeud2->get_lien_topologie();
847 MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();
848 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
849 {
850 if (aretemaille[are->get_id()]==true)
851 {
852 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud2,noeud3,TRIANGULATION);
853 are->get_lien_maillage()->ajouter(seg);
854 }
855 }
856 }
857 if (noeud3->get_lien_topologie()->get_dimension()==1)
858 if (noeud1->get_lien_topologie()->get_dimension()==0)
859 {
860 MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
861 MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();
862 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
863 {
864 if (aretemaille[are->get_id()]==true)
865 {
866 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud1,TRIANGULATION);
867 are->get_lien_maillage()->ajouter(seg);
868 }
869 }
870 }
871 if (noeud3->get_lien_topologie()->get_dimension()==1)
872 if (noeud2->get_lien_topologie()->get_dimension()==0)
873 {
874 MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
875 MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();
876 if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
877 {
878 if (aretemaille[are->get_id()]==true)
879 {
880 MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud2,TRIANGULATION);
881 are->get_lien_maillage()->ajouter(seg);
882 }
883 }
884 }*/
885 }
886 if (mode>1)
887 {
888 TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR it2;
889 for (MG_ARETE* are=listare.get_premier(it2);are!=NULL;are=listare.get_suivant(it2))
890 {
891 if (are->get_lien_maillage()->get_nb()==0)
892 {
893 unsigned long id=are->get_id();
894 MG_NOEUD* nodep=(MG_NOEUD*)are->get_cosommet1()->get_sommet()->get_lien_maillage()->get(0);
895 MG_NOEUD* noarr=(MG_NOEUD*)are->get_cosommet2()->get_sommet()->get_lien_maillage()->get(0);
896 std::map<double,NOEUDARETE,std::less<double> >::iterator it=areteamaille[id].begin();
897 MG_NOEUD* noeudcourant=nodep;
898 while (it!=areteamaille[id].end())
899 {
900 MG_NOEUD* noeud=(*it).second.no;
901 //MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeudcourant,noeud,TRIANGULATION);
902 MG_SEGMENT* seg=mai->get_mg_segment(noeudcourant->get_id(),noeud->get_id());
903 seg->change_lien_topologie(are);
904 noeudcourant=noeud;
905 it++;
906 }
907 //MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeudcourant,noarr,TRIANGULATION);
908 MG_SEGMENT* seg=mai->get_mg_segment(noeudcourant->get_id(),noarr->get_id());
909 seg->change_lien_topologie(are);
910 }
911 }
912
913 }
914 }
915
916 void OCC_SURFACE::get_liste_pole(std::vector< double> *liste_pole,double eps)
917 {
918 //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
919 Handle(Standard_Type) type=surface->DynamicType();
920 //std::cout << surface->DynamicType()->Name() <<std::endl;
921 if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
922 {
923 Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
924 Handle(Geom_Surface) basissurface = RTSurface->BasisSurface();
925 Handle(Standard_Type) type_basissurface=basissurface->DynamicType();
926 if(type_basissurface==STANDARD_TYPE(Geom_SphericalSurface))
927 {
928 type=STANDARD_TYPE(Geom_SphericalSurface);
929 }
930 }
931 if(type==STANDARD_TYPE(Geom_SphericalSurface))
932 {
933 BRepClass_FaceClassifier faceclassifier;
934 gp_Pnt2d pnt2d_pole_sud(0.0,-M_PI/2.);
935 faceclassifier.Perform(face,pnt2d_pole_sud,eps);
936 if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
937 {
938 liste_pole->push_back(0.0);
939 liste_pole->push_back(-M_PI/2.);
940 }
941 gp_Pnt2d pnt2d_pole_nord(0.0,M_PI/2.);
942 faceclassifier.Perform(face,pnt2d_pole_nord,eps);
943 if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
944 {
945 liste_pole->push_back(0.0);
946 liste_pole->push_back(M_PI/2.);
947 }
948 }
949 if(type==STANDARD_TYPE(Geom_SurfaceOfRevolution))
950 {
951 BRepClass_FaceClassifier faceclassifier;
952 gp_Pnt2d pnt2d_pole_sud(0.0,v_min);
953 faceclassifier.Perform(face,pnt2d_pole_sud,eps);
954 if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
955 {
956 liste_pole->push_back(0.0);
957 liste_pole->push_back(v_min);
958 }
959 gp_Pnt2d pnt2d_pole_nord(0.0,v_max);
960 faceclassifier.Perform(face,pnt2d_pole_nord,eps);
961 if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
962 {
963 liste_pole->push_back(0.0);
964 liste_pole->push_back(v_max);
965 }
966 }
967 // if(type==STANDARD_TYPE(Geom_ConicalSurface))
968 // {
969 // BRepClass_FaceClassifier faceclassifier;
970 // gp_Pnt2d pnt2d_pole(0.0,v_max);
971 // faceclassifier.Perform(face,pnt2d_pole,eps);
972 // if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
973 // {
974 // liste_pole->push_back(0.0);
975 // liste_pole->push_back(v_max);
976 // }
977 // }
978 }
979
980
981 void OCC_SURFACE::analyse_bspline(void)
982 {
983 Handle(Geom_BSplineSurface) bspline=Handle(Geom_BSplineSurface)::DownCast(surface);
984 //point de controle et poids
985 bool ferme=true;
986 for (int v=1; v<=bspline->NbVPoles(); v++)
987 {
988 gp_Pnt pctr1=bspline->Pole(1, v);
989 gp_Pnt pctr2=bspline->Pole(bspline->NbUPoles(), v);
990 double poids1=bspline->Weight(1, v);
991 double poids2=bspline->Weight(bspline->NbUPoles(), v);
992 double xyz1[4];
993 double xyz2[4];
994 xyz1[0]=pctr1.X();
995 xyz1[1]=pctr1.Y();
996 xyz1[2]=pctr1.Z();
997 xyz1[3]=poids1;
998 xyz2[0]=pctr2.X();
999 xyz2[1]=pctr2.Y();
1000 xyz2[2]=pctr2.Z();
1001 xyz2[3]=poids2;
1002 OT_VECTEUR_4D vec(xyz1,xyz2);
1003 if (vec.get_longueur()>fonction1->get_precision())
1004 ferme=false;
1005
1006 }
1007 if (ferme)
1008 {
1009 estperiodeu=1;
1010 periode_u=u_max-u_min;
1011 }
1012 ferme=true;
1013 for (int u=1; u<=bspline->NbUPoles(); u++)
1014 {
1015 gp_Pnt pctr1=bspline->Pole(u,1);
1016 gp_Pnt pctr2=bspline->Pole(u,bspline->NbVPoles());
1017 double poids1=bspline->Weight(u,1);
1018 double poids2=bspline->Weight(u,bspline->NbVPoles());
1019 double xyz1[4];
1020 double xyz2[4];
1021 xyz1[0]=pctr1.X();
1022 xyz1[1]=pctr1.Y();
1023 xyz1[2]=pctr1.Z();
1024 xyz1[3]=poids1;
1025 xyz2[0]=pctr2.X();
1026 xyz2[1]=pctr2.Y();
1027 xyz2[2]=pctr2.Z();
1028 xyz2[3]=poids2;
1029 OT_VECTEUR_4D vec(xyz1,xyz2);
1030 if (vec.get_longueur()>fonction1->get_precision())
1031 ferme=false;
1032
1033 }
1034 if (ferme)
1035 {
1036 estperiodev=1;
1037 periode_v=v_max-v_min;
1038 }
1039 }
1040
1041 void OCC_SURFACE::inverser2(double *uv,double *xyz,double precision)
1042 {
1043 int NUMDECOUP=20;
1044 double u1;
1045 double u2;
1046 double v1;
1047 double v2;
1048 double u,v,distmin=1e300;
1049 surface->Bounds(u1,u2,v1, v2);
1050 for (int i=0;i<NUMDECOUP+1;i++)
1051 for (int j=0;j<NUMDECOUP+1;j++)
1052 {
1053 double ut=u1+i*1.0/NUMDECOUP*(u2-u1);
1054 double vt=v1+j*1.0/NUMDECOUP*(v2-v1);
1055 gp_Pnt P;
1056 surface->D0(ut,vt,P);
1057 double xyz2[3]={P.X(),P.Y(),P.Z()};
1058 OT_VECTEUR_3D vec(xyz,xyz2);
1059 double dist=vec.get_longueur();
1060 if (dist<distmin)
1061 {
1062 distmin=dist;
1063 u=ut;
1064 v=vt;
1065 }
1066 }
1067 bool sortie=false;
1068 int nbiteration=0;
1069 while (sortie==false)
1070 {
1071 double uv2[3]={u,v,0.};
1072 double xyzuv[3],xyzuvdu[3],xyzuvdv[3],xyzuvduu[3],xyzuvduv[3],xyzuvdvv[3];
1073 deriver_seconde(uv2,xyzuvduu,xyzuvduv,xyzuvdvv,xyzuv,xyzuvdu,xyzuvdv);
1074 OT_VECTEUR_3D vxyz(xyz);
1075 OT_VECTEUR_3D vxyzuv(xyzuv);
1076 OT_VECTEUR_3D P=vxyzuv-vxyz;
1077 OT_VECTEUR_3D vxyzuvdu(xyzuvdu);
1078 OT_VECTEUR_3D vxyzuvdv(xyzuvdv);
1079 OT_VECTEUR_3D vxyzuvduu(xyzuvduu);
1080 OT_VECTEUR_3D vxyzuvdvv(xyzuvdvv);
1081 OT_VECTEUR_3D vxyzuvduv(xyzuvduv);
1082 double Fdu=P*vxyzuvduu+vxyzuvdu*vxyzuvdu;
1083 double Fdv=P*vxyzuvduv+vxyzuvdv*vxyzuvdu;
1084 double Gdu=P*vxyzuvduv+vxyzuvdu*vxyzuvdv;
1085 double Gdv=P*vxyzuvdvv+vxyzuvdv*vxyzuvdv;
1086 double F=P*xyzuvdu;
1087 double G=P*xyzuvdv;
1088 double det=Fdu*Gdv-Gdu*Fdv;
1089 double deltau=(-F*Gdv+Gdu*G)/det;
1090 double deltav=(-Fdu*G+F*Fdv)/det;
1091 double new_u=u+deltau;
1092 double new_v=v+deltav;
1093 if (new_u<u1) new_u=u1;
1094 if (new_u>u2) new_u=u2;
1095 if (new_v<v1) new_v=v1;
1096 if (new_v>v2) new_v=v2;
1097 if (fabs(new_u-u)<precision)
1098 if (fabs(new_v-v)<precision)
1099 sortie=true;
1100 u=new_u;
1101 v=new_v;
1102 nbiteration++;
1103 if (nbiteration>1000) sortie=true;
1104 }
1105 uv[0]=u;
1106 uv[1]=v;
1107 }
1108
1109 #endif