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;
}

BOITE_3D MG_ARETE::get_boite_3D(int pas_echantillon)
{
  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==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:	911
Committed:	Wed Jan 10 17:25:24 2018 UTC (7 years, 4 months ago) by couturad
File size:	18436 byte(s)
Log Message:	* Modifications pour MG_CG * Modification de MG_SOUS_GEOMETRIE pour la fonction de compression
#	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	couturad	911	BOITE_3D MG_ARETE::get_boite_3D(int pas_echantillon)
120			{
121			double xmin=std::numeric_limits<double>::max();
122			double ymin=std::numeric_limits<double>::max();
123			double zmin=std::numeric_limits<double>::max();
124			double xmax=std::numeric_limits<double>::min();
125			double ymax=std::numeric_limits<double>::min();
126			double zmax=std::numeric_limits<double>::min();
127			double tmin=get_tmin();
128			double tmax=get_tmax();
129			double pas_t = (tmax-tmin)/pas_echantillon;
130			for(int i=0;i<pas_echantillon;i++)
131			{
132			double t=tmin+i*pas_t;
133			double xyz[3];
134			evaluer(t,xyz);
135			if (xyz[0]<xmin) xmin=xyz[0];
136			if (xyz[0]>xmax) xmax=xyz[0];
137			if (xyz[1]<ymin) ymin=xyz[1];
138			if (xyz[1]>ymax) ymax=xyz[1];
139			if (xyz[2]<zmin) zmin=xyz[2];
140			if (xyz[2]>zmax) zmax=xyz[2];
141			}
142			BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
143			return boite;
144			}
145
146	francois	283	void MG_ARETE::evaluer(double t,double *xyz)
147			{
148			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
149			courbe->evaluer(t,xyz);
150			}
151
152			void MG_ARETE::deriver(double t,double *xyz)
153			{
154			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
155			courbe->deriver(t,xyz);
156			if (orientation!=MEME_SENS)
157			{
158			xyz[0]=-xyz[0];
159			xyz[1]=-xyz[1];
160			xyz[2]=-xyz[2];
161			}
162			}
163
164			void MG_ARETE::deriver_seconde(double t,double ddxyz,double dxyz,double* xyz)
165			{
166			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
167			courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
168			if (orientation!=MEME_SENS)
169			{
170			dxyz[0]=-dxyz[0];
171			dxyz[1]=-dxyz[1];
172			dxyz[2]=-dxyz[2];
173			}
174			}
175
176			void MG_ARETE::inverser(double& t,double *xyz,double precision)
177			{
178			courbe->inverser(t,xyz,precision);
179			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
180			}
181
182			double MG_ARETE::get_tmin(void)
183			{
184			return cosommet1->get_t();
185			}
186
187			double MG_ARETE::get_tmax(void)
188			{
189			if (!courbe->est_periodique()) return cosommet2->get_t();
190			double tmin=cosommet1->get_t();
191			double tmax=cosommet2->get_t();
192			if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
193			return tmax;
194			}
195
196			double MG_ARETE::get_longueur(double t1,double t2,double precis)
197			{
198			if (orientation!=MEME_SENS)
199			{
200			t1=courbe->get_tmin()+courbe->get_tmax()-t1;
201			t2=courbe->get_tmin()+courbe->get_tmax()-t2;
202			return (-courbe->get_longueur(t1,t2,precis));
203			}
204			return courbe->get_longueur(t1,t2,precis);
205
206			}
207
208			double MG_ARETE::get_M(double t)
209			{
210			if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
211			return courbe->get_M(t);
212			}
213
214
215
216			void MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
217			{
218			std::vector<MG_COARETE*>::iterator i;
219			for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
220			{
221			if ((*i)==coarete)
222			{
223			lst_coarete.erase(i);
224			return;
225			}
226			}
227			}
228
229
230
231			MG_COARETE* MG_ARETE::get_mg_coarete(int num)
232			{
233			return lst_coarete[num];
234			}
235
236	couturad	906	int MG_ARETE::get_type(void)
237			{
238			return TYPE_ELEMENT_TOPOLOGIQUE::ARETE;
239			}
240	francois	283
241	couturad	906
242	francois	283	int MG_ARETE::get_dimension(void)
243			{
244			return 1;
245			}
246
247
248			VCT& MG_ARETE::get_vectorisation(void)
249			{
250			if (vect==NULL) vect=new VCT_ARETE(this);
251			return *vect;
252			}
253
254	francois	763	void MG_ARETE::enregistrer(std::ostream& o,double version)
255	francois	283	{
256			int nb=get_nb_ccf();
257	francois	763	o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," ;
258			if (version<2)
259			{
260			o<< nb;
261	francois	283	if (nb!=0)
262			{
263			o << ",(";
264			for (int i=0;i<nb;i++)
265			{
266			char nom[3];
267			get_type_ccf(i,nom);
268			o << "(" << nom << "," << get_valeur_ccf(i) << ")";
269			if (i!=nb-1) o << "," ;
270			}
271			o << ")";
272			}
273	francois	763	}
274			else enregistrer_ccf(o,version);
275	francois	283	o << ");" << std::endl;
276			}
277
278
279			void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
280			{
281			OT_VECTEUR_3D P1_INTERSECTION;
282
283			double U,W;
284
285			double xyz2[3];
286			double xyz1[3];
287			double drv1[3];
288			double drv2[3];
289
290			evaluer(t1,xyz1);
291			deriver(t1,drv1);
292
293			evaluer(t2,xyz2);
294			deriver(t2,drv2);
295
296			OT_VECTEUR_3D a(xyz1) ;
297			OT_VECTEUR_3D b(drv1) ;
298			OT_VECTEUR_3D c(xyz2) ;
299			OT_VECTEUR_3D d(drv2) ;
300			OT_VECTEUR_3D f=c-a;
301			d=-d;
302			double det1=b[0]d[1]-b[1]d[0];
303			double det2=b[0]d[2]-b[2]d[0];
304			double det3=b[1]d[2]-b[2]d[1];
305
306			bool c_fait=false;
307
308			if (fabs(det1)>=1e-10)
309			{
310			double u,w;
311			u=(d[1]f[0]-d[0]f[1])/det1;
312			w=(b[0]f[1]-b[1]f[0])/det1;
313			double verif_x=a[2] +u*b[2];
314			double verif_y=c[2]-w*d[2];
315			if (fabs(verif_x-verif_y)<0.000001)
316			{
317			U=u ;
318			W=w ;
319			c_fait=true;
320			}
321
322			}
323
324			if (!c_fait&&fabs(det2)>=1e-10)
325			{
326			double u,w;
327			u=(d[2]f[0]-d[0]f[2])/det2;
328			w=(b[0]f[2]-b[2]f[0])/det2;
329			double verif_x=a[1]+u*b[1];
330			double verif_y=c[1]-w*d[1];
331			if (fabs(verif_x-verif_y)<0.000001)
332			{
333			U=u ;
334			W=w ;
335			c_fait=true;
336			}
337			}
338
339			if (!c_fait&&fabs(det3)>=1e-10)
340			{
341			double u,w;
342			u=(d[2]f[1]-d[1]f[2])/det3;
343			w=(b[1]f[2]-b[2]f[1])/det3;
344			double verif_x=a[0] +u*b[0];
345			double verif_y=c[0]-w*d[0];
346			if (fabs(verif_x-verif_y)<0.000001)
347			{
348			U=u ;
349			W=w ;
350			c_fait=true;
351			}
352			}
353
354			P1_INTERSECTION=a+U*b;
355
356			point_iners[0]=P1_INTERSECTION[0];
357			point_iners[1]=P1_INTERSECTION[1];
358			point_iners[2]=P1_INTERSECTION[2];
359
360			}
361
362
363
364			void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
365			{
366
367			TPL_LISTE_ENTITE<double> param1;
368			int type=courbe->get_type_geometrique(param1);
369
370			double tcs1=get_tmin();
371			double tcs2=get_tmax();
372			double angle=tcs1-tcs2;
373
374			if ((type==MGCo_CIRCLE\|\|type==MGCo_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*PI)
375			{
376
377			double xyz[3];
378			param.ajouter(1);
379			param.ajouter(4);
380			param.ajouter(0);
381			param.ajouter(7);
382			param.ajouter(0);
383
384			double arc_lengh= get_longueur(tcs1,tcs2);
385			double rayon=arc_lengh/fabs(angle);
386
387			double dif=fabs(fabs(angle)-PI);
388			if (fabs(fabs(angle)-PI)<=1e-6\|\|fabs(angle)-PI<0.) //ouverture inf ou egale a PI
389			{
390
391
392			double request_lengh=arc_lengh/2;
393
394			param.ajouter(0);
395			param.ajouter(0);
396			param.ajouter(0);
397			param.ajouter(0.5);
398			param.ajouter(0.5);
399			param.ajouter(1);
400			param.ajouter(1);
401			param.ajouter(1);
402
403
404			double xyz1[3];
405			double xyz1_inters[3];
406			double xyz2_inters[3];
407			double xyz2[3];
408			double xyz3[3];
409			double t_mid=get_t(tcs1,tcs2,request_lengh);
410
411			evaluer(tcs1,xyz1);
412
413			param.ajouter(xyz1[0]);
414			param.ajouter(xyz1[1]);
415			param.ajouter(xyz1[2]);
416			param.ajouter(1);
417
418			get_intersection(tcs1, t_mid,xyz1_inters);
419
420			param.ajouter(xyz1_inters[0]);
421			param.ajouter(xyz1_inters[1]);
422			param.ajouter(xyz1_inters[2]);
423			param.ajouter(0.5);
424
425			get_intersection(t_mid, tcs2,xyz2_inters);
426
427			param.ajouter(xyz2_inters[0]);
428			param.ajouter(xyz2_inters[1]);
429			param.ajouter(xyz2_inters[2]);
430			param.ajouter(0.5);
431
432			evaluer(tcs2,xyz3);
433
434			param.ajouter(xyz3[0]);
435			param.ajouter(xyz3[1]);
436			param.ajouter(xyz3[2]);
437			param.ajouter(1);
438
439			indx_premier_ptctr=13;
440
441			}
442
443			if (fabs(fabs(angle)-PI)>1e-6&&fabs(angle)-PI>0.)
444			{
445
446			double xyz[3];
447
448	francois	290	double lengh_at_mid=arc_lengh;
449			double request_lengh=lengh_at_mid/2.;
450			double tmid=get_t(tcs1,tcs2,request_lengh);
451			double tmid1=get_t(tcs1,tmid,request_lengh/2.);
452			double tmid2=get_t(tmid,tcs2,request_lengh/2.);
453	francois	283
454			param.ajouter(0);
455			param.ajouter(0);
456			param.ajouter(0);
457			param.ajouter(0.25);
458			param.ajouter(0.25);
459			param.ajouter(0.5);
460			param.ajouter(0.5);
461			param.ajouter(0.75);
462			param.ajouter(0.75);
463			param.ajouter(1);
464			param.ajouter(1);
465			param.ajouter(1);
466
467
468
469			evaluer(tcs1,xyz);
470
471			param.ajouter(xyz[0]);
472			param.ajouter(xyz[1]);
473			param.ajouter(xyz[2]);
474			param.ajouter(1);
475
476	francois	290
477			get_intersection(tcs1,tmid1,xyz);
478	francois	283
479			param.ajouter(xyz[0]);
480			param.ajouter(xyz[1]);
481			param.ajouter(xyz[2]);
482			param.ajouter(0.5);
483
484	francois	290	get_intersection(tmid1, tmid,xyz);
485	francois	283
486			param.ajouter(xyz[0]);
487			param.ajouter(xyz[1]);
488			param.ajouter(xyz[2]);
489			param.ajouter(0.5);
490
491	francois	290	evaluer(tmid,xyz);
492	francois	283
493			param.ajouter(xyz[0]);
494			param.ajouter(xyz[1]);
495			param.ajouter(xyz[2]);
496			param.ajouter(1);
497
498	francois	290
499			get_intersection(tmid,tmid2,xyz);
500	francois	283
501			param.ajouter(xyz[0]);
502			param.ajouter(xyz[1]);
503			param.ajouter(xyz[2]);
504			param.ajouter(0.5);
505
506
507	francois	290	get_intersection(tmid2,tcs2,xyz);
508	francois	283
509			param.ajouter(xyz[0]);
510			param.ajouter(xyz[1]);
511			param.ajouter(xyz[2]);
512			param.ajouter(0.5);
513
514			evaluer(tcs2,xyz);
515
516			param.ajouter(xyz[0]);
517			param.ajouter(xyz[1]);
518			param.ajouter(xyz[2]);
519			param.ajouter(1);
520
521			indx_premier_ptctr=17;
522			}
523
524			}
525
526
527			if ((type==MGCo_ELLIPSE\|\|type==MGCo_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*PI,0.000001))
528			courbe->get_param_NURBS(indx_premier_ptctr,param);
529
530			if (type==MGCo_BSPLINE)
531			courbe->get_param_NURBS(indx_premier_ptctr,param);
532
533			if (type==MGCo_LINE)
534			{
535
536			param.ajouter(1);
537			param.ajouter(2);
538			param.ajouter(0);
539			param.ajouter(2);
540			param.ajouter(0);
541
542
543			param.ajouter(0);
544			param.ajouter(0);
545			param.ajouter(1);
546			param.ajouter(1);
547
548			double xyz1[3];
549			double xyz2[3];
550
551			double tcs1=get_tmin();
552			double tcs2=get_tmax();
553
554			evaluer(tcs1,xyz1);
555
556			param.ajouter(xyz1[0]);
557			param.ajouter(xyz1[1]);
558			param.ajouter(xyz1[2]);
559			param.ajouter(1);
560
561			evaluer(tcs2,xyz2);
562
563			param.ajouter(xyz2[0]);
564			param.ajouter(xyz2[1]);
565			param.ajouter(xyz2[2]);
566			param.ajouter(1);
567
568			indx_premier_ptctr=9;
569			}
570
571			}
572
573			struct integrale {
574			double ti;
575			double li;
576			};
577
578			double MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
579			{
580			double pt[3];
581			double t,ti,tii;
582
583			//double longueur_voulue,longueur=0;
584			double longueur=0;
585			std::vector<integrale> tab_integrale;
586			//longueur_voulue= alpha *get_longueur(t1,t2);
587
588			integrale pas;
589			ti=t1;
590			do
591			{
592			tii=ti+(t2-t1)*dt;
593
594			t=0.7886751345ti+0.2113248654tii;
595			deriver(t,pt);
596
597			double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
598			longueur=longueur+0.5(tii-ti)sqrt(facteur);
599			t=0.7886751345tii+0.2113248654ti;
600			deriver(t,pt);
601			facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
602			longueur=longueur+0.5(tii-ti)sqrt(facteur);
603			integrale pas;
604			pas.ti=ti;
605			pas.li=longueur;
606			tab_integrale.insert(tab_integrale.end(),pas);
607			ti=tii;
608			}
609			while (fabs(tab_integrale.back().li)<fabs(longueur_voulue));
610
611			//double t_trouve=0.7886751345tab_integrale.back().ti+0.2113248654tab_integrale.back().ti+(t2-t1)*dt;
612			double tn_1=tab_integrale.back().ti;
613			double ln_1=tab_integrale[tab_integrale.size()-2].li;
614			double tn=ti;
615			double ln=tab_integrale.back().li;
616
617
618			double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
619			return t_trouve ;
620			}
621
622	francois	604	void MG_ARETE::recupere_resultat(class FEM_SOLUTION* sol,int numchamps,char* fichier)
623			{
624			FEM_MAILLAGE *mai=sol->get_maillage();
625			MG_SOMMET* som1=cosommet1->get_sommet();
626			MG_SOMMET* som2=cosommet2->get_sommet();
627			sol->active_solution(numchamps);
628
629			FEM_NOEUD* noeudcourant=NULL;
630			for (int i=0;i<som1->get_lien_fem_maillage()->get_nb();i++)
631			{
632			FEM_NOEUD* no=(FEM_NOEUD*)som1->get_lien_fem_maillage()->get(i);
633			if (mai->get_fem_noeudid(no->get_id())==no) noeudcourant=no;
634			}
635			FEM_NOEUD* noeudfin=NULL;
636			for (int i=0;i<som2->get_lien_fem_maillage()->get_nb();i++)
637			{
638			FEM_NOEUD* no=(FEM_NOEUD*)som2->get_lien_fem_maillage()->get(i);
639			if (mai->get_fem_noeudid(no->get_id())==no) noeudfin=no;
640			}
641			int n1=0;
642			int n2;
643			if (mai->get_degre()==1) n2=1;
644			if (mai->get_degre()==2) n2=2;
645			double *xyz=noeudcourant->get_coord();
646			double t;
647			inverser(t,xyz);
648			std::vector<double> resultat;
649			double tprec=t;
650			resultat.push_back(t);
651			resultat.push_back(0.);
652			resultat.push_back(noeudcourant->get_solution());
653			while (noeudcourant!=noeudfin)
654			{
655			FEM_ELEMENT1 *segcourant=NULL;
656			for (int i=0;i<get_lien_fem_maillage()->get_nb();i++)
657			{
658			FEM_ELEMENT1* seg=(FEM_ELEMENT1*)get_lien_fem_maillage()->get(i);
659			if (mai->get_fem_element1id(seg->get_id())==seg)
660			{
661			if ((seg->get_fem_noeud(n1)==noeudcourant)\|\|(seg->get_fem_noeud(n2)==noeudcourant))
662			if (seg!=segcourant)
663			segcourant=seg;
664			}
665			}
666			FEM_NOEUD* autrenoeud=NULL;
667			if (segcourant->get_fem_noeud(n1)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n2);
668			if (segcourant->get_fem_noeud(n2)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n1);
669			FEM_NOEUD *noprec=noeudcourant;
670			for (int i=1;i<segcourant->get_nb_fem_noeud();i++)
671			{
672			FEM_NOEUD* no=segcourant->get_fem_noeud(i);
673			double *xyz=no->get_coord();
674			double t2;
675			inverser(t2,xyz);
676	francois	672	if (get_courbe()->est_periodique())
677			if (t2<tprec) t2=t2+get_courbe()->get_periode();
678	francois	604	double s=resultat[resultat.size()-2];
679			int pas=32;
680			for (int j=0;j<pas;j++)
681			{
682			double ti=tprec+1.0j(t2-tprec)/pas;
683			double tii=tprec+1.0(j+1)(t2-tprec)/pas;
684			double tgauss1=0.7886751345tii+0.2113248654ti;
685			double dxyz[3];
686			deriver(tgauss1,dxyz);
687			s=s+0.5(tii-ti)sqrt(dxyz[0]dxyz[0]+dxyz[1]dxyz[1]+dxyz[2]*dxyz[2]);
688			double tgauss2=0.7886751345ti+0.2113248654tii;
689			deriver(tgauss2,dxyz);
690			s=s+0.5(tii-ti)sqrt(dxyz[0]dxyz[0]+dxyz[1]dxyz[1]+dxyz[2]*dxyz[2]);
691			}
692			resultat.push_back(t2);
693			resultat.push_back(s);
694			resultat.push_back(no->get_solution());
695			noprec=no;
696			tprec=t2;
697			}
698			noeudcourant=autrenoeud;
699			}
700			FILE *in=fopen(fichier,"wt");
701			char mess[500];
702			sprintf(mess,"%s_%s",(char)sol->get_nom().c_str(),(char)sol->get_legende(numchamps).c_str());
703			fprintf(in,"t;s;%s;\n",mess);
704			for (int i=0;i<resultat.size();i=i+3)
705			fprintf(in,"%lf;%lf;%lf;\n",resultat[i],resultat[i+1],resultat[i+2]);
706			fclose(in);
707			}