/branches/FREEFEM/freefem/lib/adapt/cad.cxx

/*
 * ADAPT2D : a software for automatic mesh adaptation in 2D
 *
 * AUTHOR : Manuel J. Castro Diaz(e-mail:castro@gamba.cie.uma.es)
 * ADAPTED FOR FREEFEM : Prud'homme Christophe (e-mail:prudhomm@ann.jussieu.fr)
 *
 * this code is public domain
 * 
 * You may copy freely these files and use it for    
 * teaching or research. These or part of these may   
 * not be sold or used for a commercial purpose without
 * our consent
 * 
 * Any problems should be reported to the AUTHOR
 * at the following address : castro@gamba.cie.uma.es
 */



#include <cad.hxx>
   
/*    
________________________________________________________________________
                       
                         SUBRRUTINA TESTEO      

    Testea si se puden ir uniendo trozos de frontera para formar una
    curva cerrada. Si los puede unir los unira en un sentido y otro
    segun pueda. El resultado final se una lista (CAD)
    cuyos elementos seran a su vez una lista de aristas, estas ordenadas
    y de tal manera que forman curvas. El numero de elementos de la
    lista (CAD) nos da el numero de componentes conexas
    de la frontera.
    
    Entrada: Una lista desordenada de aristas...
    Salida:  Las aristas ordenadas formando curvas. 
*/
Boolean CAD::testeo() {
  FraristaCAD_link* aux,*aux1,*sig;
  Arista_CAD_dlist *l0_aux;
  Arista_CAD_dlink *lk_aux,*prev0;
  R2 c0,c1;
  int ref;
  Boolean modif=FALSE;
  aux=begin;
  while(aux) {
    aux1=aux->suc;
    while (aux1) {
      if ((*aux)<(*aux1)) {   
        sig=aux1;
        aux1=aux1->suc;
        enlaza(aux,sig,1);
        delete sig;
        modif=TRUE;
      }
      else {
        if ((*aux1)<(*aux)) {
          sig=aux1;
          aux1=aux1->suc;
          enlaza (aux,sig,-1);
          delete sig;
          modif=TRUE;
        }
        else {
          c0=aux->Aretelist->last->a->c[1];
          c1=aux1->Aretelist->last->a->c[1];
          if (c0==c1) {
            l0_aux=aux1->Aretelist;
            prev0=NIL;
            lk_aux=l0_aux->last;
            l0_aux->last=l0_aux->begin;
            l0_aux->begin=lk_aux;
            while (lk_aux) {
              c0=lk_aux->a->c[0];
              lk_aux->a->c[0]=lk_aux->a->c[1];
              lk_aux->a->c[1]=c0;
              ref=lk_aux->a->refv[0];
              lk_aux->a->refv[0]=lk_aux->a->refv[1];
              lk_aux->a->refv[1]=ref;
              lk_aux->suc=lk_aux->prev;
              lk_aux->prev=prev0;
              prev0=lk_aux;
              lk_aux=lk_aux->sig();
            }
            sig=aux1;
            aux1=aux1->suc;
            enlaza (aux,sig,1);
            delete sig;
            modif=TRUE;
          }
          else {
            c0=aux->Aretelist->begin->a->c[0];
            c1=aux1->Aretelist->begin->a->c[0];
            if (c0==c1) {
              l0_aux=aux->Aretelist;
              prev0=NIL;
              lk_aux=l0_aux->last;
              l0_aux->last=l0_aux->begin;
              l0_aux->begin=lk_aux;
              while (lk_aux) {
                c0=lk_aux->a->c[0];
                lk_aux->a->c[0]=lk_aux->a->c[1];
                lk_aux->a->c[1]=c0;
                ref=lk_aux->a->refv[0];
                lk_aux->a->refv[0]=lk_aux->a->refv[1];
                lk_aux->a->refv[1]=ref;
                lk_aux->suc=lk_aux->prev;
                lk_aux->prev=prev0;
                prev0=lk_aux;
                lk_aux=lk_aux->sig();
              }
              sig=aux1;
              aux1=aux1->suc;
              enlaza (aux,sig,1);
              delete sig;
              modif=TRUE;
            }
            else aux1=aux1->suc;
          }
        }
      }
    }
    aux=aux->suc;
  }
  return modif;
}

