mtu/src/occ_surface.cpp

//---------------------------------------------------------------------------
//------------------------------------------------------------
//------------------------------------------------------------
//  MAGiC
//  Jean Christophe Cuilli�re et Vincent FRANCOIS
//  D�partement de G�nie M�canique - UQTR
//------------------------------------------------------------
//  Le projet MAGIC est un projet de recherche du d�partement
//  de g�nie m�canique de l'Universit� du Qu�bec �
//  Trois Rivi�res
//  Les librairies ne peuvent �tre utilis�es sans l'accord
//  des auteurs (contact : francois@uqtr.ca)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       OCC_Surface.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 � 11H22
//------------------------------------------------------------
//------------------------------------------------------------

#pragma hdrstop
#include "gestionversion.h"
#ifdef BREP_OCC

#include "occ_surface.h"
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include "GeomAPI_ProjectPointOnSurf.hxx"
#include <ShapeAnalysis_Surface.hxx>
#include <BRep_Tool.hxx>
#include <BRepAdaptor_Surface.hxx>
//***********************************
#include <BRep_Tool.hxx>
#include <Poly_Triangulation.hxx>
#include <Geom_Plane.hxx>
#include <gp_Pln.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <gp_Cylinder.hxx>
#include <Geom_ConicalSurface.hxx>
#include <gp_Cone.hxx>
#include <Geom_SphericalSurface.hxx>
#include <IGESSolid_SphericalSurface.hxx>
#include <gp_Sphere.hxx>
#include "Geom_ToroidalSurface.hxx"
#include <gp_Torus.hxx>
#include <Geom_BSplineSurface.hxx>
#include <Geom_BezierSurface.hxx>
#include <GeomConvert.hxx>
#include <BRepBuilderAPI_NurbsConvert.hxx>
#include <BRepLib_FindSurface.hxx>
#include "mg_gestionnaire.h"
#include "constantegeo.h"
#include "ot_mathematique.h"
#include <Geom_RectangularTrimmedSurface.hxx>
#include <BRepClass_FaceClassifier.hxx>
#include <gp_Pnt2d.hxx>
#include <Geom_SurfaceOfRevolution.hxx>

#pragma package(smart_init)


class NOEUDARETE
{
public:
  MG_NOEUD* no;
  MG_ARETE* are;
  double t;
};


OCC_SURFACE::OCC_SURFACE(unsigned long num, TopoDS_Face srf, OCC_FONCTION* fonc):MG_SURFACE(num),face(srf), fonction1(fonc)
{
    surface=BRep_Tool::Surface(face);
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    Handle(Standard_Type) type=surface->DynamicType(); 
    if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
        {
        Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
        surface = RTSurface->BasisSurface();
        }
    double u1;
    double u2;
    double v1;
    double v2;
    surface->Bounds(u1,u2,v1, v2);
    u_min=u1;
    u_max=u2;
    v_min=v1;
    v_max=v2;
    estperiodeu=surface->IsUClosed();
    estperiodev=surface->IsVClosed();
    if (surface->IsUPeriodic()) periode_u=surface->UPeriod();
    else if (surface->IsUClosed())   periode_u= u_max-u_min;
    else periode_u=0;
    if (surface->IsVPeriodic()) periode_v=surface->VPeriod();
    else if (surface->IsVClosed())   periode_v= v_max-v_min;
    else periode_v=0;
    if (type==STANDARD_TYPE(Geom_BSplineSurface)) 
          analyse_bspline();
    
    //**************************************************************
//     ShapeAnalysis::GetFaceUVBounds(face,u_min,u_max,v_min,v_max);
    //**************************************************************

}

OCC_SURFACE::OCC_SURFACE(TopoDS_Face srf, OCC_FONCTION* fonc):MG_SURFACE(),face(srf), fonction1(fonc)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    surface=BRep_Tool::Surface(face);
    Handle(Standard_Type) type=surface->DynamicType(); 
    if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
        {
        Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
        surface = RTSurface->BasisSurface();
        }
    double u1;
    double u2;
    double v1;
    double v2;
    surface->Bounds(u1,u2,v1, v2);
    u_min=u1;
    u_max=u2;
    v_min=v1;
    v_max=v2;
    estperiodeu=surface->IsUClosed();
    estperiodev=surface->IsVClosed();
    if (surface->IsUPeriodic()) periode_u=surface->UPeriod();
    else if (surface->IsUClosed())   periode_u= u_max-u_min;
    else periode_u=0;
    if (surface->IsVPeriodic()) periode_v=surface->VPeriod();
    else if (surface->IsVClosed())   periode_v= v_max-v_min;
    else periode_v=0;
    if (type==STANDARD_TYPE(Geom_BSplineSurface)) 
          analyse_bspline();
    //**************************************************************
//     ShapeAnalysis::GetFaceUVBounds(face,u_min,u_max,v_min,v_max);
    //**************************************************************
    
    
