ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/occ_surface.cpp
Revision: 339
Committed: Wed May 30 18:20:47 2012 UTC (12 years, 11 months ago) by francois
File size: 15402 byte(s)
Log Message:
Correction bug mailleur bloc + correction bug inversion avec open cascade + preparation pour element  XFEM

File Contents

# User Rev Content
1 francois 283 //---------------------------------------------------------------------------
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 <gp_Sphere.hxx>
47     #include "Geom_ToroidalSurface.hxx"
48     #include <gp_Torus.hxx>
49     #include <Geom_BSplineSurface.hxx>
50     #include <Geom_BezierSurface.hxx>
51     #include <GeomConvert.hxx>
52     #include <BRepBuilderAPI_NurbsConvert.hxx>
53     #include <BRepLib_FindSurface.hxx>
54     #include "mg_gestionnaire.h"
55     #include "constantegeo.h"
56 francois 295 #include "ot_mathematique.h"
57 francois 283
58     #pragma package(smart_init)
59     OCC_SURFACE::OCC_SURFACE(unsigned long num, TopoDS_Face srf, OCC_FONCTION1& fonc):MG_SURFACE(num),face(srf), fonction1(fonc)
60     {
61    
62     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
63     double u1;
64     double u2;
65     double v1;
66     double v2;
67     surface->Bounds(u1,u2,v1, v2);
68     u_min=u1;
69     u_max=u2;
70     v_min=v1;
71     v_max=v2;
72    
73     }
74    
75     OCC_SURFACE::OCC_SURFACE(TopoDS_Face srf, OCC_FONCTION1& fonc):MG_SURFACE(),face(srf), fonction1(fonc)
76     {
77     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
78     double u1;
79     double u2;
80     double v1;
81     double v2;
82     surface->Bounds(u1,u2,v1, v2);
83     u_min=u1;
84     u_max=u2;
85     v_min=v1;
86     v_max=v2;
87     }
88    
89     OCC_SURFACE::OCC_SURFACE(OCC_SURFACE& mdd):MG_SURFACE(mdd),face(mdd.face), fonction1(mdd.fonction1)
90     {
91     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
92     u_min=mdd.get_umin();
93     u_max=mdd.get_umax();
94     v_min=mdd.get_vmin();
95     v_max=mdd.get_vmax();
96    
97     }
98     OCC_SURFACE::~OCC_SURFACE()
99     {
100     }
101     void OCC_SURFACE::evaluer(double *uv,double *xyz)
102     {
103    
104     const Handle(Geom_Surface) &surface=BRep_Tool::Surface(face);
105     gp_Pnt P;
106     double u=uv[0];
107     double v=uv[1];
108    
109     surface->D0(u,v,P);
110    
111     xyz[0]=P.X();
112     xyz[1]=P.Y();
113     xyz[2]=P.Z();
114     }
115     void OCC_SURFACE::deriver(double *uv,double *xyzdu, double *xyzdv)
116     {
117     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
118     double u=uv[0];
119     double v=uv[1];
120     gp_Vec D1U;
121     gp_Vec D1V;
122     gp_Pnt P;
123    
124     surface->D1(u,v,P,D1U,D1V);
125     xyzdu[0]=D1U.X();
126     xyzdu[1]=D1U.Y();
127     xyzdu[2]=D1U.Z();
128    
129     xyzdv[0]=D1V.X();
130     xyzdv[1]=D1V.Y();
131     xyzdv[2]=D1V.Z();
132    
133     }
134     void OCC_SURFACE::deriver_seconde(double *uv,double* xyzduu,double* xyzduv,double* xyzdvv,double *xyz, double *xyzdu, double *xyzdv)
135     {
136     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
137     double u=uv[0];
138     double v=uv[1];
139     gp_Pnt P;
140     gp_Vec D1U;
141     gp_Vec D1V;
142     gp_Vec D2U;
143     gp_Vec D2V;
144     gp_Vec D2UV;
145    
146     surface->D2(u,v,P,D1U,D1V,D2U,D2V,D2UV);
147    
148     xyz[0]=P.X();
149     xyz[1]=P.Y();
150     xyz[2]=P.Z();
151    
152     xyzdu[0]=D1U.X();
153     xyzdu[1]=D1U.Y();
154     xyzdu[2]=D1U.Z();
155    
156     xyzdv[0]=D1V.X();
157     xyzdv[1]=D1V.Y();
158     xyzdv[2]=D1V.Z();
159    
160     xyzduu[0]=D2U.X();
161     xyzduu[1]=D2U.Y();
162     xyzduu[2]=D2U.Z();
163    
164     xyzdvv[0]=D2V.X();
165     xyzdvv[1]=D2V.Y();
166     xyzdvv[2]=D2V.Z();
167    
168     xyzduv[0]=D2UV.X();
169     xyzduv[1]=D2UV.Y();
170     xyzduv[2]=D2UV.Z();
171    
172     }
173    
174     void OCC_SURFACE::inverser(double *uv,double *xyz,double precision)
175     {
176     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
177     double u=xyz[0];
178     double v=xyz[1];
179     double w=xyz[2];
180     gp_Pnt P(u,v,w);
181    
182     //ShapeAnalysis_Surface SAS(surface);
183     //gp_Pnt2d pnt2d=SAS.ValueOfUV(P, precision);
184     GeomAPI_ProjectPointOnSurf PPS(P,surface, precision);
185 francois 339 if (PPS.NbPoints() < 1)
186     {
187     uv[0]=1e308;
188     uv[1]=1e308;
189     return;
190     }
191 francois 283 //PPS.Perform(P);
192     double UU, VV;
193     PPS.LowerDistanceParameters(UU,VV);
194     uv[0]=UU;
195     uv[1]=VV;
196     }
197     int OCC_SURFACE::est_periodique_u(void)
198     {
199     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
200     return surface->IsUClosed();
201     }
202    
203     int OCC_SURFACE::est_periodique_v(void)
204     {
205     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
206     return surface->IsVPeriodic();
207     }
208     double OCC_SURFACE::get_periode_u(void)
209     {
210     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
211     if (surface->IsUPeriodic()==0) return 0.;
212     return surface->UPeriod();
213     }
214    
215     double OCC_SURFACE::get_periode_v(void)
216     {
217     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
218     if (surface->IsVPeriodic()==0) return 0.;
219     return surface->VPeriod();
220     }
221     void OCC_SURFACE::enregistrer(std::ostream& o)
222     {
223     o <<"%"<<get_id()<< "=SURFACE_OCC("<< fonction1.GetID(face)<< ");" << std::endl;
224     }
225     int OCC_SURFACE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
226     {
227     Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
228     Handle(Standard_Type) type=surface->DynamicType();
229     //******plan
230     if (type==STANDARD_TYPE(Geom_Plane))
231     {
232     Handle(Geom_Plane) Pln=Handle(Geom_Plane)::DownCast(surface);
233     gp_Pln plan=Pln->Pln();
234    
235     double origine[3];
236     gp_Pnt centre=plan.Location();
237    
238     origine[0]=centre.X();
239     origine[1]=centre.Y();
240     origine[2]=centre.Z();
241    
242     double normal[3];
243     gp_Ax1 axe=plan.Axis();
244     gp_Dir direction=axe.Direction();
245    
246     normal[0]=direction.X();
247     normal[1]=direction.Y();
248     normal[2]=direction.Z();
249    
250     param.ajouter(origine[0]);
251     param.ajouter(origine[1]);
252     param.ajouter(origine[2]);
253     param.ajouter(normal[0]);
254     param.ajouter(normal[1]);
255     param.ajouter(normal[2]);
256    
257     return MGCo_PLAN;
258     }
259     //******cylindre
260     if (type==STANDARD_TYPE(Geom_CylindricalSurface))
261     {
262     Handle(Geom_CylindricalSurface) cylinder=Handle(Geom_CylindricalSurface)::DownCast(surface);
263     gp_Cylinder cylin=cylinder->Cylinder();
264    
265     double origine[3];
266     gp_Pnt centre=cylin.Location();
267     origine[0]=centre.X();
268     origine[1]=centre.Y();
269     origine[2]=centre.Z();
270     double direction[3];
271     gp_Ax1 axe=cylin.Axis();
272     gp_Dir dir=axe.Direction();
273     direction[0]=dir.X();
274     direction[1]=dir.Y();
275     direction[2]=dir.Z();
276     double rayon;
277     rayon=cylin.Radius();
278    
279     param.ajouter(origine[0]);
280     param.ajouter(origine[1]);
281     param.ajouter(origine[2]);
282     param.ajouter(direction[0]);
283     param.ajouter(direction[1]);
284     param.ajouter(direction[2]);
285     param.ajouter(rayon);
286    
287     return MGCo_CYLINDRE;
288     }
289     //******Cone
290     if (type==STANDARD_TYPE(Geom_ConicalSurface))
291     {
292     Handle(Geom_ConicalSurface) cone=Handle(Geom_ConicalSurface)::DownCast(surface);
293     gp_Cone con=cone->Cone();
294    
295     double origine[3];
296     gp_Pnt centre=con.Location();
297     origine[0]=centre.X();
298     origine[1]=centre.Y();
299     origine[2]=centre.Z();
300     double direction[3];
301     gp_Ax1 axe=con.Axis();
302     gp_Dir dir=axe.Direction();
303     direction[0]=dir.X();
304     direction[1]=dir.Y();
305     direction[2]=dir.Z();
306     double rayon;
307     rayon=con.RefRadius();
308     double angle;
309     angle=con.SemiAngle();
310    
311     param.ajouter(origine[0]);
312     param.ajouter(origine[1]);
313     param.ajouter(origine[2]);
314     param.ajouter(direction[0]);
315     param.ajouter(direction[1]);
316     param.ajouter(direction[2]);
317     param.ajouter(rayon);
318     param.ajouter(angle);
319    
320     return MGCo_CONE;
321    
322     }
323     //*****Sphere
324     if (type==STANDARD_TYPE(Geom_SphericalSurface))
325     {
326     Handle(Geom_SphericalSurface) sphere=Handle(Geom_SphericalSurface)::DownCast(surface);
327     gp_Sphere sph=sphere->Sphere();
328    
329     double origine[3];
330     gp_Pnt centre=sph.Location();
331     origine[0]=centre.X();
332     origine[1]=centre.Y();
333     origine[2]=centre.Z();
334     double rayon;
335     rayon=sph.Radius();
336     param.ajouter(origine[0]);
337     param.ajouter(origine[1]);
338     param.ajouter(origine[2]);
339     param.ajouter(rayon);
340    
341     return MGCo_SPHERE;
342    
343     }
344     //****** Tore
345     if (type==STANDARD_TYPE(Geom_ToroidalSurface))
346     {
347     Handle(Geom_ToroidalSurface) tore=Handle(Geom_ToroidalSurface)::DownCast(surface);
348     gp_Torus tro=tore->Torus();
349     double origine[3];
350     gp_Pnt centre=tro.Location();
351     origine[0]=centre.X();
352     origine[1]=centre.Y();
353     origine[2]=centre.Z();
354     double direction[3];
355     gp_Ax1 axe=tro.Axis();
356     gp_Dir dir=axe.Direction();
357     direction[0]=dir.X();
358     direction[1]=dir.Y();
359     direction[2]=dir.Z();
360     double Grayon;
361     Grayon=tro.MajorRadius();
362     double Prayon;
363     Prayon=tro.MinorRadius();
364    
365     param.ajouter(origine[0]);
366     param.ajouter(origine[1]);
367     param.ajouter(origine[2]);
368     param.ajouter(direction[0]);
369     param.ajouter(direction[1]);
370     param.ajouter(direction[2]);
371     param.ajouter(Grayon);
372     param.ajouter(Prayon);
373    
374     return MGCo_TORE;
375    
376     }
377     //*******BSpline
378     if (type==STANDARD_TYPE(Geom_BSplineSurface))
379     {
380     Handle(Geom_BSplineSurface) bspline=Handle(Geom_BSplineSurface)::DownCast(surface);
381    
382     //nombre des noeuds suivant Udirection
383     int nb_Uknot=bspline->NbUKnots();
384     //valeur de Unoeud
385     for (int i=1; i<=nb_Uknot; i++)
386     {
387     double Uvaleur=bspline->UKnot(i);
388     param.ajouter(Uvaleur);
389     }
390     //nombre des noeuds suivants Vdirection
391     int nb_Vknot=bspline->NbVKnots();
392     //valeur de Vnoeud
393     for (int j=1; j<=nb_Vknot; j++)
394     {
395     double Vvaleur=bspline->VKnot(j);
396     param.ajouter(Vvaleur);
397     }
398     //point de controle et poids
399     gp_Pnt pctr;
400     double poids;
401     for (int u=1; u<=bspline->NbUPoles(); u++)
402     for (int v=1; v<=bspline->NbVPoles(); v++)
403     {
404     pctr=bspline->Pole(u, v);
405    
406     param.ajouter(pctr.X());
407     param.ajouter(pctr.Y());
408     param.ajouter(pctr.Z());
409    
410     poids=bspline->Weight(u, v);
411     param.ajouter(poids);
412     }
413    
414     double uDegree=bspline->UDegree();
415     param.ajouter(uDegree);
416     double vDegree=bspline->VDegree();
417     param.ajouter(vDegree);
418     return MGCo_BSPLINES;
419     }
420    
421     }
422     void OCC_SURFACE::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
423     {
424     //Conversion of the complete geometry of a shape into
425     //NURBS geometry
426     BRepBuilderAPI_NurbsConvert NURBS(face);
427     Handle(Geom_Surface) surface=BRepLib_FindSurface(NURBS).Surface();
428     Handle(Geom_BSplineSurface) bspline=GeomConvert::SurfaceToBSplineSurface(surface) ;
429    
430     // The first parameter indicate the code access
431     param.ajouter(2);
432     //The follewing two parameters of the list indicate the orders of the net points
433     param.ajouter( bspline->UDegree()+1);
434     param.ajouter(bspline->VDegree()+1);
435    
436     //The follewing two parameters indicate the number of rows and colons of the control points
437     //respectively to the two parameters directions
438     param.ajouter(bspline->NbUPoles());
439     param.ajouter(bspline->NbVPoles());
440    
441     // this present the knot vector in the u-direction
442     for (unsigned int i=1;i<=bspline->NbUKnots();i++)
443     {
444     param.ajouter(bspline->UKnot(i));
445     }
446     //This present the knot vector in the v-direction
447     for (unsigned int j=1;j<=bspline->NbVKnots();j++)
448     {
449     param.ajouter(bspline->VKnot(j));
450     }
451    
452     for (int u=1;u<=bspline->NbUPoles();u++)
453     {
454     for (int v=1;v<=bspline->NbVPoles();v++)
455     {
456     double w=bspline->Weight(u,v);
457     gp_Pnt point=bspline->Pole(u, v);
458     double x=point.X();
459     double y=point.Y();
460     double z=point.Z();
461     param.ajouter(x);
462     param.ajouter(y);
463     param.ajouter(z);
464     param.ajouter(w);
465     }
466    
467     }
468     indx_premier_ptctr=5+bspline->NbUKnots()+bspline->NbVKnots();
469    
470    
471    
472    
473    
474     }
475    
476 francois 295 void OCC_SURFACE::get_triangulation(MG_MAILLAGE* mai,MG_FACE* mgface,std::multimap<double,MG_NOEUD*,std::less<double> >& tabnoeudfus,double eps,int mode)
477 francois 283 {
478     TopLoc_Location L;
479     Handle (Poly_Triangulation) pt=BRep_Tool::Triangulation(face,L);
480     int nbnoeud=pt->NbNodes();
481     int nbmaille=pt->NbTriangles();
482     const TColgp_Array1OfPnt& nodes = pt->Nodes();
483     const Poly_Array1OfTriangle& triangles = pt->Triangles();
484     const TColgp_Array1OfPnt2d& uvNodes = pt->UVNodes();
485 francois 295 std::vector<MG_NOEUD*> tabnoeud;
486 francois 283 for ( Standard_Integer i = 0; i < nbnoeud; i++ )
487     {
488     gp_Pnt p1=nodes(i+1);
489     double xx=p1.X();
490     double yy=p1.Y();
491     double zz=p1.Z();
492     double key=fabs(xx)+fabs(yy)+fabs(zz);
493     MG_NOEUD* nvnoeud=NULL;
494     if (mode>1)
495     {
496 francois 295 std::multimap<double,MG_NOEUD*,std::less<double> >::iterator it,itbas,ithaut;
497 francois 283 itbas=tabnoeudfus.lower_bound(key*0.99);
498     ithaut=tabnoeudfus.upper_bound(key*1.1010101);
499     for ( it=itbas ; it != ithaut; it++ )
500     {
501     MG_NOEUD* ntmp=(*it).second;
502     double xtmp=ntmp->get_x();
503     double ytmp=ntmp->get_y();
504     double ztmp=ntmp->get_z();
505     OT_VECTEUR_3D vec(xtmp-xx,ytmp-yy,ztmp-zz);
506     if (vec.get_longueur()<1e-6*eps) {
507     nvnoeud=ntmp;
508     break;
509     }
510     }
511     }
512     if (nvnoeud==NULL)
513     {
514     nvnoeud=new MG_NOEUD(mgface,xx,yy,zz,TRIANGULATION);
515     mai->ajouter_mg_noeud(nvnoeud);
516     std::pair<double,MG_NOEUD*> tmp(key,nvnoeud);
517     tabnoeudfus.insert(tmp);
518     }
519    
520    
521     tabnoeud.insert(tabnoeud.end(),nvnoeud);
522     }
523     for ( Standard_Integer i = 0; i < nbmaille; i++ )
524     {
525     int n1,n2,n3;
526     Poly_Triangle triangle = triangles( i + 1 );
527     bool face_reversed = (face.Orientation() == TopAbs_REVERSED);
528     if ( face_reversed )
529     triangle.Get( n1, n3, n2 );
530     else
531     triangle.Get( n1, n2, n3 );
532     MG_NOEUD* noeud1=tabnoeud[n1-1];
533     MG_NOEUD* noeud2=tabnoeud[n2-1];
534     MG_NOEUD* noeud3=tabnoeud[n3-1];
535     mai->ajouter_mg_triangle(mgface,noeud1,noeud2,noeud3,TRIANGULATION);
536    
537    
538     }
539     }
540    
541     #endif