/*  --------------------------------------------------------------------

                       SUBRRUTINA D_FRONTERA

   Proposito:  Dado un punto y la Frontera de un dominio, busco la
   ---------   arista y por tanto, el elemento mas proximo al punto
               dado. Calcularemos tambien el proyectado sobre la arista.

   Entrada:    pt= Punto dado.
   -------     - Frontera apuntada por this.
               dista=distancia inicial, si no se conoce
                     poner infinito.
   Salida:     proyec= proyectado sobre la frontera.
   ------      aret=  arista mas cercana.
               trian= elemento que contine la arista.


Scalar CAD::d_frontera(R2 pt, R2& proyec,\
         Arista_T0& arete, Triangulo_T0& triangle,const Scalar& dista) {

         FraristaT0_link* aux1;
         AristaT0_list* comp_conex;
         AristaT0_link* aux;
         Arista_T0* aret;
         Scalar d0,dist=dista;
         R2 pr_pt(0,0);
         aux1=begin;

         while (aux1) {
           comp_conex=aux1->Aretelist;
           aux=comp_conex->principio();
           while(aux) {
             aret=aux->a;
             d0=d_arete(*aret,pt,pr_pt);
             if (d0<dist) {
                 dist=d0;
                 arete=*aret;
                 triangle=*(aret->tr);
                 proyec=pr_pt;
             }
             aux=aux->sig();
           }
           aux1=aux1->sig();
         }
         return dist;
    }

    */
void CAD::write (ostream& os) {
  int num=0,i=0;
  FraristaCAD_link* aux;
  num=this->num_elem();
  os <<endl<<"Numero de curvas de la frontera:"<<num<<endl;
  aux=begin;
  while (aux) {
    i++;
    os <<aux<<" Curva:"<<i<<" constituida por las Aristas"<<endl;
    os <<(*(aux->Aretelist))<<endl;
    aux=aux->suc;
  }
}

ostream& operator<<(ostream& os, CAD& fl)
{ fl.write (os); return os;}

CAD* crea_frontera(Mallado_T0* malla,int angulo) {
  CAD* fr_aux,*fr;
  FraristaCAD_link* comp_conex,*cc;
  Arista_CAD_dlist* arete_list,*al;
  Arista_CAD_dlink* arete,*arete_sig,*aa;
  Triangulo_T0* t0;
  Arista_T0* a0;
  Arista_CAD* a_cad;
  int nba=malla->nbaa();
  int i,ref0,ref1,pos0,nbcc,ang0;
  Boolean modif,fin;
  fr_aux=new CAD;
  if (fr_aux==NIL) ERROR();
  for (i=0; i<nba; i++) {
    a0=malla->saca_arete(i);
    if (a0->front==TRUE) {
      a_cad=new Arista_CAD;
      if (a_cad==NIL) ERROR();
      a_cad->set(a0->s[0]->c,a0->s[1]->c,a0->s[0]->ref,a0->s[1]->ref,a0->ref,1);
      
      arete=new Arista_CAD_dlink;
      if (arete==NIL) ERROR();
      arete->set(a_cad,NIL,NIL);
      
      arete_list=new Arista_CAD_dlist;
      if (arete_list==NIL) ERROR();
      arete_list->append(arete);
      
      comp_conex=new FraristaCAD_link;
      if (comp_conex==NIL) ERROR();
      comp_conex->set(arete_list,NIL,NIL);
      
      fr_aux->append(comp_conex);
      modif=fr_aux->testeo();
    }
  }
  for (i=0; i<nba; i++) {
    a0=malla->saca_arete(i);
    if (a0->front==FALSE) {
      t0=a0->tr;
      ref0=t0->ref;
      pos0=t0->arista(a0);
      if (pos0<0) {
        cerr<<"Edge is not in triangle."<<endl;
        cerr<<a0<<" "<<t0<<endl;
        cerr<<"Error in crea_frontera. (Frontera_CAD.C)"<<endl;
        exit(-1);
      }
      ref1=t0->t[pos0]->ref;
      if (ref0!=ref1) {
        a_cad=new Arista_CAD;
        if (a_cad==NIL) ERROR();
        a_cad->set(a0->s[0]->c,a0->s[1]->c,a0->s[0]->ref,a0->s[1]->ref,a0->ref,2);
        arete=new Arista_CAD_dlink;
        if (arete==NIL) ERROR();
        arete->set(a_cad,NIL,NIL);
        
        arete_list=new Arista_CAD_dlist;
        if (arete_list==NIL) ERROR();
        arete_list->append(arete);
        
        comp_conex=new FraristaCAD_link;
        if (comp_conex==NIL) ERROR();
        comp_conex->set(arete_list,NIL,NIL);
        
        fr_aux->append(comp_conex);
        modif=fr_aux->testeo();
      }
    }
  }
  modif=TRUE;
  while (modif) {modif=fr_aux->testeo();}

  fr=new CAD;
  if (fr==NIL) ERROR();
  nbcc=fr_aux->num_elem();
  comp_conex=fr_aux->principio();//cojo la primera componente
                                // conexa de la frontera.    
  for (i=0; i<nbcc; i++) {
    arete_list=comp_conex->Aretelist; //cojo la lista de aristas de la
                                // i-esima comp. conexa.
    arete=arete_list->principio();  //cojo el principio de la lista de
                                // aristas.
    a_cad=new Arista_CAD;
    if (a_cad==NIL) ERROR();
    *a_cad=*arete->a;
    aa=new Arista_CAD_dlink;
    if (aa==NIL) ERROR();
    aa->set(a_cad,NIL,NIL);
    
    al=new Arista_CAD_dlist;
    if (al==NIL) ERROR();
    al->append(aa);
    
    cc=new FraristaCAD_link;
    if (cc==NIL) ERROR();
    cc->set(al,NIL,NIL);
    
    fr->append(cc);
    fin=FALSE;
    while (arete) {
      arete_sig=arete->sig();
      if (arete_sig!=NIL) {
        ang0=int(angle(arete->a,arete_sig->a));
        if (ang0<=angulo) {
          a_cad=new Arista_CAD;
          if (a_cad==NIL) ERROR();
          *a_cad=*arete_sig->a;
          aa=new Arista_CAD_dlink;
          if (aa==NIL) ERROR();
          aa->set(a_cad,NIL,NIL);
          al->append(aa);
        } 
        else {
          a_cad=new Arista_CAD;
          if (a_cad==NIL) ERROR();
          *a_cad=*arete_sig->a;
          aa=new Arista_CAD_dlink;
          if ( aa==NIL) ERROR();
          aa->set(a_cad,NIL,NIL);
          al=new Arista_CAD_dlist;
          if (al==NIL) ERROR();
          al->append(aa);        
          cc=new FraristaCAD_link;
          if (cc==NIL) ERROR();
          cc->set(al,NIL,NIL);
          fr->append(cc);
        }
      }
      arete=arete->sig();
    }
    comp_conex=comp_conex->sig();
  }
  delete fr_aux;
  return fr;
}
CAD* build(Triangulation* t_cad, int angulo, Scalar Scale) {
  Mallado_T0* m_cad;
  CAD* cad=NIL;
  m_cad=new Mallado_T0;
  if (m_cad==NIL) ERROR();

  m_cad->build(int(t_cad->np),int(t_cad->nt),t_cad->rp,t_cad->ng,t_cad->tr,t_cad->ngt,Scale);
  if (m_cad->m) {delete[] m_cad->m; m_cad->m=NIL;}
  cad=crea_frontera(m_cad,angulo);
  cad->set_angulo(angulo);
  delete m_cad;
  cad->P_control();
  return cad;
}

  
/*
-----------------------------------------------------------------------------

                             SUBRRUTINA TANGENTES

    A cada arista que esta en la frontera asocia en cada vertice dos
 tangentes. Estas se utilizaran a la hora de definir la geometria, para
 hacer pasar una curva C1 por el borde.

   Parametros de entrada: angulo de conservacion de las esquinas.
   Salida las tangentes en cada arista. Ademas los vertices que sean
   el vertice de un angulo seran marcados con una referencia negativa.
*/