//     std::cout << " " <<std::endl;
//     Handle(Standard_Type) type=surface->DynamicType();
//     if (type==STANDARD_TYPE(Geom_Plane)) std::cout << "plan" << std::endl;
//     else if (type==STANDARD_TYPE(Geom_CylindricalSurface)) std::cout << "Cylindre" << std::endl; 
//     else if (type==STANDARD_TYPE(Geom_ConicalSurface)) std::cout << "Cone" << std::endl; 
//     else if (type==STANDARD_TYPE(Geom_SphericalSurface))  std::cout << "Sphere" << std::endl;
//     else if (type==STANDARD_TYPE(Geom_ToroidalSurface)) std::cout << "Tore" << std::endl;
//     else if (type==STANDARD_TYPE(Geom_BSplineSurface)) std::cout << "Bspline surf" << std::endl;
//     else std::cout << type->Name() << std::endl;
//        
//     std::cout <<  surface->IsUClosed()<< "   " <<surface->IsVClosed() <<std::endl;;
//     std::cout <<  surface->IsUPeriodic() <<"   " <<surface->IsVPeriodic() << std::endl;;
//     if (surface->IsUPeriodic()) std::cout <<  "Pu=" << surface->UPeriod() << std::endl;;
//     if (surface->IsVPeriodic()) std::cout <<  "Pv=" << surface->VPeriod() << std::endl;;
}

OCC_SURFACE::OCC_SURFACE(OCC_SURFACE& mdd):MG_SURFACE(mdd),face(mdd.face), fonction1(mdd.fonction1)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    surface=BRep_Tool::Surface(face);
    Handle(Standard_Type) type=surface->DynamicType(); 
    if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
        {
        Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
        surface = RTSurface->BasisSurface();
        }
    u_min=mdd.get_umin();
    u_max=mdd.get_umax();
    v_min=mdd.get_vmin();
    v_max=mdd.get_vmax();
    estperiodeu=mdd.estperiodeu;
    estperiodev=mdd.estperiodev;
    periode_u=mdd.periode_u;
    periode_v=mdd.periode_v;
}
OCC_SURFACE::~OCC_SURFACE()
{
}
void  OCC_SURFACE::evaluer(double *uv,double *xyz)
{

    //const Handle(Geom_Surface) &surface=BRep_Tool::Surface(face);
    gp_Pnt P;
    double u=uv[0];
    double v=uv[1];

    surface->D0(u,v,P);

    xyz[0]=P.X();
    xyz[1]=P.Y();
    xyz[2]=P.Z();
}
void  OCC_SURFACE::deriver(double *uv,double *xyzdu, double *xyzdv)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    double u=uv[0];
    double v=uv[1];
    gp_Vec D1U;
    gp_Vec D1V;
    gp_Pnt P;

    surface->D1(u,v,P,D1U,D1V);
    xyzdu[0]=D1U.X();
    xyzdu[1]=D1U.Y();
    xyzdu[2]=D1U.Z();

    xyzdv[0]=D1V.X();
    xyzdv[1]=D1V.Y();
    xyzdv[2]=D1V.Z();

}
void  OCC_SURFACE::deriver_seconde(double *uv,double* xyzduu,double* xyzduv,double* xyzdvv,double *xyz, double *xyzdu, double *xyzdv)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    double u=uv[0];
    double v=uv[1];
    gp_Pnt P;
    gp_Vec D1U;
    gp_Vec D1V;
    gp_Vec D2U;
    gp_Vec D2V;
    gp_Vec D2UV;

    surface->D2(u,v,P,D1U,D1V,D2U,D2V,D2UV);

    xyz[0]=P.X();
    xyz[1]=P.Y();
    xyz[2]=P.Z();

    xyzdu[0]=D1U.X();
    xyzdu[1]=D1U.Y();
    xyzdu[2]=D1U.Z();

    xyzdv[0]=D1V.X();
    xyzdv[1]=D1V.Y();
    xyzdv[2]=D1V.Z();

    xyzduu[0]=D2U.X();
    xyzduu[1]=D2U.Y();
    xyzduu[2]=D2U.Z();

    xyzdvv[0]=D2V.X();
    xyzdvv[1]=D2V.Y();
    xyzdvv[2]=D2V.Z();

    xyzduv[0]=D2UV.X();
    xyzduv[1]=D2UV.Y();
    xyzduv[2]=D2UV.Z();

}

void OCC_SURFACE::inverser(double *uv,double *xyz,double precision)
{
   // Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    double u=xyz[0];
    double v=xyz[1];
    double w=xyz[2];
    gp_Pnt P(u,v,w);

    //ShapeAnalysis_Surface SAS(surface);
    //gp_Pnt2d pnt2d=SAS.ValueOfUV(P, precision);
    GeomAPI_ProjectPointOnSurf PPS(P,surface, precision);
    double dist=PPS.LowerDistance();
    if (dist>precision)
        inverser2(uv,xyz,precision);
    else 
        {
        double UU, VV;
        PPS.LowerDistanceParameters(UU,VV);
        uv[0]=UU;
        uv[1]=VV;
        }
}

bool OCC_SURFACE::est_sur_surface(double* xyz, double precision)
{
#ifdef ALL_OCC
  GeomAPI_ProjectPointOnSurf projecteur;
  double u0, u1;
  double v0, v1;
  ShapeAnalysis::GetFaceUVBounds(face,u0,u1,v0,v1);
  gp_Pnt point(xyz[0],xyz[1],xyz[2]);
  projecteur.Init(point,surface,u0,u1,v0,v1);
  if(projecteur.NbPoints()==0) return false;
  if(projecteur.LowerDistance()<precision) return true;
  return false;
#else
  std::cerr << "OCC_SURFACE::est_sur_surface(double* xyz, double precision) --> Non disponible" << std::endl;
#endif
  
}

