mailleur/src/mailleur_stl.cpp

#include "mailleur_stl.h"
#include "mg_maillage.h"
#include "tpl_fonctions_generiques.h"
#include "fct_taille.h"
#include <tpl_map_entite.h>
#include <mg_maillage_outils.h>
#include <fct_taille_fem_solution_generateur_gradient.h>


MAILLEUR_STL::MAILLEUR_STL(MG_MAILLAGE* maiori, MG_GESTIONNAIRE* gt, FCT_TAILLE* carte,double limite,double angle1,double angle2):MAILLEUR(false),maiin(maiori),gest(gt),metrique(carte),limite_discretisation(limite),angle_arete(angle1),angle_limite(angle2),pas(4),nx(25),ny(25),nz(25),qualmaxaopt(0.5),nbcoucheopt(1)
{
}

MAILLEUR_STL::MAILLEUR_STL(MAILLEUR_STL& mdd):MAILLEUR(mdd),maiin(mdd.maiin),mai(mdd.mai),metrique(mdd.metrique),pas(mdd.pas),angle_arete(mdd.angle_arete),angle_limite(mdd.angle_limite),limite_discretisation(mdd.limite_discretisation),nx(mdd.nx),ny(mdd.ny),nz(mdd.nz),qualmaxaopt(mdd.qualmaxaopt),nbcoucheopt(mdd.nbcoucheopt)
{
}


MAILLEUR_STL::~MAILLEUR_STL()
{
delete gesttmp;
}


int MAILLEUR_STL::maille(class MG_GROUPE_TOPOLOGIQUE* mggt)
{
    
affiche((char*)"  Parametre du mailleur");
char mess[255];
sprintf(mess,(char*)"       Ratio limite de discrétisation = %f",limite_discretisation);
affiche(mess);
sprintf(mess,(char*)"       Ratio limite angle arete = %f",angle_arete);
affiche(mess);
sprintf(mess,(char*)"       Ratio limite angle de la cavité de Delaunay = %f",angle_limite);
affiche(mess);
affiche((char*)"  Duplication du maillage original");
//mai=maiin->dupliquer(gest);
mai=maiin->detacher(gest);
affiche((char*)"  Adaptation de la carte de taille");
adapte_carte();
affiche((char*)"  Prepartion de la projection");
prepare_projection();
affiche((char*)"  Determination des arêtes");
determine_arete();
affiche((char*)"  Maillage");
cree_maillage();
affiche((char*)"  Optimisation de maillage");
int oldnb=mai->get_nb_mg_triangle();
for (int i=0;i<nbcoucheopt;i++)
    {
    sprintf(mess,"    Passe numero %d",i+1);
    affiche(mess);
    int nb=optimise_maillage();
    if ((nb==1)||(nb==oldnb)) i=nbcoucheopt;
    oldnb=nb;  
    }
return OK;
}