void CAD::Tangentes (int ang) {
  FraristaCAD_link* comp_conex;
  Arista_CAD_dlist* arete_list;
  Arista_CAD_dlink* arete,*arete_sig,*arete_ant;
  R2 c0,c1;
  R2 tang0,tang1;
  R2 rc0,rc1,rc2,rc3;
  int i,nbcc;
  int ang0,ang1;
  Boolean fin=FALSE,cerrada;
  if (begin!=NIL) {
    nbcc=this->num_elem();
    comp_conex=this->principio();//cojo la primera componente
                                // conexa de la frontera.    
    for (i=0; i<nbcc; i++) {
      arete_list=comp_conex->Aretelist; //cojo la lista de aristas de la
                                // i-esima comp. conexa.
      arete=arete_list->principio();  //cojo el principio de la lista de
                                // aristas.
      arete_ant=arete_list->fin();
      c0=arete->a->c[0];
      c1=arete_ant->a->c[1];
      if (c0==c1) { // la curva es cerrada
        ang0=int(angle(arete_ant->a,arete->a));
        cerrada=TRUE;
      }
      else{
        cerrada=FALSE;
        ang0=360;
      }
      fin=FALSE;
      while (fin==FALSE) {
        arete_sig=arete->sig();
        if (arete_sig==NIL) {
          arete_sig=arete_list->principio();
          fin=TRUE;
        }    
        rc1=arete->a->c[0];
        rc2=arete->a->c[1];
        rc0=arete_ant->a->c[0];
        rc3=arete_sig->a->c[1];
        if (fin==TRUE && cerrada==FALSE) 
          ang1=360;
        else 
          ang1=int(angle(arete->a,arete_sig->a));
        if (ang0<=ang) 
          tang0=rc2-rc0;
        else { 
          tang0=rc2-rc1;
        }
        if (ang1<=ang)
          tang1=rc3-rc1;
        else {
          tang1=rc2-rc1;
        }
//        tang=new Tangente;
//        if (tang==NIL) ERROR();
//        tang->set(tang0,tang1);
        arete->a->tg.set(tang0,tang1);
        arete_ant=arete;
        arete=arete_sig;
        ang0=ang1;
      }
      comp_conex=comp_conex->sig();
    }
  }
}
void CAD::P_control () {
  FraristaCAD_link* comp_conex;
  Arista_CAD_dlist* arete_list;
  Arista_CAD_dlink* arete,*aux;
  R2 p_ini,p_fin,p,psig;
  Scalar m[4][4],cox,coy;
  int i,j,k,l,nbcc,nba,cont;
  Boolean cerrada,modif;

  m[0][0]=-1;m[0][1]=2;m[0][2]=-1;m[0][3]=0;
  m[1][0]=3;m[1][1]=-5;m[1][2]=0;m[1][3]=2;
  m[2][0]=-3;m[2][1]=4;m[2][2]=1;m[2][3]=0;
  m[3][0]=1;m[3][1]=-1;m[3][2]=0;m[3][3]=0;

  if (begin!=NIL) {
    nbcc=this->num_elem();
    comp_conex=this->principio();//cojo la primera componente
                                // conexa de la frontera.    
    for (i=0; i<nbcc; i++) {
      arete_list=comp_conex->Aretelist; //cojo la lista de aristas de la
                                        // i-esima comp. conexa.      
      arete=arete_list->principio();    //cojo el principio de la lista de
      cerrada=FALSE;                     // aristas.
      aux=arete_list->fin();
      if (arete->a->c[0]==aux->a->c[1]) cerrada=TRUE;
      if (cerrada==TRUE) {
        p_ini=aux->a->c[0];
        p_fin=arete->a->c[1];
      }
      else {
        p_ini=arete->a->c[0]*2-arete->a->c[1];
        p_fin=aux->a->c[1]*2-aux->a->c[0];
      }
      nba=arete_list->num_elem();
      for (j=0;j<nba;j++) {
//        cout<<" j:"<<j<<"  --------------------------------------------------"<<endl;
        for (k=0; k<4; k++) {
          cox=coy=0.0;
          aux=arete->ant();
          if (aux==NIL){ 
            p=p_ini;
            psig=arete->a->c[0];
          }
          else {
            p=aux->a->c[0];
            psig=aux->a->c[1];
          }
          cont=0;
          for (l=0; l<4; l++) {
//            if (k==0) cout<<p<<endl;
            cox+=m[l][k]*p.x;
            coy+=m[l][k]*p.y;
            if (arete->ant()==NIL && aux==NIL) {
              aux=arete;
            }
            else {
              if (aux) aux=aux->sig();
            }
            if (aux) {
              p=aux->a->c[0];
              psig=aux->a->c[1];
            }
            else {
              if (cont==0){ 
                p=psig;
                cont=1;
              }
              else 
                p=p_fin;
            }
          }
          if (arete) arete->a->cf[k].set(cox,coy);
        }
        arete=arete->sig();
      }
      comp_conex=comp_conex->sig();
    }
    modif=TRUE;
    while (modif) {modif=this->testeo();}
  }   
}    