int OCC_SURFACE::est_periodique_u(void)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    //return surface->IsUClosed();
  return estperiodeu;
}

int OCC_SURFACE::est_periodique_v(void)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
//    return surface->IsVPeriodic();
//    return surface->IsVClosed();
  return estperiodev;
}
double OCC_SURFACE::get_periode_u(void)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    //if (surface->IsUPeriodic()) return surface->UPeriod();
    //if (surface->IsUClosed())   return u_max-u_min;
    //return 1;
  return periode_u;
}

double OCC_SURFACE::get_periode_v(void)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    //if (surface->IsVPeriodic()) return surface->VPeriod();
    //if (surface->IsVClosed()) return v_max-v_min;
    //return 0;
    return periode_v;
    
}

void  OCC_SURFACE::enregistrer(std::ostream& o,double version)
{
  if(fonction1->get_version()=="OCCV2017")
  {
    o <<"%"<<get_id()<< "=SURFACE_OCC("<< get_idoriginal() << ");" << std::endl;
  }
  else
  {
    o <<"%"<<get_id()<< "=SURFACE_OCC("<< fonction1->GetID(face)<< ");" << std::endl;
  }
}
int OCC_SURFACE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
{
    //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
    Handle(Standard_Type) type=surface->DynamicType();
//******plan
    if (type==STANDARD_TYPE(Geom_Plane))
    {
        Handle(Geom_Plane) Pln=Handle(Geom_Plane)::DownCast(surface);
        gp_Pln plan=Pln->Pln();

        double origine[3];
        gp_Pnt centre=plan.Location();

        origine[0]=centre.X();
        origine[1]=centre.Y();
        origine[2]=centre.Z();

        double normal[3];
        gp_Ax1 axe=plan.Axis();
        gp_Dir direction=axe.Direction();

        normal[0]=direction.X();
        normal[1]=direction.Y();
        normal[2]=direction.Z();

        param.ajouter(origine[0]);
        param.ajouter(origine[1]);
        param.ajouter(origine[2]);
        param.ajouter(normal[0]);
        param.ajouter(normal[1]);
        param.ajouter(normal[2]);

        return GEOMETRIE::CONST::Co_PLAN;
    }
//******cylindre
    if (type==STANDARD_TYPE(Geom_CylindricalSurface))
    {
        Handle(Geom_CylindricalSurface) cylinder=Handle(Geom_CylindricalSurface)::DownCast(surface);
        gp_Cylinder cylin=cylinder->Cylinder();

        double origine[3];
        gp_Pnt centre=cylin.Location();
        origine[0]=centre.X();
        origine[1]=centre.Y();
        origine[2]=centre.Z();
        double directionX[3];
        double direction[3];
        double rayon;
        rayon=cylin.Radius();

        gp_Ax3 repere=cylin.Position();
        gp_Dir dirX=repere.XDirection();
        gp_Dir dirZ=repere.Direction();
        directionX[0]=dirX.X();
        directionX[1]=dirX.Y();
        directionX[2]=dirX.Z();
        direction[0]=dirZ.X();
        direction[1]=dirZ.Y();
        direction[2]=dirZ.Z();
        
        
        param.ajouter(origine[0]);
        param.ajouter(origine[1]);
        param.ajouter(origine[2]);
        param.ajouter(directionX[0]);
        param.ajouter(directionX[1]);
        param.ajouter(directionX[2]);
        param.ajouter(direction[0]);
        param.ajouter(direction[1]);
        param.ajouter(direction[2]);
        param.ajouter(rayon);
        param.ajouter(rayon);

        return GEOMETRIE::CONST::Co_CYLINDRE;
    }
//******Cone
    if (type==STANDARD_TYPE(Geom_ConicalSurface))
    {
        Handle(Geom_ConicalSurface) cone=Handle(Geom_ConicalSurface)::DownCast(surface);
        gp_Cone con=cone->Cone();

        double origine[3];
        gp_Pnt centre=con.Location();
        origine[0]=centre.X();
        origine[1]=centre.Y();
        origine[2]=centre.Z();
        double direction[3];
        double directionX[3];
        double rayon;
        rayon=con.RefRadius();
        double angle;
        angle=con.SemiAngle();
        
        gp_Ax3 repere=con.Position();
        gp_Dir dirX=repere.XDirection();
        gp_Dir dirZ=repere.Direction();
        directionX[0]=dirX.X();
        directionX[1]=dirX.Y();
        directionX[2]=dirX.Z();
        direction[0]=dirZ.X();
        direction[1]=dirZ.Y();
        direction[2]=dirZ.Z();
        param.ajouter(origine[0]);
        param.ajouter(origine[1]);
        param.ajouter(origine[2]);
        param.ajouter(directionX[0]);
        param.ajouter(directionX[1]);
        param.ajouter(directionX[2]);
        param.ajouter(direction[0]);
        param.ajouter(direction[1]);
        param.ajouter(direction[2]);
        param.ajouter(rayon);
        param.ajouter(rayon);
        param.ajouter(cos(angle));
        param.ajouter(sin(angle));

        return GEOMETRIE::CONST::Co_CONE;

    }
//*****Sphere
    if (type==STANDARD_TYPE(Geom_SphericalSurface))
    {
        Handle(Geom_SphericalSurface) sphere=Handle(Geom_SphericalSurface)::DownCast(surface);
        gp_Sphere sph=sphere->Sphere();

        double origine[3];
        gp_Pnt centre=sph.Location();
        origine[0]=centre.X();
        origine[1]=centre.Y();
        origine[2]=centre.Z();
        double rayon;
        rayon=sph.Radius();
        param.ajouter(origine[0]);
        param.ajouter(origine[1]);
        param.ajouter(origine[2]);
        param.ajouter(rayon);

        return GEOMETRIE::CONST::Co_SPHERE;

    }
//****** Tore
    if (type==STANDARD_TYPE(Geom_ToroidalSurface))
    {
        Handle(Geom_ToroidalSurface) tore=Handle(Geom_ToroidalSurface)::DownCast(surface);
        gp_Torus tro=tore->Torus();
        double origine[3];
        gp_Pnt centre=tro.Location();
        origine[0]=centre.X();
        origine[1]=centre.Y();
        origine[2]=centre.Z();
        double direction[3];
        double directionX[3];
        double Grayon;
        Grayon=tro.MajorRadius();
        double Prayon;
        Prayon=tro.MinorRadius();
        gp_Ax3 repere=tro.Position();
        gp_Dir dirX=repere.XDirection();
        gp_Dir dirZ=repere.Direction();
        directionX[0]=dirX.X();
        directionX[1]=dirX.Y();
        directionX[2]=dirX.Z();
        direction[0]=dirZ.X();
        direction[1]=dirZ.Y();
        direction[2]=dirZ.Z();
        param.ajouter(origine[0]);
        param.ajouter(origine[1]);
        param.ajouter(origine[2]);
        param.ajouter(directionX[0]);
        param.ajouter(directionX[1]);
        param.ajouter(directionX[2]);
        param.ajouter(direction[0]);
        param.ajouter(direction[1]);
        param.ajouter(direction[2]);
        param.ajouter(Grayon);
        param.ajouter(Prayon);

        return GEOMETRIE::CONST::Co_TORE;

    }
//*******BSpline
    if (type==STANDARD_TYPE(Geom_BSplineSurface))
    {
        Handle(Geom_BSplineSurface) bspline=Handle(Geom_BSplineSurface)::DownCast(surface);

//nombre des noeuds suivant Udirection
        int nb_Uknot=bspline->NbUKnots();
//valeur de Unoeud
        for (int i=1; i<=nb_Uknot; i++)
        {
            double Uvaleur=bspline->UKnot(i);
            param.ajouter(Uvaleur);
        }
//nombre des noeuds suivants Vdirection
        int nb_Vknot=bspline->NbVKnots();
//valeur de Vnoeud
        for (int j=1; j<=nb_Vknot; j++)
        {
            double Vvaleur=bspline->VKnot(j);
            param.ajouter(Vvaleur);
        }
//point de controle et poids
        gp_Pnt pctr;
        double poids;
        for (int u=1; u<=bspline->NbUPoles(); u++)
            for (int v=1; v<=bspline->NbVPoles(); v++)
            {
                pctr=bspline->Pole(u, v);

                param.ajouter(pctr.X());
                param.ajouter(pctr.Y());
                param.ajouter(pctr.Z());

                poids=bspline->Weight(u, v);
                param.ajouter(poids);
            }

        double uDegree=bspline->UDegree();
        param.ajouter(uDegree);
        double vDegree=bspline->VDegree();
        param.ajouter(vDegree);
        return GEOMETRIE::CONST::Co_BSPLINES;
    }
return 0;
}
void OCC_SURFACE::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
{
    //Conversion of the complete geometry of a shape into
    //NURBS geometry
    BRepBuilderAPI_NurbsConvert NURBS(face);
    Handle(Geom_Surface) surface=BRepLib_FindSurface(NURBS).Surface();
    Handle(Geom_BSplineSurface) bspline=GeomConvert::SurfaceToBSplineSurface(surface) ;

// The first parameter indicate the code access
    param.ajouter(2);
//The  follewing two parameters of the list indicate the orders of the net points
    param.ajouter( bspline->UDegree()+1);
    param.ajouter(bspline->VDegree()+1);

//The follewing two parameters indicate the number of rows and colons of the control points
//respectively to the two parameters directions
    param.ajouter(bspline->NbUPoles());
    param.ajouter(bspline->NbVPoles());

// this present the knot vector in the u-direction
    for (unsigned int i=1;i<=bspline->NbUKnots();i++)
    {
        param.ajouter(bspline->UKnot(i));
    }
    //This present the knot vector in the v-direction
    for (unsigned int j=1;j<=bspline->NbVKnots();j++)
    {
        param.ajouter(bspline->VKnot(j));
    }
    for (int v=1;v<=bspline->NbVPoles();v++)
    {
        for (int u=1;u<=bspline->NbUPoles();u++)
        {
            double w=bspline->Weight(u,v);
            gp_Pnt point=bspline->Pole(u, v);
            double x=point.X();
            double y=point.Y();
            double z=point.Z();
            param.ajouter(x);
            param.ajouter(y);
            param.ajouter(z);
            param.ajouter(w);
        }

    }
    indx_premier_ptctr=5+bspline->NbUKnots()+bspline->NbVKnots();


}