int MAILLEUR_STL::optimise_maillage(void)
{
std::multimap<double,MG_TRIANGLE*> lsttriaopt;
LISTE_MG_TRIANGLE::iterator it;    
for (MG_TRIANGLE* tri=mai->get_premier_triangle(it);tri!=NULL;tri=mai->get_suivant_triangle(it))
    {
    double val=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
    lsttriaopt.insert(std::pair<double,MG_TRIANGLE*>(val,tri));
    tri->change_solution(val);
    }
int ok=0;
int compteur_iter=0;
do
{
compteur_iter++;
if (lsttriaopt.begin()->first<qualmaxaopt)
    {
    MG_TRIANGLE *tri=lsttriaopt.begin()->second;
    double crit1,crit2,crit3;
    double x1,y1,z1;
    double x2,y2,z2;
    double x3,y3,z3;
    int res1=bouge_point(tri->get_noeud1(),crit1,x1,y1,z1);
    int res2=bouge_point(tri->get_noeud2(),crit2,x2,y2,z2);
    int res3=bouge_point(tri->get_noeud3(),crit3,x3,y3,z3);
    double critmax=0;
    double nvx,nvy,nvz;
    MG_NOEUD* nno;
    if (res1==1)
        if (crit1>critmax) 
        {critmax=crit1;nvx=x1;nvy=y1;nvz=z1;nno=tri->get_noeud1();}
    if (res2==1)
        if (crit2>critmax) 
        {critmax=crit2;nvx=x2;nvy=y2;nvz=z2;nno=tri->get_noeud2();}
    if (res3==1)
        if (crit3>critmax) 
        {critmax=crit3;nvx=x3;nvy=y3;nvz=z3;nno=tri->get_noeud3();}
    if (critmax>lsttriaopt.begin()->first)
        {
        projete_point_sur_triangulation(nvx,nvy,nvz,nvx,nvy,nvz);
        nno->change_x(nvx);
        nno->change_y(nvy);
        nno->change_z(nvz);
        for (int i=0;i<nno->get_lien_triangle()->get_nb();i++)
            {
            MG_TRIANGLE* tri=nno->get_lien_triangle()->get(i);
            double *xyz1=tri->get_noeud1()->get_coord();
            double *xyz2=tri->get_noeud2()->get_coord();
            double *xyz3=tri->get_noeud3()->get_coord();
            double val=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3); 
            std::multimap<double,MG_TRIANGLE*>:: iterator it=lsttriaopt.lower_bound(tri->get_solution());
            int trouve=0;
            while (trouve==0) 
                if (it->first>tri->get_solution()) trouve=-1;
                else if (it->second==tri) trouve=1;
                    else it++;
            if (trouve==1) lsttriaopt.erase(it);
            //lsttriaopt.insert(std::pair<double,MG_TRIANGLE*>(val,tri));
            //tri->change_solution(val);
            }
            
        }
    else lsttriaopt.erase(lsttriaopt.begin());
             
    
    }
else 
    ok=1;
}
while (ok==0);
 char chaine[500];
sprintf(chaine,"     %d iterations effectuées",compteur_iter);
affiche(chaine); 
return(compteur_iter);
}


void MAILLEUR_STL::adapte_carte(void)
{
LISTE_MG_SEGMENT::iterator it;
std::vector<double> lstraf;
for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
    {
    double val=metrique->calcul_distance_metrique(seg,pas);
    if (val<1.)
        {
        double l=seg->get_longueur()/limite_discretisation/1.01;    
        double *xyz=seg->get_noeud1()->get_coord();
        lstraf.push_back(xyz[0]);
        lstraf.push_back(xyz[1]);
        lstraf.push_back(xyz[2]);
        lstraf.push_back(l);
        xyz=seg->get_noeud2()->get_coord();
        lstraf.push_back(xyz[0]);
        lstraf.push_back(xyz[1]);
        lstraf.push_back(xyz[2]);
        lstraf.push_back(l);
        }
    }
gesttmp=new MG_GESTIONNAIRE();
FCT_TAILLE_FEM_SOLUTION_GENERATEUR_GRADIENT* metriquetmp=new FCT_TAILLE_FEM_SOLUTION_GENERATEUR_GRADIENT(gesttmp,mai,metrique,2,nx,ny,nz,(char*)"carteadapte.sol",lstraf);
metriquetmp->construit();        
metrique2=metrique;
metrique=metriquetmp;
}


void MAILLEUR_STL::prepare_projection()
{
LISTE_MG_NOEUD::iterator it;
MG_NOEUD* no=maiin->get_premier_noeud(it);
BOITE_3D boite(no->get_boite_3D());
TPL_LISTE_ENTITE<MG_NOEUD*> lstini;
for (MG_NOEUD* no=maiin->get_premier_noeud(it);no!=NULL;no=maiin->get_suivant_noeud(it))
    {
    boite=no->get_boite_3D()+boite;
    lstini.ajouter(no);
    }
boite.change_grosseur(1.2);
octree.initialiser(&lstini,1,boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax());
LISTE_MG_TRIANGLE::iterator itt;
for (MG_TRIANGLE* tri=maiin->get_premier_triangle(itt);tri!=NULL;tri=maiin->get_suivant_triangle(itt))
        octree.inserer(tri);
}


