* fichier :  recirc.dgibi
*
*     Exemple test : recirculation dans une cavite semi circulaire
*                    Eclt laminaire incompressible
*                    On impose une vitesse tg sur la surface libre
*                    V normale nulle sur le fond
*
*     Post  : debit sur le contour , fonction de courant , pression
*     Teste vitesse normale imposee en 2D plan algorithme Taylor Hood
*
*     Voir aussi RECIRCP.DGIBI pour algorithme de projection

opti dime 2 elem QUA8 ;
DISCR = 'QUAF';
KPRES= 'CENTREP1' ;

GRAPH=FAUX ;
*GRAPH=VRAI ;
nbit = 5 ;

p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;

n1=  5 ; n2= 5 ;
ab= p1 d n1 p2 ;
slib=  (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));

bc= p2 c p0 n2 p03 c  p0 n2 p1;
*bc= p2 d n2 p3 d n2 p4 d n2 p1 ;

cnt=ab et bc ;

mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;

$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;

* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt KPRES ;
bcp=elem  prep1 POI1 (lect 1) ;
*

*NU= 1. / 61. ;
NU= 1.e-2 ;
U0=1. ;

rv = eqex 'ITMA' 0 'NITER' nbit  'OMEGA' 0.8
  OPTI EF 'CENTREE' 'IMPL' KPRES
  'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU  'INCO' 'UN'
  'ZONE' $MT 'OPER' 'KBBT' 1.  'INCO' 'UN' 'PRES'
  'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
  'CLIM'
         'UN' UIMP slib 1. 'UN' VIMP ab 0.
         'PRES' TIMP bcp 0. ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES  0. ;
rv.'INCO'.'VN'=kcht $bc scal sommet 0. ;

exec rv ;

un = rv.inco.'UN' ;
ung= vect un 0.11  ux uy jaune ;

Si GRAPH ;trace ung Mmt ; Finsi ;

cnt1=chan cnt seg2 ;                                                    
Mcnt1= chan cnt1 quaf;                                                  
$cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
q=dbit (rv.inco.'UN') $cnt ;                                            
Mess ' Debit sur le contour ' q ;
si ( q  > 1.5e-15 ) ; erreur 5 ; finsi ;

*
* Calcul de la fonction de courant
*

sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
     'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
     'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
     'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;

opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .26820) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-5 ) ; erreur 5 ; finsi ;

*
* Calcul de la pression
*
  pn= elno (rv.inco.'PRES') $mt KPRES;
Si GRAPH ;  trace pn Mmt ; Finsi ;


*
*  Deuxieme cas : Meme probleme geometrie rectangulaire
*                 MACRO CENTRE
*
*
*
*
*
*

DISCR = 'MACRO';
KPRES= 'CENTRE' ;

p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;

n1=  5 ; n2= 5 ;
ab= p1 d n1 p2 ;
slib=  (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));

bc= p2 d n2 p3 d n2 p4 d n2 p1 ;

cnt=ab et bc ;

mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;

$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;

* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt KPRES ;
bcp=elem  prep1 POI1 (lect 1) ;
*

NU= 1.e-2 ;
U0=1. ;

rv = eqex 'ITMA' 0 'NITER' nbit  'OMEGA' 0.8
  OPTI EF 'CENTREE' 'IMPL' KPRES
  'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  'ZONE' $MT 'OPER' 'KBBT' 1.  'INCO' 'UN' 'PRES'
  'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
  'CLIM'
         'UN' UIMP slib 1. 'UN' VIMP ab 0.
         'PRES' TIMP bcp 0. ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
rv.'INCO'.'VN'=0. ;

exec rv ;


un = rv.inco.'UN' ;
ung= vect un 0.11  ux uy jaune ;

Si GRAPH ;trace ung Mmt ; Finsi ;

cnt1=chan cnt seg2 ;                                                    
Mcnt1= chan cnt1 quaf;                                                  
$cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
q=dbit (rv.inco.'UN') $cnt ;                                            
Mess ' Debit sur le contour ' q ;
si ( q  > 1.5e-15 ) ; erreur 5 ; finsi ;

sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
     'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
     'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
     'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;

opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
mess ' MINI PSI = ' (mini psi) ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .22) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;


*
* Calcul de la pression
*
  pn= elno (rv.inco.'PRES') $mt KPRES;
Si GRAPH ;  trace pn Mmt ; Finsi ;
FIN ;