void OCC_SURFACE::get_echantillonnage(int numechantillon,std::vector<double> &tab,double eps,double angle_dev)
{
    BRepMesh_IncrementalMesh(face,eps,angle_dev); 
    TopLoc_Location L;
    Handle (Poly_Triangulation) occtri=BRep_Tool::Triangulation(face,L);
    if (occtri.IsNull())
        {
        BRepMesh_IncrementalMesh(face,eps); 
        occtri=BRep_Tool::Triangulation(face,L);
        }
    int nbmaille=occtri->NbTriangles();
    const Poly_Array1OfTriangle& triangles = occtri->Triangles();
    const TColgp_Array1OfPnt2d& uvNodes = occtri->UVNodes();
    int numechantillonpartri=sqrt(numechantillon*1.0/nbmaille)+1;
    for ( Standard_Integer i = 0; i < nbmaille; i++ )
    {
        int n1,n2,n3;
        Poly_Triangle triangle = triangles( i + 1 );
        bool face_reversed = (face.Orientation() == TopAbs_REVERSED);
        if ( face_reversed )
            triangle.Get( n1, n3, n2 );
        else
            triangle.Get( n1, n2, n3 );
        gp_Pnt2d p1=uvNodes(n1);
        double u1=p1.X();
        double v1=p1.Y();
        gp_Pnt2d p2=uvNodes(n2);
        double u2=p2.X();
        double v2=p2.Y();
        gp_Pnt2d p3=uvNodes(n3);
        double u3=p3.X();
        double v3=p3.Y();
        for (int ii=0;ii<numechantillonpartri;ii++)
            for (int jj=0;jj<numechantillonpartri;jj++)
                {
                double ksi=(ii*1.0+0.5)/(double)numechantillonpartri;    
                double eta=(jj*1.0+0.5)/(double)numechantillonpartri; 
                if (ksi+eta>0.999999999999) jj=numechantillonpartri;
                else 
                    {
                    double u=ksi*u1+eta*u2+(1.-ksi-eta)*u3;
                    double v=ksi*v1+eta*v2+(1.-ksi-eta)*v3;
                    gp_Pnt P;
                    surface->D0(u,v,P);
                    tab.push_back(P.X());
                    tab.push_back(P.Y());
                    tab.push_back(P.Z());
                    tab.push_back(u);
                    tab.push_back(v);
                    }
                }
            
        
    }
    
}