void MAILLEUR_STL::determine_arete(void)
{
LISTE_MG_SEGMENT::iterator it;
int nbsegmentarete=0;
for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
    {
    int nbtri=seg->get_lien_triangle()->get_nb();
        if (nbtri!=2) 
            {
            char chaine[500];
            sprintf(chaine,"       Maillage STL invalide : segment à %d voisin : noeud 1 %lu noeud 2 %lu ",nbtri,seg->get_noeud1()->get_id(),seg->get_noeud2()->get_id());   
            affiche(chaine);
            }
    if (nbtri==2)
        {
        MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);    
        MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);    
        double angle=get_angle_par_noeud<MG_NOEUD*>(tri1->get_noeud1(),tri1->get_noeud2(),tri1->get_noeud3(),tri2->get_noeud1(),tri2->get_noeud2(),tri2->get_noeud3());
        if ((angle<M_PI-angle_arete) || (angle > M_PI+angle_arete))
            {
            seg->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
            seg->get_noeud1()->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
            seg->get_noeud2()->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
            nbsegmentarete++;
            }
        ajouter_seg_a_discretiser(seg);
        }
    }
char mess[255];
sprintf(mess,(char*)"       Nombre de segments aretes = %d",nbsegmentarete);
affiche(mess);
}


void MAILLEUR_STL::calcul_nouveau_noeud_et_etoile(MG_SEGMENT* seg,double& x, double& y, double& z,OT_VECTEUR_3D &normale, TPL_LISTE_ENTITE<MG_TRIANGLE*> &listuniquetri,std::map<unsigned long,std::pair<MG_SEGMENT*,int> > &lstsegment)
{
MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);    
MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);   
double xyz1[3],xyz2[3];
tri1->normal(xyz1);
tri2->normal(xyz2);
OT_VECTEUR_3D nor1(xyz1);
OT_VECTEUR_3D nor2(xyz2);
normale=0.5*nor1+0.5*nor2;
MG_NOEUD* no1=seg->get_noeud1();
MG_NOEUD* no2=seg->get_noeud2();
x=0.5*(no1->get_x()+no2->get_x());
y=0.5*(no1->get_y()+no2->get_y());
z=0.5*(no1->get_z()+no2->get_z());
listuniquetri.ajouter(tri1);
listuniquetri.ajouter(tri2);
for (int i=0;i<listuniquetri.get_nb();i++)
     {
     MG_TRIANGLE* tri=listuniquetri.get(i);
     for (int cote=1;cote<4;cote++)
     {
        MG_TRIANGLE* triv=get_voisin(tri,cote);
        if (!(listuniquetri.est_dans_la_liste(triv)))    
            if (triv!=NULL)
            {
            double xyzv[3];
            triv->normal(xyzv);
            OT_VECTEUR_3D norv(xyzv);
            if (norv*normale>cos(angle_limite))
                {
                bool res=respecte_delaunay(x,y,z,triv); 
                if (!res)
                    listuniquetri.ajouter(triv);    
                }
            }
     }
     

    }
for (int i=0;i<listuniquetri.get_nb();i++)
     {
     MG_TRIANGLE* tri=listuniquetri.get(i);
     std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp;
     MG_SEGMENT* seg=tri->get_segment1();
     ptmp=lstsegment.find(seg->get_id());
     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
     else ptmp->second.second=ptmp->second.second+1;
     seg=tri->get_segment2();
     ptmp=lstsegment.find(seg->get_id());
     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
     else ptmp->second.second=ptmp->second.second+1;
     seg=tri->get_segment3();
     ptmp=lstsegment.find(seg->get_id());
     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
     else ptmp->second.second=ptmp->second.second+1;
     }
 
 
projete_point_sur_triangulation(x,y,z,x,y,z); 
 
 
bool convergecoquille; 
 
