mtu/src/mg_arete.cpp

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mg_arete.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:55 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------


#include "gestionversion.h"
#include "mg_arete.h"
#include "vct_arete.h"
#include "mg_definition.h"
#include <math.h>
#include "ot_mathematique.h"
#include "constantegeo.h"
#include "fem_solution.h"
#include "fem_maillage.h"
#include <limits>

MG_ARETE::MG_ARETE(std::string idori,unsigned long num,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
{
}

MG_ARETE::MG_ARETE(std::string idori,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
{
}

MG_ARETE::MG_ARETE(MG_ARETE& mdd):MG_ELEMENT_TOPOLOGIQUE(mdd),cosommet1(mdd.cosommet1),cosommet2(mdd.cosommet2),courbe(mdd.courbe),orientation(mdd.orientation),vect(NULL)
{
}

MG_ARETE::MG_ARETE(std::string idori,unsigned long num,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),courbe(crb),orientation(sens),vect(NULL)
{
  cosommet1=NULL;
  cosommet2=NULL;
}
MG_ARETE::MG_ARETE(std::string idori,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),courbe(crb),orientation(sens),vect(NULL)
{
  cosommet1=NULL;
  cosommet2=NULL;
}
void MG_ARETE::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lst)
{
    MG_SOMMET* som1=cosommet1->get_sommet();
    MG_SOMMET* som2=cosommet2->get_sommet();
    lst->ajouter(som1);
    lst->ajouter(som2);
    som1->get_topologie_sousjacente(lst);
    som2->get_topologie_sousjacente(lst);
}

MG_ARETE::~MG_ARETE()
{
    // if (lst_coarete.size()!=0) afficheur << WARCOARETE << this->get_id()<< enderr;
    if (vect!=NULL) delete vect;
}

void MG_ARETE::changer_cosommet1(class MG_COSOMMET* cosom)
{
    cosommet1=cosom;
}
void MG_ARETE::changer_cosommet2(class MG_COSOMMET* cosom)
{
    cosommet2=cosom;
}
MG_COSOMMET* MG_ARETE::get_cosommet1(void)
{
    return cosommet1;
}
MG_COSOMMET*  MG_ARETE::get_cosommet2(void)
{
    return cosommet2;
}
MG_COURBE* MG_ARETE::get_courbe(void)
{
    return courbe;
}
int MG_ARETE::get_orientation(void)
{
    return orientation;
}


void  MG_ARETE::ajouter_mg_coarete(class MG_COARETE* coarete)
{
    lst_coarete.insert(lst_coarete.end(),coarete);
}

int  MG_ARETE::get_nb_mg_coarete(void)
{
    return lst_coarete.size();
}

bool MG_ARETE::est_une_topo_element(void)
{
return false;
}

BOITE_3D MG_ARETE::get_boite_3D(void)
{
  double xmin=std::numeric_limits<double>::max();
  double ymin=std::numeric_limits<double>::max();
  double zmin=std::numeric_limits<double>::max();
  double xmax=std::numeric_limits<double>::min();
  double ymax=std::numeric_limits<double>::min();
  double zmax=std::numeric_limits<double>::min();
  double tmin=get_tmin();
  double tmax=get_tmax();
  double pas_t = (tmax-tmin)/pas_echantillon;
  for(int i=0;i<pas_echantillon;i++)
  {
    double t=tmin+i*pas_t; 
    double xyz[3];
    evaluer(t,xyz);
    if (xyz[0]<xmin) xmin=xyz[0];
    if (xyz[0]>xmax) xmax=xyz[0];
    if (xyz[1]<ymin) ymin=xyz[1];
    if (xyz[1]>ymax) ymax=xyz[1];
    if (xyz[2]<zmin) zmin=xyz[2];
    if (xyz[2]>zmax) zmax=xyz[2];    
  }
  BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
  return boite;
}

void  MG_ARETE::evaluer(double t,double *xyz)
{
    if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
    courbe->evaluer(t,xyz);
}

void  MG_ARETE::deriver(double t,double *xyz)
{
    if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
    courbe->deriver(t,xyz);
    if (orientation!=MEME_SENS)
    {
        xyz[0]=-xyz[0];
        xyz[1]=-xyz[1];
        xyz[2]=-xyz[2];
    }
}

void  MG_ARETE::deriver_seconde(double t,double *ddxyz,double* dxyz,double* xyz)
{
    if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
    courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
    if (orientation!=MEME_SENS)
    {
        dxyz[0]=-dxyz[0];
        dxyz[1]=-dxyz[1];
        dxyz[2]=-dxyz[2];
    }
}

void  MG_ARETE::inverser(double& t,double *xyz,double precision)
{
    courbe->inverser(t,xyz,precision);
    if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
}

double MG_ARETE::get_tmin(void)
{
    return cosommet1->get_t();
}

double MG_ARETE::get_tmax(void)
{
    if (!courbe->est_periodique()) return cosommet2->get_t();
    double tmin=cosommet1->get_t();
    double tmax=cosommet2->get_t();
    if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
    return tmax;
}

double MG_ARETE::get_longueur(double t1,double t2,double precis)
{
    if (orientation!=MEME_SENS)
    {
        t1=courbe->get_tmin()+courbe->get_tmax()-t1;
        t2=courbe->get_tmin()+courbe->get_tmax()-t2;
        return (-courbe->get_longueur(t1,t2,precis));
    }
    return courbe->get_longueur(t1,t2,precis);

}

double MG_ARETE::get_M(double t)
{
    if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
    return courbe->get_M(t);
}


void  MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
{
    std::vector<MG_COARETE*>::iterator i;
    for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
    {
        if ((*i)==coarete)
        {
            lst_coarete.erase(i);
            return;
        }
    }
}


MG_COARETE*  MG_ARETE::get_mg_coarete(int num)
{
    return lst_coarete[num];
}

int MG_ARETE::get_type(void)
{
  return TYPE_ELEMENT_TOPOLOGIQUE::ARETE;
}


int MG_ARETE::get_dimension(void)
{
    return 1;
}


VCT& MG_ARETE::get_vectorisation(void)
{
    if (vect==NULL) vect=new VCT_ARETE(this);
    return *vect;
}

void MG_ARETE::enregistrer(std::ostream& o,double version)
{
    int nb=get_nb_ccf();
    o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," ;
    if (version<2)
    {
    o<< nb;
    if (nb!=0)
    {
        o << ",(";
        for (int i=0;i<nb;i++)
        {
            char nom[3];
            get_type_ccf(i,nom);
            o << "(" <<  nom << "," << get_valeur_ccf(i) << ")";
            if (i!=nb-1) o << "," ;
        }
        o << ")";
    }
    }
    else enregistrer_ccf(o,version);
    o <<   ");" << std::endl;
}


void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
{
    OT_VECTEUR_3D P1_INTERSECTION;

    double U,W;

    double xyz2[3];
    double xyz1[3];
    double drv1[3];
    double drv2[3];

    evaluer(t1,xyz1);
    deriver(t1,drv1);

    evaluer(t2,xyz2);
    deriver(t2,drv2);

    OT_VECTEUR_3D a(xyz1)  ;
    OT_VECTEUR_3D b(drv1)  ;
    OT_VECTEUR_3D c(xyz2)  ;
    OT_VECTEUR_3D d(drv2)  ;
    OT_VECTEUR_3D f=c-a;
    d=-d;
    double det1=b[0]*d[1]-b[1]*d[0];
    double det2=b[0]*d[2]-b[2]*d[0];
    double det3=b[1]*d[2]-b[2]*d[1];

    bool c_fait=false;

    if (fabs(det1)>=1e-10)
    {
        double u,w;
        u=(d[1]*f[0]-d[0]*f[1])/det1;
        w=(b[0]*f[1]-b[1]*f[0])/det1;
        double verif_x=a[2] +u*b[2];
        double verif_y=c[2]-w*d[2];
        if (fabs(verif_x-verif_y)<0.000001)
        {
            U=u  ;
            W=w  ;
            c_fait=true;
        }

    }

    if (!c_fait&&fabs(det2)>=1e-10)
    {
        double u,w;
        u=(d[2]*f[0]-d[0]*f[2])/det2;
        w=(b[0]*f[2]-b[2]*f[0])/det2;
        double verif_x=a[1]+u*b[1];
        double verif_y=c[1]-w*d[1];
        if (fabs(verif_x-verif_y)<0.000001)
        {
            U=u  ;
            W=w  ;
            c_fait=true;
        }
    }

    if (!c_fait&&fabs(det3)>=1e-10)
    {
        double u,w;
        u=(d[2]*f[1]-d[1]*f[2])/det3;
        w=(b[1]*f[2]-b[2]*f[1])/det3;
        double verif_x=a[0] +u*b[0];
        double verif_y=c[0]-w*d[0];
        if (fabs(verif_x-verif_y)<0.000001)
        {
            U=u  ;
            W=w  ;
            c_fait=true;
        }
    }

    P1_INTERSECTION=a+U*b;

    point_iners[0]=P1_INTERSECTION[0];
    point_iners[1]=P1_INTERSECTION[1];
    point_iners[2]=P1_INTERSECTION[2];

}


void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
{

    TPL_LISTE_ENTITE<double> param1;
    int type=courbe->get_type_geometrique(param1);

    double tcs1=get_tmin();
    double tcs2=get_tmax();
    double angle=tcs1-tcs2;

    if ((type==GEOMETRIE::CONST::Co_CIRCLE||type==GEOMETRIE::CONST::Co_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*M_PI)
    {

        double xyz[3];
        param.ajouter(1);
        param.ajouter(4);
        param.ajouter(0);
        param.ajouter(7);
        param.ajouter(0);

        double arc_lengh= get_longueur(tcs1,tcs2);
        double rayon=arc_lengh/fabs(angle);

        double dif=fabs(fabs(angle)-M_PI);
        if (fabs(fabs(angle)-M_PI)<=1e-6||fabs(angle)-M_PI<0.)  //ouverture inf ou egale a PI
        {


            double request_lengh=arc_lengh/2;

            param.ajouter(0);
            param.ajouter(0);
            param.ajouter(0);
            param.ajouter(0.5);
            param.ajouter(0.5);
            param.ajouter(1);
            param.ajouter(1);
            param.ajouter(1);


            double xyz1[3];
            double xyz1_inters[3];
            double xyz2_inters[3];
            double xyz2[3];
            double xyz3[3];
            double t_mid=get_t(tcs1,tcs2,request_lengh);

            evaluer(tcs1,xyz1);

            param.ajouter(xyz1[0]);
            param.ajouter(xyz1[1]);
            param.ajouter(xyz1[2]);
            param.ajouter(1);

            get_intersection(tcs1, t_mid,xyz1_inters);

            param.ajouter(xyz1_inters[0]);
            param.ajouter(xyz1_inters[1]);
            param.ajouter(xyz1_inters[2]);
            param.ajouter(0.5);

            get_intersection(t_mid, tcs2,xyz2_inters);

            param.ajouter(xyz2_inters[0]);
            param.ajouter(xyz2_inters[1]);
            param.ajouter(xyz2_inters[2]);
            param.ajouter(0.5);

            evaluer(tcs2,xyz3);

            param.ajouter(xyz3[0]);
            param.ajouter(xyz3[1]);
            param.ajouter(xyz3[2]);
            param.ajouter(1);

            indx_premier_ptctr=13;

        }

        if (fabs(fabs(angle)-M_PI)>1e-6&&fabs(angle)-M_PI>0.)
        {

            double xyz[3];

            double lengh_at_mid=arc_lengh;
            double request_lengh=lengh_at_mid/2.;
            double tmid=get_t(tcs1,tcs2,request_lengh);
            double tmid1=get_t(tcs1,tmid,request_lengh/2.);
            double tmid2=get_t(tmid,tcs2,request_lengh/2.);

            param.ajouter(0);
            param.ajouter(0);
            param.ajouter(0);
            param.ajouter(0.25);
            param.ajouter(0.25);
            param.ajouter(0.5);
            param.ajouter(0.5);
            param.ajouter(0.75);
            param.ajouter(0.75);
            param.ajouter(1);
            param.ajouter(1);
            param.ajouter(1);


            evaluer(tcs1,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(1);

            
            get_intersection(tcs1,tmid1,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(0.5);

            get_intersection(tmid1, tmid,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(0.5);

            evaluer(tmid,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(1);

            
            get_intersection(tmid,tmid2,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(0.5);


            get_intersection(tmid2,tcs2,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(0.5);

            evaluer(tcs2,xyz);

            param.ajouter(xyz[0]);
            param.ajouter(xyz[1]);
            param.ajouter(xyz[2]);
            param.ajouter(1);

            indx_premier_ptctr=17;
        }

    }


    if ((type==GEOMETRIE::CONST::Co_ELLIPSE||type==GEOMETRIE::CONST::Co_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*M_PI,0.000001))
        courbe->get_param_NURBS(indx_premier_ptctr,param);

    if (type==GEOMETRIE::CONST::Co_BSPLINE)
        courbe->get_param_NURBS(indx_premier_ptctr,param);

    if (type==GEOMETRIE::CONST::Co_LINE)
    {

        param.ajouter(1);
        param.ajouter(2);
        param.ajouter(0);
        param.ajouter(2);
        param.ajouter(0);


        param.ajouter(0);
        param.ajouter(0);
        param.ajouter(1);
        param.ajouter(1);

        double xyz1[3];
        double xyz2[3];

        double tcs1=get_tmin();
        double tcs2=get_tmax();

        evaluer(tcs1,xyz1);

        param.ajouter(xyz1[0]);
        param.ajouter(xyz1[1]);
        param.ajouter(xyz1[2]);
        param.ajouter(1);

        evaluer(tcs2,xyz2);

        param.ajouter(xyz2[0]);
        param.ajouter(xyz2[1]);
        param.ajouter(xyz2[2]);
        param.ajouter(1);

        indx_premier_ptctr=9;
    }

}

struct integrale {
    double ti;
    double li;
};

double  MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
{
    double pt[3];
    double t,ti,tii;

    double longueur=0;
    std::vector<integrale> tab_integrale;

    integrale pas;
    ti=t1;
    do
    {
        tii=ti+(t2-t1)*dt;

        t=0.7886751345*ti+0.2113248654*tii;
        deriver(t,pt);

        double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
        longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
        t=0.7886751345*tii+0.2113248654*ti;
        deriver(t,pt);
        facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
        longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
        integrale pas;
        pas.ti=ti;
        pas.li=longueur;
        tab_integrale.insert(tab_integrale.end(),pas);
        ti=tii;
    }
    while (fabs(tab_integrale.back().li)<fabs(longueur_voulue));

    double tn_1=tab_integrale.back().ti;
    double ln_1=tab_integrale[tab_integrale.size()-2].li;
    double tn=ti;
    double ln=tab_integrale.back().li;


    double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
    return t_trouve ;
}

void MG_ARETE::recupere_resultat(class FEM_SOLUTION* sol,int numchamps,char* fichier)
{
FEM_MAILLAGE *mai=sol->get_maillage();
MG_SOMMET* som1=cosommet1->get_sommet();
MG_SOMMET* som2=cosommet2->get_sommet();
sol->active_solution(numchamps);

FEM_NOEUD* noeudcourant=NULL;
for (int i=0;i<som1->get_lien_fem_maillage()->get_nb();i++)
    {
    FEM_NOEUD* no=(FEM_NOEUD*)som1->get_lien_fem_maillage()->get(i);
    if (mai->get_fem_noeudid(no->get_id())==no) noeudcourant=no;
    }
FEM_NOEUD* noeudfin=NULL;
for (int i=0;i<som2->get_lien_fem_maillage()->get_nb();i++)
    {
    FEM_NOEUD* no=(FEM_NOEUD*)som2->get_lien_fem_maillage()->get(i);
    if (mai->get_fem_noeudid(no->get_id())==no) noeudfin=no;
    }
int n1=0;
int n2;
if (mai->get_degre()==1) n2=1;
if (mai->get_degre()==2) n2=2;
double *xyz=noeudcourant->get_coord();
double t;
inverser(t,xyz);
std::vector<double> resultat;
double tprec=t;
resultat.push_back(t);
resultat.push_back(0.);
resultat.push_back(noeudcourant->get_solution());
while (noeudcourant!=noeudfin)
  {
  FEM_ELEMENT1 *segcourant=NULL;  
  for (int i=0;i<get_lien_fem_maillage()->get_nb();i++)
    {
    FEM_ELEMENT1* seg=(FEM_ELEMENT1*)get_lien_fem_maillage()->get(i);
    if (mai->get_fem_element1id(seg->get_id())==seg)
        {
        if ((seg->get_fem_noeud(n1)==noeudcourant)||(seg->get_fem_noeud(n2)==noeudcourant))
        if (seg!=segcourant)  
          segcourant=seg;
        }
    }
  FEM_NOEUD* autrenoeud=NULL;
  if (segcourant->get_fem_noeud(n1)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n2);
  if (segcourant->get_fem_noeud(n2)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n1);
  FEM_NOEUD *noprec=noeudcourant;
  for (int i=1;i<segcourant->get_nb_fem_noeud();i++)
      {
      FEM_NOEUD* no=segcourant->get_fem_noeud(i);
      double *xyz=no->get_coord();
      double t2;
      inverser(t2,xyz);
      if (get_courbe()->est_periodique())
        if (t2<tprec) t2=t2+get_courbe()->get_periode();
      double s=resultat[resultat.size()-2];
      int pas=32;
      for (int j=0;j<pas;j++)
          {
          double ti=tprec+1.0*j*(t2-tprec)/pas;
          double tii=tprec+1.0*(j+1)*(t2-tprec)/pas;
          double tgauss1=0.7886751345*tii+0.2113248654*ti;
          double dxyz[3];
          deriver(tgauss1,dxyz);
          s=s+0.5*(tii-ti)*sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
          double tgauss2=0.7886751345*ti+0.2113248654*tii;
          deriver(tgauss2,dxyz);
          s=s+0.5*(tii-ti)*sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
          }
      resultat.push_back(t2);
      resultat.push_back(s);
      resultat.push_back(no->get_solution());
      noprec=no;
      tprec=t2;
      }
  noeudcourant=autrenoeud;
  }
FILE *in=fopen(fichier,"wt");
char mess[500];
sprintf(mess,"%s_%s",(char*)sol->get_nom().c_str(),(char*)sol->get_legende(numchamps).c_str());
fprintf(in,"t;s;%s;\n",mess);
for (int i=0;i<resultat.size();i=i+3)
  fprintf(in,"%lf;%lf;%lf;\n",resultat[i],resultat[i+1],resultat[i+2]);
fclose(in);
}
Revision:	1180
Committed:	Fri Sep 20 20:10:22 2024 UTC (8 months ago) by francois
File size:	18329 byte(s)
Log Message:	Correction dans la quadratisation sans geometrie
#	User	Rev	Content
1	francois	1158	//####//------------------------------------------------------------
2			//####//------------------------------------------------------------
3			//####// MAGiC
4			//####// Jean Christophe Cuilliere et Vincent FRANCOIS
5			//####// Departement de Genie Mecanique - UQTR
6			//####//------------------------------------------------------------
7			//####// MAGIC est un projet de recherche de l equipe ERICCA
8			//####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9			//####// http://www.uqtr.ca/ericca
10			//####// http://www.uqtr.ca/
11			//####//------------------------------------------------------------
12			//####//------------------------------------------------------------
13			//####//
14			//####// mg_arete.cpp
15			//####//
16			//####//------------------------------------------------------------
17			//####//------------------------------------------------------------
18			//####// COPYRIGHT 2000-2024
19			//####// jeu 13 jun 2024 11:58:55 EDT
20			//####//------------------------------------------------------------
21			//####//------------------------------------------------------------
22	francois	283
23
24			#include "gestionversion.h"
25			#include "mg_arete.h"
26			#include "vct_arete.h"
27	francois	375	#include "mg_definition.h"
28	francois	283	#include <math.h>
29			#include "ot_mathematique.h"
30			#include "constantegeo.h"
31	francois	604	#include "fem_solution.h"
32			#include "fem_maillage.h"
33	francois	1158	#include <limits>
34	francois	283
35			MG_ARETE::MG_ARETE(std::string idori,unsigned long num,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
36			{
37			}
38
39			MG_ARETE::MG_ARETE(std::string idori,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
40			{
41			}
42
43			MG_ARETE::MG_ARETE(MG_ARETE& mdd):MG_ELEMENT_TOPOLOGIQUE(mdd),cosommet1(mdd.cosommet1),cosommet2(mdd.cosommet2),courbe(mdd.courbe),orientation(mdd.orientation),vect(NULL)
44			{
45			}
46
47			MG_ARETE::MG_ARETE(std::string idori,unsigned long num,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),courbe(crb),orientation(sens),vect(NULL)
48			{
49	couturad	906	cosommet1=NULL;
50			cosommet2=NULL;
51	francois	283	}
52			MG_ARETE::MG_ARETE(std::string idori,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),courbe(crb),orientation(sens),vect(NULL)
53			{
54	couturad	906	cosommet1=NULL;
55			cosommet2=NULL;
56	francois	283	}
57			void MG_ARETE::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE> lst)
58			{
59			MG_SOMMET* som1=cosommet1->get_sommet();
60			MG_SOMMET* som2=cosommet2->get_sommet();
61			lst->ajouter(som1);
62			lst->ajouter(som2);
63			som1->get_topologie_sousjacente(lst);
64			som2->get_topologie_sousjacente(lst);
65			}
66
67			MG_ARETE::~MG_ARETE()
68			{
69			// if (lst_coarete.size()!=0) afficheur << WARCOARETE << this->get_id()<< enderr;
70			if (vect!=NULL) delete vect;
71			}
72
73			void MG_ARETE::changer_cosommet1(class MG_COSOMMET* cosom)
74			{
75			cosommet1=cosom;
76			}
77			void MG_ARETE::changer_cosommet2(class MG_COSOMMET* cosom)
78			{
79			cosommet2=cosom;
80			}
81			MG_COSOMMET* MG_ARETE::get_cosommet1(void)
82			{
83			return cosommet1;
84			}
85			MG_COSOMMET* MG_ARETE::get_cosommet2(void)
86			{
87			return cosommet2;
88			}
89			MG_COURBE* MG_ARETE::get_courbe(void)
90			{
91			return courbe;
92			}
93			int MG_ARETE::get_orientation(void)
94			{
95			return orientation;
96			}
97
98
99
100
101			void MG_ARETE::ajouter_mg_coarete(class MG_COARETE* coarete)
102			{
103			lst_coarete.insert(lst_coarete.end(),coarete);
104			}
105
106			int MG_ARETE::get_nb_mg_coarete(void)
107			{
108			return lst_coarete.size();
109			}
110
111	francois	1180	bool MG_ARETE::est_une_topo_element(void)
112	francois	576	{
113			return false;
114			}
115	francois	283
116	francois	1095	BOITE_3D MG_ARETE::get_boite_3D(void)
117	couturad	911	{
118			double xmin=std::numeric_limits<double>::max();
119			double ymin=std::numeric_limits<double>::max();
120			double zmin=std::numeric_limits<double>::max();
121			double xmax=std::numeric_limits<double>::min();
122			double ymax=std::numeric_limits<double>::min();
123			double zmax=std::numeric_limits<double>::min();
124			double tmin=get_tmin();
125			double tmax=get_tmax();
126			double pas_t = (tmax-tmin)/pas_echantillon;
127			for(int i=0;i<pas_echantillon;i++)
128			{
129			double t=tmin+i*pas_t;
130			double xyz[3];
131			evaluer(t,xyz);
132			if (xyz[0]<xmin) xmin=xyz[0];
133			if (xyz[0]>xmax) xmax=xyz[0];
134			if (xyz[1]<ymin) ymin=xyz[1];
135			if (xyz[1]>ymax) ymax=xyz[1];
136			if (xyz[2]<zmin) zmin=xyz[2];
137			if (xyz[2]>zmax) zmax=xyz[2];
138			}
139			BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
140			return boite;
141			}
142
143	francois	283	void MG_ARETE::evaluer(double t,double *xyz)
144			{
145			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
146			courbe->evaluer(t,xyz);
147			}
148
149			void MG_ARETE::deriver(double t,double *xyz)
150			{
151			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
152			courbe->deriver(t,xyz);
153			if (orientation!=MEME_SENS)
154			{
155			xyz[0]=-xyz[0];
156			xyz[1]=-xyz[1];
157			xyz[2]=-xyz[2];
158			}
159			}
160
161			void MG_ARETE::deriver_seconde(double t,double ddxyz,double dxyz,double* xyz)
162			{
163			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
164			courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
165			if (orientation!=MEME_SENS)
166			{
167			dxyz[0]=-dxyz[0];
168			dxyz[1]=-dxyz[1];
169			dxyz[2]=-dxyz[2];
170			}
171			}
172
173			void MG_ARETE::inverser(double& t,double *xyz,double precision)
174			{
175			courbe->inverser(t,xyz,precision);
176			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
177			}
178
179			double MG_ARETE::get_tmin(void)
180			{
181			return cosommet1->get_t();
182			}
183
184			double MG_ARETE::get_tmax(void)
185			{
186			if (!courbe->est_periodique()) return cosommet2->get_t();
187			double tmin=cosommet1->get_t();
188			double tmax=cosommet2->get_t();
189			if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
190			return tmax;
191			}
192
193			double MG_ARETE::get_longueur(double t1,double t2,double precis)
194			{
195			if (orientation!=MEME_SENS)
196			{
197			t1=courbe->get_tmin()+courbe->get_tmax()-t1;
198			t2=courbe->get_tmin()+courbe->get_tmax()-t2;
199			return (-courbe->get_longueur(t1,t2,precis));
200			}
201			return courbe->get_longueur(t1,t2,precis);
202
203			}
204
205			double MG_ARETE::get_M(double t)
206			{
207			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
208			return courbe->get_M(t);
209			}
210
211
212
213			void MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
214			{
215			std::vector<MG_COARETE*>::iterator i;
216			for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
217			{
218			if ((*i)==coarete)
219			{
220			lst_coarete.erase(i);
221			return;
222			}
223			}
224			}
225
226
227
228			MG_COARETE* MG_ARETE::get_mg_coarete(int num)
229			{
230			return lst_coarete[num];
231			}
232
233	couturad	906	int MG_ARETE::get_type(void)
234			{
235			return TYPE_ELEMENT_TOPOLOGIQUE::ARETE;
236			}
237	francois	283
238	couturad	906
239	francois	283	int MG_ARETE::get_dimension(void)
240			{
241			return 1;
242			}
243
244
245			VCT& MG_ARETE::get_vectorisation(void)
246			{
247			if (vect==NULL) vect=new VCT_ARETE(this);
248			return *vect;
249			}
250
251	francois	763	void MG_ARETE::enregistrer(std::ostream& o,double version)
252	francois	283	{
253			int nb=get_nb_ccf();
254	francois	763	o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," ;
255			if (version<2)
256			{
257			o<< nb;
258	francois	283	if (nb!=0)
259			{
260			o << ",(";
261			for (int i=0;i<nb;i++)
262			{
263			char nom[3];
264			get_type_ccf(i,nom);
265			o << "(" << nom << "," << get_valeur_ccf(i) << ")";
266			if (i!=nb-1) o << "," ;
267			}
268			o << ")";
269			}
270	francois	763	}
271			else enregistrer_ccf(o,version);
272	francois	283	o << ");" << std::endl;
273			}
274
275
276			void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
277			{
278			OT_VECTEUR_3D P1_INTERSECTION;
279
280			double U,W;
281
282			double xyz2[3];
283			double xyz1[3];
284			double drv1[3];
285			double drv2[3];
286
287			evaluer(t1,xyz1);
288			deriver(t1,drv1);
289
290			evaluer(t2,xyz2);
291			deriver(t2,drv2);
292
293			OT_VECTEUR_3D a(xyz1) ;
294			OT_VECTEUR_3D b(drv1) ;
295			OT_VECTEUR_3D c(xyz2) ;
296			OT_VECTEUR_3D d(drv2) ;
297			OT_VECTEUR_3D f=c-a;
298			d=-d;
299			double det1=b[0]d[1]-b[1]d[0];
300			double det2=b[0]d[2]-b[2]d[0];
301			double det3=b[1]d[2]-b[2]d[1];
302
303			bool c_fait=false;
304
305			if (fabs(det1)>=1e-10)
306			{
307			double u,w;
308			u=(d[1]f[0]-d[0]f[1])/det1;
309			w=(b[0]f[1]-b[1]f[0])/det1;
310			double verif_x=a[2] +u*b[2];
311			double verif_y=c[2]-w*d[2];
312			if (fabs(verif_x-verif_y)<0.000001)
313			{
314			U=u ;
315			W=w ;
316			c_fait=true;
317			}
318
319			}
320
321			if (!c_fait&&fabs(det2)>=1e-10)
322			{
323			double u,w;
324			u=(d[2]f[0]-d[0]f[2])/det2;
325			w=(b[0]f[2]-b[2]f[0])/det2;
326			double verif_x=a[1]+u*b[1];
327			double verif_y=c[1]-w*d[1];
328			if (fabs(verif_x-verif_y)<0.000001)
329			{
330			U=u ;
331			W=w ;
332			c_fait=true;
333			}
334			}
335
336			if (!c_fait&&fabs(det3)>=1e-10)
337			{
338			double u,w;
339			u=(d[2]f[1]-d[1]f[2])/det3;
340			w=(b[1]f[2]-b[2]f[1])/det3;
341			double verif_x=a[0] +u*b[0];
342			double verif_y=c[0]-w*d[0];
343			if (fabs(verif_x-verif_y)<0.000001)
344			{
345			U=u ;
346			W=w ;
347			c_fait=true;
348			}
349			}
350
351			P1_INTERSECTION=a+U*b;
352
353			point_iners[0]=P1_INTERSECTION[0];
354			point_iners[1]=P1_INTERSECTION[1];
355			point_iners[2]=P1_INTERSECTION[2];
356
357			}
358
359
360
361			void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
362			{
363
364			TPL_LISTE_ENTITE<double> param1;
365			int type=courbe->get_type_geometrique(param1);
366
367			double tcs1=get_tmin();
368			double tcs2=get_tmax();
369			double angle=tcs1-tcs2;
370
371	francois	1158	if ((type==GEOMETRIE::CONST::Co_CIRCLE\|\|type==GEOMETRIE::CONST::Co_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*M_PI)
372	francois	283	{
373
374			double xyz[3];
375			param.ajouter(1);
376			param.ajouter(4);
377			param.ajouter(0);
378			param.ajouter(7);
379			param.ajouter(0);
380
381			double arc_lengh= get_longueur(tcs1,tcs2);
382			double rayon=arc_lengh/fabs(angle);
383
384	francois	1158	double dif=fabs(fabs(angle)-M_PI);
385			if (fabs(fabs(angle)-M_PI)<=1e-6\|\|fabs(angle)-M_PI<0.) //ouverture inf ou egale a PI
386	francois	283	{
387
388
389			double request_lengh=arc_lengh/2;
390
391			param.ajouter(0);
392			param.ajouter(0);
393			param.ajouter(0);
394			param.ajouter(0.5);
395			param.ajouter(0.5);
396			param.ajouter(1);
397			param.ajouter(1);
398			param.ajouter(1);
399
400
401			double xyz1[3];
402			double xyz1_inters[3];
403			double xyz2_inters[3];
404			double xyz2[3];
405			double xyz3[3];
406			double t_mid=get_t(tcs1,tcs2,request_lengh);
407
408			evaluer(tcs1,xyz1);
409
410			param.ajouter(xyz1[0]);
411			param.ajouter(xyz1[1]);
412			param.ajouter(xyz1[2]);
413			param.ajouter(1);
414
415			get_intersection(tcs1, t_mid,xyz1_inters);
416
417			param.ajouter(xyz1_inters[0]);
418			param.ajouter(xyz1_inters[1]);
419			param.ajouter(xyz1_inters[2]);
420			param.ajouter(0.5);
421
422			get_intersection(t_mid, tcs2,xyz2_inters);
423
424			param.ajouter(xyz2_inters[0]);
425			param.ajouter(xyz2_inters[1]);
426			param.ajouter(xyz2_inters[2]);
427			param.ajouter(0.5);
428
429			evaluer(tcs2,xyz3);
430
431			param.ajouter(xyz3[0]);
432			param.ajouter(xyz3[1]);
433			param.ajouter(xyz3[2]);
434			param.ajouter(1);
435
436			indx_premier_ptctr=13;
437
438			}
439
440	francois	1158	if (fabs(fabs(angle)-M_PI)>1e-6&&fabs(angle)-M_PI>0.)
441	francois	283	{
442
443			double xyz[3];
444
445	francois	290	double lengh_at_mid=arc_lengh;
446			double request_lengh=lengh_at_mid/2.;
447			double tmid=get_t(tcs1,tcs2,request_lengh);
448			double tmid1=get_t(tcs1,tmid,request_lengh/2.);
449			double tmid2=get_t(tmid,tcs2,request_lengh/2.);
450	francois	283
451			param.ajouter(0);
452			param.ajouter(0);
453			param.ajouter(0);
454			param.ajouter(0.25);
455			param.ajouter(0.25);
456			param.ajouter(0.5);
457			param.ajouter(0.5);
458			param.ajouter(0.75);
459			param.ajouter(0.75);
460			param.ajouter(1);
461			param.ajouter(1);
462			param.ajouter(1);
463
464
465
466			evaluer(tcs1,xyz);
467
468			param.ajouter(xyz[0]);
469			param.ajouter(xyz[1]);
470			param.ajouter(xyz[2]);
471			param.ajouter(1);
472
473	francois	290
474			get_intersection(tcs1,tmid1,xyz);
475	francois	283
476			param.ajouter(xyz[0]);
477			param.ajouter(xyz[1]);
478			param.ajouter(xyz[2]);
479			param.ajouter(0.5);
480
481	francois	290	get_intersection(tmid1, tmid,xyz);
482	francois	283
483			param.ajouter(xyz[0]);
484			param.ajouter(xyz[1]);
485			param.ajouter(xyz[2]);
486			param.ajouter(0.5);
487
488	francois	290	evaluer(tmid,xyz);
489	francois	283
490			param.ajouter(xyz[0]);
491			param.ajouter(xyz[1]);
492			param.ajouter(xyz[2]);
493			param.ajouter(1);
494
495	francois	290
496			get_intersection(tmid,tmid2,xyz);
497	francois	283
498			param.ajouter(xyz[0]);
499			param.ajouter(xyz[1]);
500			param.ajouter(xyz[2]);
501			param.ajouter(0.5);
502
503
504	francois	290	get_intersection(tmid2,tcs2,xyz);
505	francois	283
506			param.ajouter(xyz[0]);
507			param.ajouter(xyz[1]);
508			param.ajouter(xyz[2]);
509			param.ajouter(0.5);
510
511			evaluer(tcs2,xyz);
512
513			param.ajouter(xyz[0]);
514			param.ajouter(xyz[1]);
515			param.ajouter(xyz[2]);
516			param.ajouter(1);
517
518			indx_premier_ptctr=17;
519			}
520
521			}
522
523
524	francois	1158	if ((type==GEOMETRIE::CONST::Co_ELLIPSE\|\|type==GEOMETRIE::CONST::Co_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*M_PI,0.000001))
525	francois	283	courbe->get_param_NURBS(indx_premier_ptctr,param);
526
527	francois	1149	if (type==GEOMETRIE::CONST::Co_BSPLINE)
528	francois	283	courbe->get_param_NURBS(indx_premier_ptctr,param);
529
530	francois	1149	if (type==GEOMETRIE::CONST::Co_LINE)
531	francois	283	{
532
533			param.ajouter(1);
534			param.ajouter(2);
535			param.ajouter(0);
536			param.ajouter(2);
537			param.ajouter(0);
538
539
540			param.ajouter(0);
541			param.ajouter(0);
542			param.ajouter(1);
543			param.ajouter(1);
544
545			double xyz1[3];
546			double xyz2[3];
547
548			double tcs1=get_tmin();
549			double tcs2=get_tmax();
550
551			evaluer(tcs1,xyz1);
552
553			param.ajouter(xyz1[0]);
554			param.ajouter(xyz1[1]);
555			param.ajouter(xyz1[2]);
556			param.ajouter(1);
557
558			evaluer(tcs2,xyz2);
559
560			param.ajouter(xyz2[0]);
561			param.ajouter(xyz2[1]);
562			param.ajouter(xyz2[2]);
563			param.ajouter(1);
564
565			indx_premier_ptctr=9;
566			}
567
568			}
569
570			struct integrale {
571			double ti;
572			double li;
573			};
574
575			double MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
576			{
577			double pt[3];
578			double t,ti,tii;
579
580			double longueur=0;
581			std::vector<integrale> tab_integrale;
582
583			integrale pas;
584			ti=t1;
585			do
586			{
587			tii=ti+(t2-t1)*dt;
588
589			t=0.7886751345ti+0.2113248654tii;
590			deriver(t,pt);
591
592			double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
593			longueur=longueur+0.5(tii-ti)sqrt(facteur);
594			t=0.7886751345tii+0.2113248654ti;
595			deriver(t,pt);
596			facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
597			longueur=longueur+0.5(tii-ti)sqrt(facteur);
598			integrale pas;
599			pas.ti=ti;
600			pas.li=longueur;
601			tab_integrale.insert(tab_integrale.end(),pas);
602			ti=tii;
603			}
604			while (fabs(tab_integrale.back().li)<fabs(longueur_voulue));
605
606			double tn_1=tab_integrale.back().ti;
607			double ln_1=tab_integrale[tab_integrale.size()-2].li;
608			double tn=ti;
609			double ln=tab_integrale.back().li;
610
611
612			double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
613			return t_trouve ;
614			}
615
616	francois	604	void MG_ARETE::recupere_resultat(class FEM_SOLUTION* sol,int numchamps,char* fichier)
617			{
618			FEM_MAILLAGE *mai=sol->get_maillage();
619			MG_SOMMET* som1=cosommet1->get_sommet();
620			MG_SOMMET* som2=cosommet2->get_sommet();
621			sol->active_solution(numchamps);
622
623			FEM_NOEUD* noeudcourant=NULL;
624			for (int i=0;i<som1->get_lien_fem_maillage()->get_nb();i++)
625			{
626			FEM_NOEUD* no=(FEM_NOEUD*)som1->get_lien_fem_maillage()->get(i);
627			if (mai->get_fem_noeudid(no->get_id())==no) noeudcourant=no;
628			}
629			FEM_NOEUD* noeudfin=NULL;
630			for (int i=0;i<som2->get_lien_fem_maillage()->get_nb();i++)
631			{
632			FEM_NOEUD* no=(FEM_NOEUD*)som2->get_lien_fem_maillage()->get(i);
633			if (mai->get_fem_noeudid(no->get_id())==no) noeudfin=no;
634			}
635			int n1=0;
636			int n2;
637			if (mai->get_degre()==1) n2=1;
638			if (mai->get_degre()==2) n2=2;
639			double *xyz=noeudcourant->get_coord();
640			double t;
641			inverser(t,xyz);
642			std::vector<double> resultat;
643			double tprec=t;
644			resultat.push_back(t);
645			resultat.push_back(0.);
646			resultat.push_back(noeudcourant->get_solution());
647			while (noeudcourant!=noeudfin)
648			{
649			FEM_ELEMENT1 *segcourant=NULL;
650			for (int i=0;i<get_lien_fem_maillage()->get_nb();i++)
651			{
652			FEM_ELEMENT1* seg=(FEM_ELEMENT1*)get_lien_fem_maillage()->get(i);
653			if (mai->get_fem_element1id(seg->get_id())==seg)
654			{
655			if ((seg->get_fem_noeud(n1)==noeudcourant)\|\|(seg->get_fem_noeud(n2)==noeudcourant))
656			if (seg!=segcourant)
657			segcourant=seg;
658			}
659			}
660			FEM_NOEUD* autrenoeud=NULL;
661			if (segcourant->get_fem_noeud(n1)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n2);
662			if (segcourant->get_fem_noeud(n2)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n1);
663			FEM_NOEUD *noprec=noeudcourant;
664			for (int i=1;i<segcourant->get_nb_fem_noeud();i++)
665			{
666			FEM_NOEUD* no=segcourant->get_fem_noeud(i);
667			double *xyz=no->get_coord();
668			double t2;
669			inverser(t2,xyz);
670	francois	672	if (get_courbe()->est_periodique())
671			if (t2<tprec) t2=t2+get_courbe()->get_periode();
672	francois	604	double s=resultat[resultat.size()-2];
673			int pas=32;
674			for (int j=0;j<pas;j++)
675			{
676			double ti=tprec+1.0j(t2-tprec)/pas;
677			double tii=tprec+1.0(j+1)(t2-tprec)/pas;
678			double tgauss1=0.7886751345tii+0.2113248654ti;
679			double dxyz[3];
680			deriver(tgauss1,dxyz);
681			s=s+0.5(tii-ti)sqrt(dxyz[0]dxyz[0]+dxyz[1]dxyz[1]+dxyz[2]*dxyz[2]);
682			double tgauss2=0.7886751345ti+0.2113248654tii;
683			deriver(tgauss2,dxyz);
684			s=s+0.5(tii-ti)sqrt(dxyz[0]dxyz[0]+dxyz[1]dxyz[1]+dxyz[2]*dxyz[2]);
685			}
686			resultat.push_back(t2);
687			resultat.push_back(s);
688			resultat.push_back(no->get_solution());
689			noprec=no;
690			tprec=t2;
691			}
692			noeudcourant=autrenoeud;
693			}
694			FILE *in=fopen(fichier,"wt");
695			char mess[500];
696			sprintf(mess,"%s_%s",(char)sol->get_nom().c_str(),(char)sol->get_legende(numchamps).c_str());
697			fprintf(in,"t;s;%s;\n",mess);
698			for (int i=0;i<resultat.size();i=i+3)
699			fprintf(in,"%lf;%lf;%lf;\n",resultat[i],resultat[i+1],resultat[i+2]);
700			fclose(in);
701	francois	1095	}