void OCC_SURFACE::get_triangulation(MG_MAILLAGE* mai,MG_FACE* mgface,std::multimap<double,MG_NOEUD*,std::less<double> >& tabnoeudfus,double eps,int mode)
{
    TPL_MAP_ENTITE<MG_SOMMET*> listsom;
    TPL_MAP_ENTITE<MG_ARETE*> listare;
    //std::map<unsigned long,bool> aretemaille;
    std::map<unsigned long,std::map<double,NOEUDARETE,std::less<double> >,std::less<unsigned long> >      areteamaille;
    int nbboucle=mgface->get_nb_mg_boucle();
    for (int i=0;i<nbboucle;i++)
      {
      MG_BOUCLE* bou=mgface->get_mg_boucle(i);
      int nbarete=bou->get_nb_mg_coarete();
      for (int j=0;j<nbarete;j++)
        {
        MG_ARETE* are=bou->get_mg_coarete(j)->get_arete();
        listare.ajouter(are);
        //bool amailler=true;
        //if (are->get_lien_maillage()->get_nb()>0) amailler=false;
        //aretemaille[are->get_id()]=amailler;
        listsom.ajouter(bou->get_mg_coarete(j)->get_arete()->get_cosommet1()->get_sommet());  
        listsom.ajouter(bou->get_mg_coarete(j)->get_arete()->get_cosommet2()->get_sommet());
        std::map<double,NOEUDARETE,std::less<double> > tmp;
        std::pair<unsigned long,std::map<double,NOEUDARETE,std::less<double> > > maptmp(are->get_id(),tmp);
        areteamaille.insert(maptmp);
        }
      }


    TopLoc_Location L;
    Handle (Poly_Triangulation) pt=BRep_Tool::Triangulation(face,L);
    int nbnoeud=pt->NbNodes();
    int nbmaille=pt->NbTriangles();
    const TColgp_Array1OfPnt& nodes = pt->Nodes();
    const Poly_Array1OfTriangle& triangles = pt->Triangles();
    const TColgp_Array1OfPnt2d& uvNodes = pt->UVNodes();
    std::vector<MG_NOEUD*> tabnoeud;
    for ( Standard_Integer i = 0; i < nbnoeud; i++ )
    {
        gp_Pnt p1=nodes(i+1);
        double xx=p1.X();
        double yy=p1.Y();
        double zz=p1.Z();
        double key=fabs(xx)+fabs(yy)+fabs(zz);
        MG_NOEUD* nvnoeud=NULL;
        if (mode>1)
        {
                std::multimap<double,MG_NOEUD*,std::less<double> >::iterator it,itbas,ithaut;
                itbas=tabnoeudfus.lower_bound(key*0.99);
                ithaut=tabnoeudfus.upper_bound(key*1.1010101);
                for ( it=itbas ; it != ithaut; it++ )
                {
                        MG_NOEUD* ntmp=(*it).second;
                        double xtmp=ntmp->get_x();
                        double ytmp=ntmp->get_y();
                        double ztmp=ntmp->get_z();
                        OT_VECTEUR_3D vec(xtmp-xx,ytmp-yy,ztmp-zz);
                        if (vec.get_longueur()<1e-6*eps) {
                                nvnoeud=ntmp;
                                break;
                        }
                }
        }
        if (nvnoeud==NULL)
        {
                MG_ELEMENT_TOPOLOGIQUE *topo=mgface;
                TPL_MAP_ENTITE<MG_SOMMET*>::ITERATEUR it1;
                for (MG_SOMMET* som=listsom.get_premier(it1);som!=NULL;som=listsom.get_suivant(it1))
                {
                        double xyz[3];
                        som->get_point()->evaluer(xyz);
                        OT_VECTEUR_3D vec(xyz[0]-xx,xyz[1]-yy,xyz[2]-zz);
                        if (vec.get_longueur()<1e-6*eps)
                        {
                                topo=som;
                                break;

                        }
                }
                double param_t;
                if (topo==mgface)
                {
                        TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR it2;
                        for (MG_ARETE* are=listare.get_premier(it2);are!=NULL;are=listare.get_suivant(it2))
                        {
                                double t;
                                double xyz[3]={xx,yy,zz};
                                are->inverser(t,xyz);
                                if (are->get_courbe()->est_periodique())
                                        if (t< are->get_tmin()) t=t+are->get_courbe()->get_periode();
                                double xyztmp[3];
                                are->evaluer(t,xyztmp );
                                OT_VECTEUR_3D vec(xyz,xyztmp);
                                if (vec.get_longueur()<1e-6*eps)
                                        if ((t>are->get_tmin()) && (t<are->get_tmax()))
                                        {
                                                topo=are;
                                                param_t=t;
                                                break;
                                        }
                        }
                }
                nvnoeud=new MG_NOEUD(topo,xx,yy,zz,MAGIC::ORIGINE::TRIANGULATION);
                mai->ajouter_mg_noeud(nvnoeud);
                std::pair<double,MG_NOEUD*> tmp(key,nvnoeud);
                tabnoeudfus.insert(tmp);
                if (topo->get_dimension()==1)
                {
                        NOEUDARETE na;
                        na.no=nvnoeud;
                        na.are=(MG_ARETE*)topo;
                        na.t=param_t;
                        std::pair<double,NOEUDARETE> tmp(na.t,na);
                        areteamaille[topo->get_id()].insert(tmp);
                }
        }


        tabnoeud.insert(tabnoeud.end(),nvnoeud);
    }
    for ( Standard_Integer i = 0; i < nbmaille; i++ )
    {
        int n1,n2,n3;
        Poly_Triangle triangle = triangles( i + 1 );
        bool face_reversed = (face.Orientation() == TopAbs_REVERSED);
        if ( face_reversed )
            triangle.Get( n1, n3, n2 );
        else
            triangle.Get( n1, n2, n3 );
        MG_NOEUD* noeud1=tabnoeud[n1-1];
        MG_NOEUD* noeud2=tabnoeud[n2-1];
        MG_NOEUD* noeud3=tabnoeud[n3-1];
        if (noeud1==noeud2) continue;
        if (noeud1==noeud3) continue;
        if (noeud2==noeud3) continue;
        mai->ajouter_mg_triangle(mgface,noeud1,noeud2,noeud3,MAGIC::ORIGINE::TRIANGULATION);
        /*if (noeud1->get_lien_topologie()->get_dimension()==0) 
            {
            MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();  
            som->get_lien_maillage()->ajouter(noeud1);
            }
        if (noeud2->get_lien_topologie()->get_dimension()==0) 
            {
            MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();  
            som->get_lien_maillage()->ajouter(noeud2);
            }
        if (noeud3->get_lien_topologie()->get_dimension()==0) 
            {
            MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();  
            som->get_lien_maillage()->ajouter(noeud3);
            }*/
        /*if (noeud1->get_lien_topologie()==noeud2->get_lien_topologie())
            if (noeud1->get_lien_topologie()->get_dimension()==1)
              {
              MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
              if (aretemaille[are->get_id()]==true) 
                {
                MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud2,TRIANGULATION);
                are->get_lien_maillage()->ajouter(seg);
                }
              }
        if (noeud1->get_lien_topologie()==noeud3->get_lien_topologie())
            if (noeud1->get_lien_topologie()->get_dimension()==1)
              {
              MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
              if (aretemaille[are->get_id()]==true) 
                {
                MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud3,TRIANGULATION);
                are->get_lien_maillage()->ajouter(seg);
                }
              }#include <../gmsh/tutorial/t8.geo>
          if (noeud3->get_lien_topologie()==noeud2->get_lien_topologie())
            if (noeud3->get_lien_topologie()->get_dimension()==1)
              {
              MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
              if (aretemaille[are->get_id()]==true) 
                {
                MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud2,TRIANGULATION);
                are->get_lien_maillage()->ajouter(seg);
                }
              }
          if (noeud1->get_lien_topologie()->get_dimension()==1)
            if (noeud2->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud2,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              }
           if (noeud1->get_lien_topologie()->get_dimension()==1)
            if (noeud3->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud1->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud1,noeud3,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              }  
           if (noeud2->get_lien_topologie()->get_dimension()==1)
            if (noeud1->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud2->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud2,noeud1,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              }
           if (noeud2->get_lien_topologie()->get_dimension()==1)
            if (noeud3->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud2->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud3->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud2,noeud3,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              } 
           if (noeud3->get_lien_topologie()->get_dimension()==1)
            if (noeud1->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud1->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud1,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              }
           if (noeud3->get_lien_topologie()->get_dimension()==1)
            if (noeud2->get_lien_topologie()->get_dimension()==0)
              {
              MG_ARETE* are=(MG_ARETE*)noeud3->get_lien_topologie();
              MG_SOMMET* som=(MG_SOMMET*)noeud2->get_lien_topologie();
              if ((are->get_cosommet1()->get_sommet()==som)||(are->get_cosommet2()->get_sommet()==som))
                  {
                  if (aretemaille[are->get_id()]==true) 
                    {
                    MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeud3,noeud2,TRIANGULATION);
                    are->get_lien_maillage()->ajouter(seg);
                    }  
                  }
              }*/   
    }
  if (mode>1)
  {
  TPL_MAP_ENTITE<MG_ARETE*>::ITERATEUR it2;
  for (MG_ARETE* are=listare.get_premier(it2);are!=NULL;are=listare.get_suivant(it2))  
    {
    if (are->get_lien_maillage()->get_nb()==0)
      {
      unsigned long id=are->get_id();
      MG_NOEUD* nodep=(MG_NOEUD*)are->get_cosommet1()->get_sommet()->get_lien_maillage()->get(0);
      MG_NOEUD* noarr=(MG_NOEUD*)are->get_cosommet2()->get_sommet()->get_lien_maillage()->get(0);
      std::map<double,NOEUDARETE,std::less<double> >::iterator it=areteamaille[id].begin();
      MG_NOEUD* noeudcourant=nodep;
      while (it!=areteamaille[id].end())
          {
          MG_NOEUD* noeud=(*it).second.no;
          //MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeudcourant,noeud,TRIANGULATION);
          MG_SEGMENT* seg=mai->get_mg_segment(noeudcourant->get_id(),noeud->get_id());
          seg->change_lien_topologie(are);
          noeudcourant=noeud;
          it++;
          }
      //MG_SEGMENT* seg=mai->ajouter_mg_segment(are,noeudcourant,noarr,TRIANGULATION);
      MG_SEGMENT* seg=mai->get_mg_segment(noeudcourant->get_id(),noarr->get_id());
      seg->change_lien_topologie(are);
      }
    }
    
  }
}