do
{ 
convergecoquille=true;
TPL_LISTE_ENTITE<MG_SEGMENT*> contour;
for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
    if (ptmp->second.second==1)
        {
        MG_SEGMENT* seg=ptmp->second.first;   
        contour.ajouter(seg);
        }
        
for (int i=0;i<contour.get_nb();i++)
    {
    MG_SEGMENT* seg=contour.get(i);
    MG_TRIANGLE *triin,*triext;
    if (listuniquetri.est_dans_la_liste(seg->get_lien_triangle()->get(0)))
        {
        triin=seg->get_lien_triangle()->get(0);    
        triext=seg->get_lien_triangle()->get(1);    
        }
        else
        {
        triin=seg->get_lien_triangle()->get(1);    
        triext=seg->get_lien_triangle()->get(0);    
        }
    MG_NOEUD* no1=seg->get_noeud1();
    MG_NOEUD* no2=seg->get_noeud2();
    int sens=-1;
    if ((triin->get_noeud1()==no1) && (triin->get_noeud2()==no2)) sens=1;
    if ((triin->get_noeud2()==no1) && (triin->get_noeud3()==no2)) sens=1;
    if ((triin->get_noeud3()==no1) && (triin->get_noeud1()==no2)) sens=1;
    double *xyz1,*xyz2;
    if (sens==-1)
        {
        xyz1=no2->get_coord();
        xyz2=no1->get_coord();
        }
    else
        {
        xyz1=no1->get_coord();
        xyz2=no2->get_coord();
        }
    double xyz3[3]={x,y,z};
    OT_VECTEUR_3D vec1(xyz1,xyz2);
    OT_VECTEUR_3D vec2(xyz1,xyz3);
    OT_VECTEUR_3D nor=vec1&vec2;
    nor.norme();
    double xyzin[3];
    triin->normal(xyzin);
    OT_VECTEUR_3D norin(xyzin);
    if (nor*norin<0.) 
        {
        affiche("etoile non conforme");
        listuniquetri.ajouter(triext);    
        std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp;
        MG_SEGMENT* seg=triext->get_segment1();
        ptmp=lstsegment.find(seg->get_id());
        ptmp=lstsegment.find(seg->get_id());
        if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
        else ptmp->second.second=ptmp->second.second+1;
        seg=triext->get_segment2();
        ptmp=lstsegment.find(seg->get_id());
        if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
        else ptmp->second.second=ptmp->second.second+1;
        seg=triext->get_segment3();
        ptmp=lstsegment.find(seg->get_id());
        if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
        else ptmp->second.second=ptmp->second.second+1;
        convergecoquille=false;
        }
    }
}
while (convergecoquille==false);
        
 /*       
MG_SEGMENT* segcourant=contour.get(0);
MG_NOEUD* n1=segcourant->get_noeud1();
MG_NOEUD* n2=segcourant->get_noeud2();
MG_NOEUD* ns;
MG_NOEUD* premier=n2;

do
{
MG_SEGMENT* segsuivant=NULL;
int i=0;
while (segsuivant==NULL)
    {
    if (contour.get(i)!=segcourant)    
        if ((contour.get(i)->get_noeud1()==n2) || ((contour.get(i)->get_noeud2()==n2)))
            {   
            segsuivant=contour.get(i);
            if (segsuivant->get_noeud1()==n2) ns=segsuivant->get_noeud2(); else ns=segsuivant->get_noeud1();
            }
     i++;
     }
OT_VECTEUR_3D vecbase(n1->get_coord(),n2->get_coord());
OT_VECTEUR_3D vecsuiv(n2->get_coord(),ns->get_coord());
OT_VECTEUR_3D vecori(n1->get_coord(),ns->get_coord());
OT_VECTEUR_3D sens(x-0.5*(n1->get_x()+n2->get_x()),y-0.5*(n1->get_y()+n2->get_y()),z-0.5*(n1->get_z()+n2->get_z()));
vecbase.norme();
vecsuiv.norme();
vecori.norme();
sens.norme();
OT_VECTEUR_3D norm=vecbase&sens;
sens=norm&vecbase;
double angle=acos((-vecbase)*vecsuiv);
double signe=vecori*sens;
if (signe<0.) angle=2*M_PI-angle;
if (angle>M_PI) 
    affiche("etoile concave");
segcourant=segsuivant;
n1=n2;
n2=ns;
}
while (n2!=premier);
*/
}