void CAD::Pinta (int np) {
  FraristaCAD_link* comp_conex;
  Arista_CAD_dlist* arete_list;
  Arista_CAD_dlink* arete;
  Arista_CAD* a0;
  R2 pt;
  Scalar temp;
  int j,i,nbcc;
  char path[72]="spline";

  if (begin!=NIL) {
    OPEN(escritura, path);
    escritura.setf(ios::scientific,ios::floatfield);

    nbcc=this->num_elem();
    comp_conex=this->principio();//cojo la primera componente
                                // conexa de la frontera.    
    for (i=0; i<nbcc; i++) {
      arete_list=comp_conex->Aretelist; //cojo la lista de aristas de la
                                // i-esima comp. conexa.
      arete=arete_list->principio();  //cojo el principio de la lista de
                                // aristas.
      
      while (arete) {
        a0=arete->a;
        for (j=0; j<np; j++) {
          temp=j/Scalar(np-1);
          pt=(a0->cf[3]+a0->cf[2]*temp+a0->cf[1]*(temp*temp)+a0->cf[0]*(temp*temp*temp))*0.5;
          escritura<<""<<setw(15)<<pt.x<<setw(15)<<pt.y<<endl;
        }
        arete=arete->sig();
      }
      comp_conex=comp_conex->sig();
    }
    escritura.close();
  }
}