void OCC_SURFACE::get_liste_pole(std::vector< double> *liste_pole,double eps)
{  
  //Handle(Geom_Surface) surface=BRep_Tool::Surface(face);
  Handle(Standard_Type) type=surface->DynamicType(); 
  //std::cout << surface->DynamicType()->Name()    <<std::endl;
  if(type==STANDARD_TYPE(Geom_RectangularTrimmedSurface))
  {
     Handle(Geom_RectangularTrimmedSurface) RTSurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
     Handle(Geom_Surface) basissurface = RTSurface->BasisSurface();
     Handle(Standard_Type) type_basissurface=basissurface->DynamicType();
     if(type_basissurface==STANDARD_TYPE(Geom_SphericalSurface))
     {
       type=STANDARD_TYPE(Geom_SphericalSurface);
    }
  }
  if(type==STANDARD_TYPE(Geom_SphericalSurface))
  {    
    BRepClass_FaceClassifier faceclassifier;
    gp_Pnt2d pnt2d_pole_sud(0.0,-M_PI/2.);
    faceclassifier.Perform(face,pnt2d_pole_sud,eps);
    if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
    {
      liste_pole->push_back(0.0);
      liste_pole->push_back(-M_PI/2.);
    }    
    gp_Pnt2d pnt2d_pole_nord(0.0,M_PI/2.);
    faceclassifier.Perform(face,pnt2d_pole_nord,eps);
    if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
    {
      liste_pole->push_back(0.0);
      liste_pole->push_back(M_PI/2.);
    }
  }
  if(type==STANDARD_TYPE(Geom_SurfaceOfRevolution))
  {
    BRepClass_FaceClassifier faceclassifier;
    gp_Pnt2d pnt2d_pole_sud(0.0,v_min);
    faceclassifier.Perform(face,pnt2d_pole_sud,eps);
    if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
    {
      liste_pole->push_back(0.0);
      liste_pole->push_back(v_min);
    }    
    gp_Pnt2d pnt2d_pole_nord(0.0,v_max);
    faceclassifier.Perform(face,pnt2d_pole_nord,eps);
    if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
    {
      liste_pole->push_back(0.0);
      liste_pole->push_back(v_max);
    }
  }
//   if(type==STANDARD_TYPE(Geom_ConicalSurface))
//   {    
//     BRepClass_FaceClassifier faceclassifier;
//     gp_Pnt2d pnt2d_pole(0.0,v_max);
//     faceclassifier.Perform(face,pnt2d_pole,eps);
//     if(faceclassifier.State()==TopAbs_IN || faceclassifier.State()==TopAbs_ON)
//     {
//       liste_pole->push_back(0.0);
//       liste_pole->push_back(v_max);
//     }    
//   }
}