int MAILLEUR_STL::bouge_point(MG_NOEUD* mg_noeud,double& crit,double& x,double& y, double& z)
{
if (mg_noeud->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE) 
    return 0;
TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
double critdep=1.;
for (int i=0;i<mg_noeud->get_lien_triangle()->get_nb();i++)
    {
    MG_TRIANGLE* tri=mg_noeud->get_lien_triangle()->get(i);
    double val=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
    if (val<critdep) critdep=val;
    if (tri->get_noeud1()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud1());
    if (tri->get_noeud2()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud2());
    if (tri->get_noeud3()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud3());
    }
TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it;
double xi=0.;
double yi=0.;
double zi=0.;
for (MG_NOEUD* no=lstnoeud.get_premier(it);no!=NULL;no=lstnoeud.get_suivant(it))
    {
    xi=xi+no->get_x();    
    yi=yi+no->get_y();    
    zi=zi+no->get_z();    
    }
xi=xi/lstnoeud.get_nb();
yi=yi/lstnoeud.get_nb();
zi=zi/lstnoeud.get_nb();
OT_VECTEUR_3D dir(xi-mg_noeud->get_x(),yi-mg_noeud->get_y(),zi-mg_noeud->get_z());
if (dir.get_longueur()<1e-10) return 0;
double alpha=0.1;
bool converge=false;
double critencours=critdep;
int passe=0;
while (converge==false)
    {
    double xyz[3];
    xyz[0]=mg_noeud->get_x()+alpha*dir.get_x();
    xyz[1]=mg_noeud->get_y()+alpha*dir.get_y();
    xyz[2]=mg_noeud->get_z()+alpha*dir.get_z();
    double critcourant=1.;
    for (int i=0;i<mg_noeud->get_lien_triangle()->get_nb();i++)
        {
        MG_TRIANGLE* tri=mg_noeud->get_lien_triangle()->get(i);
        double *xyz1=tri->get_noeud1()->get_coord();
        double *xyz2=tri->get_noeud2()->get_coord();
        double *xyz3=tri->get_noeud3()->get_coord();
        if (tri->get_noeud1()==mg_noeud) xyz1=xyz;
        if (tri->get_noeud2()==mg_noeud) xyz2=xyz;
        if (tri->get_noeud3()==mg_noeud) xyz3=xyz;
        double val=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3);
        if (val<critcourant) critcourant=val;
        }
    if (critencours>critcourant)
        {
        converge=true;
        return passe;
        }
    else 
        {
        alpha=alpha+0.1;
        x=xyz[0];
        y=xyz[1];
        z=xyz[2];
        critencours=critcourant;
        crit=critencours;
        passe=1;
        }
        
    }
return 0;
}

void MAILLEUR_STL::maillage_et_segadiscretiser_update(MG_SEGMENT* segbase,double x, double y, double z, OT_VECTEUR_3D  normale,TPL_LISTE_ENTITE<MG_TRIANGLE *>& listuniquetri, std::map<unsigned long, std::pair<MG_SEGMENT *, int> >& lstsegment)
{
projete_point_sur_triangulation(x,y,z,x,y,z);
MG_NOEUD *no=mai->ajouter_mg_noeud(NULL,x,y,z,segbase->get_origine());
MG_NOEUD* nobase1=NULL;
MG_NOEUD* nobase2=NULL;
if (segbase->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)
    {
    nobase1=segbase->get_noeud1();
    nobase2=segbase->get_noeud2();
    }
for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
    if (ptmp->second.second==2)
        supprimer_seg_a_discretiser(ptmp->second.first);
for (int i=0;i<listuniquetri.get_nb();i++)
     {
     MG_TRIANGLE* tri=listuniquetri.get(i);
     mai->supprimer_mg_triangleid(tri->get_id());
     }
for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
    if (ptmp->second.second==1)
    {
    MG_SEGMENT* seg=ptmp->second.first;
    MG_SEGMENT* mgsegment1=mai->get_mg_segment(seg->get_noeud1()->get_id(),seg->get_noeud2()->get_id());
    MG_SEGMENT* mgsegment2=mai->get_mg_segment(seg->get_noeud2()->get_id(),no->get_id());
    MG_SEGMENT* mgsegment3=mai->get_mg_segment(no->get_id(),seg->get_noeud1()->get_id());
    if (mgsegment1==NULL) 
        {
        mgsegment1=mai->ajouter_mg_segment(NULL,seg->get_noeud1(),seg->get_noeud2(),MAGIC::ORIGINE::TRIANGULATION);
        ajouter_seg_a_discretiser(mgsegment1);
        }
    if (mgsegment2==NULL) 
        {
        int origine=MAGIC::ORIGINE::TRIANGULATION;
        if ((seg->get_noeud2()==nobase1)|| (seg->get_noeud2()==nobase2))
           origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
        mgsegment2=mai->ajouter_mg_segment(NULL,seg->get_noeud2(),no,origine);
        ajouter_seg_a_discretiser(mgsegment2);
        }
    if (mgsegment3==NULL) 
        {
        int origine=MAGIC::ORIGINE::TRIANGULATION;
        if ((seg->get_noeud1()==nobase1)|| (seg->get_noeud1()==nobase2))
          origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
        mgsegment3=mai->ajouter_mg_segment(NULL,no,seg->get_noeud1(),origine);
        ajouter_seg_a_discretiser(mgsegment3);
        }
    MG_TRIANGLE* nvtri=new MG_TRIANGLE(NULL,seg->get_noeud1(),seg->get_noeud2(),no,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::TRIANGULATION);
    mai->ajouter_mg_triangle(nvtri);
    double xyz[3];
    nvtri->normal(xyz);
    OT_VECTEUR_3D nvnor(xyz);
    if (nvnor*normale<0) 
       nvtri->inverse_sens();
    
    MG_TRIANGLE* tri1=get_voisin(nvtri,1);
    MG_TRIANGLE* tri2=get_voisin(nvtri,2);
    MG_TRIANGLE* tri3=get_voisin(nvtri,3);
    double nor[3],nor1[3],nor2[3],nor3[3];
    nvtri->normal(nor);
    OT_VECTEUR_3D vnor(nor);
    if (tri1!=NULL)
        {   
        tri1->normal(nor1);
        OT_VECTEUR_3D vnor1(nor1);
        if (vnor*vnor1<0.) 
            
            affiche("Triangle retourne v1");
        
        }
    if (tri2!=NULL)
        {   
        tri2->normal(nor2);
        OT_VECTEUR_3D vnor2(nor2);
        if (vnor*vnor2<0.) 
            affiche("Triangle retourne v2");
        }
    if (tri3!=NULL)
        {   
        tri3->normal(nor3);
        OT_VECTEUR_3D vnor3(nor3);
        if (vnor*vnor3<0.) 
            affiche("Triangle retourne v3");
        }
    }
}


