geometrie/src/mg_arete.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)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       mg_arete.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 � 11H22
//------------------------------------------------------------
//------------------------------------------------------------


#include "gestionversion.h"
#include "mg_arete.h"
#include "vct_arete.h"
//#include "message.h"
//#include "affiche.h"
#include "mg_definition.h"
#include <math.h>
#include "ot_mathematique.h"
#include "constantegeo.h"
#include "fem_solution.h"
#include "fem_maillage.h"
#define PI 3.1415926535897932384626433832795

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_arete_element(void)
{
return false;
}

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==MGCo_CIRCLE||type==MGCo_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*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)-PI);
        if (fabs(fabs(angle)-PI)<=1e-6||fabs(angle)-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)-PI)>1e-6&&fabs(angle)-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==MGCo_ELLIPSE||type==MGCo_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*PI,0.000001))
        courbe->get_param_NURBS(indx_premier_ptctr,param);

    if (type==MGCo_BSPLINE)
        courbe->get_param_NURBS(indx_premier_ptctr,param);

    if (type==MGCo_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_voulue,longueur=0;
    double longueur=0;
    std::vector<integrale> tab_integrale;
//longueur_voulue= alpha *get_longueur(t1,t2);

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