void OCC_SURFACE::analyse_bspline(void)
{
  Handle(Geom_BSplineSurface) bspline=Handle(Geom_BSplineSurface)::DownCast(surface);
//point de controle et poids
        bool ferme=true;
            for (int v=1; v<=bspline->NbVPoles(); v++)
            {
                gp_Pnt pctr1=bspline->Pole(1, v);
                gp_Pnt pctr2=bspline->Pole(bspline->NbUPoles(), v);
                double poids1=bspline->Weight(1, v);
                double poids2=bspline->Weight(bspline->NbUPoles(), v);
                double xyz1[4];
                double xyz2[4];
                xyz1[0]=pctr1.X();
                xyz1[1]=pctr1.Y();
                xyz1[2]=pctr1.Z();
                xyz1[3]=poids1;
                xyz2[0]=pctr2.X();
                xyz2[1]=pctr2.Y();
                xyz2[2]=pctr2.Z();
                xyz2[3]=poids2;
                OT_VECTEUR_4D vec(xyz1,xyz2);
                 if (vec.get_longueur()>fonction1->get_precision()) 
                   ferme=false;
              
            }
        if (ferme)
            {
            estperiodeu=1;
            periode_u=u_max-u_min;
            }
            ferme=true;
            for (int u=1; u<=bspline->NbUPoles(); u++)
            {
                gp_Pnt pctr1=bspline->Pole(u,1);
                gp_Pnt pctr2=bspline->Pole(u,bspline->NbVPoles());
                double poids1=bspline->Weight(u,1);
                double poids2=bspline->Weight(u,bspline->NbVPoles());
                double xyz1[4];
                double xyz2[4];
                xyz1[0]=pctr1.X();
                xyz1[1]=pctr1.Y();
                xyz1[2]=pctr1.Z();
                xyz1[3]=poids1;
                xyz2[0]=pctr2.X();
                xyz2[1]=pctr2.Y();
                xyz2[2]=pctr2.Z();
                xyz2[3]=poids2;
                OT_VECTEUR_4D vec(xyz1,xyz2);
                 if (vec.get_longueur()>fonction1->get_precision())
                   ferme=false;
              
            }
        if (ferme)
            {
            estperiodev=1;
            periode_v=v_max-v_min;
            }
}