void MAILLEUR_STL::cree_maillage()
{
int ok=0;
int compteur_iter=0;
do
{
compteur_iter++;
if (lstsegadiscretiser.begin()->first>limite_discretisation)
    {
    MG_SEGMENT* seg=lstsegadiscretiser.begin()->second;
    double x,y,z;
    OT_VECTEUR_3D normale;
    TPL_LISTE_ENTITE<MG_TRIANGLE*> listtrietoile;
    std::map<unsigned long,std::pair<MG_SEGMENT*,int> > lstsegmentetoile;
    calcul_nouveau_noeud_et_etoile(seg,x,y,z,normale,listtrietoile,lstsegmentetoile);
    maillage_et_segadiscretiser_update(seg,x,y,z,normale,listtrietoile,lstsegmentetoile);
    }
else
    ok=1;
}
while (ok==0);
char chaine[500];
sprintf(chaine,"     %d iterations effectuées",compteur_iter);
affiche(chaine);
}


bool MAILLEUR_STL::respecte_delaunay(double x, double y, double z, MG_TRIANGLE* tri)
{
OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
OT_VECTEUR_3D vec3=vec1&vec2;
vec1.norme();
vec3.norme();
vec2=vec3&vec1;
OT_MATRICE_3D base(vec1,vec2,vec3);
OT_MATRICE_3D baset=base.transpose();;
OT_VECTEUR_3D pt1(tri->get_noeud1()->get_coord());
OT_VECTEUR_3D pt2(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
OT_VECTEUR_3D pt3(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
OT_VECTEUR_3D pt1l(0.,0.,0.);
OT_VECTEUR_3D pt2l=baset*pt2;
OT_VECTEUR_3D pt3l=baset*pt3;

OT_VECTEUR_3D colx(pt1l.get_x(),pt2l.get_x(),pt3l.get_x());
OT_VECTEUR_3D coly(pt1l.get_y(),pt2l.get_y(),pt3l.get_y());
OT_VECTEUR_3D colz(1.,1.,1.);
OT_VECTEUR_3D colx2y2(pt1l.get_x()*pt1l.get_x()+pt1l.get_y()*pt1l.get_y(),pt2l.get_x()*pt2l.get_x()+pt2l.get_y()*pt2l.get_y(),pt3l.get_x()*pt3l.get_x()+pt3l.get_y()*pt3l.get_y());
OT_MATRICE_3D sys1(colx,coly,colz);
double det=2.*sys1.get_determinant();
OT_MATRICE_3D sys2(colx2y2,coly,colz);
double xcl=sys2.get_determinant()/det;
OT_MATRICE_3D sys3(colx2y2,colx,colz);
double ycl=-sys3.get_determinant()/det;
OT_VECTEUR_3D centrel(xcl,ycl,0.);
OT_VECTEUR_3D centre=pt1+base*centrel;

OT_VECTEUR_3D r1(centre.get_xyz(),tri->get_noeud1()->get_coord());
double rayon=r1.get_longueur();
OT_VECTEUR_3D dist(x,y,z);
dist=dist-centre;
bool rep=true;
if (dist.get_longueur()<rayon) 
    rep=false;
return rep;
}


void MAILLEUR_STL::projete_point_sur_triangulation(double x, double y, double z, double& nvx, double& nvy, double& nvz)
{
double mat[9];
double xyz[3]={x,y,z};
metrique->evaluer(xyz,mat);
double dref=1./sqrt(mat[0]);
TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
while (lst.get_nb()==0)
    {
    octree.rechercher(x,y,z,dref,lst);
    dref=dref*1.05;
    }
MG_MAILLAGE_OUTILS ot;
TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR it;
double dis=1e308;
for (MG_TRIANGLE* tri=lst.get_premier(it);tri!=NULL;tri=lst.get_suivant(it))
    {
    double xt,yt,zt;
    ot.projetedansletriangle(tri,x,y,z,xt,yt,zt);
    OT_VECTEUR_3D vec(x-xt,y-yt,z-zt);
    double l=vec.get_longueur();
    if (l<dis)
        {
        dis=l;
        nvx=xt;
        nvy=yt;
        nvz=zt;
        }
    }
}


void MAILLEUR_STL::change_angle_arete(double val)
{
angle_arete=val;
}


double MAILLEUR_STL::get_angle_arete()
{
return angle_arete;
}

MG_TRIANGLE * MAILLEUR_STL::get_voisin(MG_TRIANGLE* tri, int cote)
{
MG_SEGMENT* seg;
if (cote==1) seg=tri->get_segment1();
if (cote==2) seg=tri->get_segment2();
if (cote==3) seg=tri->get_segment3();
if (seg->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE) return NULL;
if (seg->get_lien_triangle()->get_nb()!=2) return NULL;
MG_TRIANGLE* tri1=seg->get_lien_triangle()->get(0);
MG_TRIANGLE* tri2=seg->get_lien_triangle()->get(1);
if (tri1==tri) return tri2;
return tri1;
}

void MAILLEUR_STL::ajouter_seg_a_discretiser(MG_SEGMENT* seg)
{
double l=metrique->calcul_distance_metrique(seg,pas);
seg->change_distance_metrique(l);
lstsegadiscretiser.insert(std::pair<double,MG_SEGMENT*>(l,seg));
}

void MAILLEUR_STL::supprimer_seg_a_discretiser(MG_SEGMENT* seg)
{
double l;
seg->get_distance_metrique(&l);
std::multimap<double,MG_SEGMENT*,std::greater<double >  >:: iterator it=lstsegadiscretiser.lower_bound(l);
while (it->second!=seg) 
    it++;
lstsegadiscretiser.erase(it);
}

void MAILLEUR_STL::change_pas(int val)
{
pas=val;
}

void MAILLEUR_STL::change_nbpasseoptimisation(int nb)
{
nbcoucheopt=nb;
}


void MAILLEUR_STL::change_nx(int val)
{
nx=val;
}


void MAILLEUR_STL::change_ny(int val)
{
ny=val;
}


void MAILLEUR_STL::change_nz(int val)
{
nz=val;
}

void MAILLEUR_STL::change_nxnynz(int nvx, int nvy, int nvz)
{
nx=nvx;
ny=nvy;
nz=nvz;
}


void MAILLEUR_STL::change_qualmaxaoptimiser(double val)
{
qualmaxaopt=val;
}
Revision:	1082
Committed:	Thu Dec 23 21:17:26 2021 UTC (3 years, 4 months ago) by francois
Original Path:	*magic/lib/mailleur_auto/src/mailleur_stl.cpp*
File size:	24246 byte(s)
Log Message:	correction de bug dans le generateur de particule RSA et DCR + ajout mailleur_stl (a ameliorer)