void OCC_SURFACE::inverser2(double *uv,double *xyz,double precision)
{
int NUMDECOUP=20;
double u1;
double u2;
double v1;
double v2;
double u,v,distmin=1e300;
surface->Bounds(u1,u2,v1, v2);    
for (int i=0;i<NUMDECOUP+1;i++)
    for (int j=0;j<NUMDECOUP+1;j++)
        {
        double ut=u1+i*1.0/NUMDECOUP*(u2-u1);    
        double vt=v1+j*1.0/NUMDECOUP*(v2-v1);
        gp_Pnt P;
        surface->D0(ut,vt,P);
        double xyz2[3]={P.X(),P.Y(),P.Z()};
        OT_VECTEUR_3D vec(xyz,xyz2);
        double dist=vec.get_longueur();
        if (dist<distmin)
                {
                distmin=dist;
                u=ut;
                v=vt;
                }
        }
bool sortie=false;
int nbiteration=0;
while (sortie==false)
{
double uv2[3]={u,v,0.};
double xyzuv[3],xyzuvdu[3],xyzuvdv[3],xyzuvduu[3],xyzuvduv[3],xyzuvdvv[3];
deriver_seconde(uv2,xyzuvduu,xyzuvduv,xyzuvdvv,xyzuv,xyzuvdu,xyzuvdv);
OT_VECTEUR_3D vxyz(xyz);
OT_VECTEUR_3D vxyzuv(xyzuv);
OT_VECTEUR_3D P=vxyzuv-vxyz;
OT_VECTEUR_3D vxyzuvdu(xyzuvdu);
OT_VECTEUR_3D vxyzuvdv(xyzuvdv);
OT_VECTEUR_3D vxyzuvduu(xyzuvduu);
OT_VECTEUR_3D vxyzuvdvv(xyzuvdvv);
OT_VECTEUR_3D vxyzuvduv(xyzuvduv);
double Fdu=P*vxyzuvduu+vxyzuvdu*vxyzuvdu;
double Fdv=P*vxyzuvduv+vxyzuvdv*vxyzuvdu;
double Gdu=P*vxyzuvduv+vxyzuvdu*vxyzuvdv;
double Gdv=P*vxyzuvdvv+vxyzuvdv*vxyzuvdv;
double F=P*xyzuvdu;
double G=P*xyzuvdv;
double det=Fdu*Gdv-Gdu*Fdv;
double deltau=(-F*Gdv+Gdu*G)/det;
double deltav=(-Fdu*G+F*Fdv)/det;
double new_u=u+deltau;
double new_v=v+deltav;
if (new_u<u1) new_u=u1;
if (new_u>u2) new_u=u2;
if (new_v<v1) new_v=v1;
if (new_v>v2) new_v=v2;
if (fabs(new_u-u)<precision)
    if (fabs(new_v-v)<precision)
        sortie=true;
u=new_u;
v=new_v;
nbiteration++;
if (nbiteration>1000) sortie=true;
}
uv[0]=u;
uv[1]=v;
}

#endif
Revision:	1149
Committed:	Fri May 24 20:23:11 2024 UTC (11 months, 4 weeks ago) by francois
File size:	37405 byte(s)
Log Message:	mise de constante en namespace