* fichier :  ccar3d.dgibi
************************************************************************
* Section : Fluides Transitoire
************************************************************************
**
**  --- 12 OCTOBRE 1998 ---
**
**   TEST CAVITE CUBIQUE
**
**   teste KCCT NS en 3D + le Bi CG

GRAPH = 'N'                   ;
err1=5.e-3;

option dime 2 elem tri6 ;
opti isov suli ;
p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;

ds1=0.051  ; ds2=0.12 ;
*ds1=0.01 ; ds2=0.1 ;
 ds1=2.   ; ds2=2.  ;
 ds1=0.25 ; ds2=0.2 ;
 ds1=0.3  ; ds2=0.3 ;
ab=p1 d dini ds1 dfin ds2 p12 d dini ds2 dfin ds1 p2 ;
ab12= p12 d dini ds1 dfin ds2 (0.5 0.5) d dini ds2
dfin ds1 (0.5 1.) ;
mt= ab trans dini ds1 dfin ds2 (0 0.5) trans dini ds2
dfin ds1 (0 0.5) ;
 bc=cote 2 mt ;
 cd=cote 3 mt ;
 da=cote 4 mt ;
ct=ab et bc et cd et da ;
elim ct 1.e-3 ;
 mt=ab bc cd da daller ;

mt=orie mt ;

opti dime 3 elem cu20;
mth=mt plus (0 0 1) ;
ab= ab12 plus (0 0 1) ;
oeil = 10 10 100 ;
nph=2 ;
cav=mt volu  nph mth ;
*trace oeil cav;
macro= '   ' ;
f1=face 1 cav ;
f2=face 2 cav ;
f3=face 3 cav ;
cav= chan cav quaf ;
f1=chan f1 quaf ;
f2=chan f2 quaf ;
f3=chan f3 quaf ;
ab=chan ab quaf ;
psup=cd trans nph (0 0 1) ;
psup= chan psup quaf ;
elim (cav et f1 et f2 et f3 et psup et ab ) 1.e-5 ;
elim (f3 et psup) 1.e-5 ;
paroi=psup diff f3 ;

 $mt=mode cav 'NAVIER_STOKES' QUAF;
 $AB=mode AB 'NAVIER_STOKES'  QUAF;

MU=1. ;
RO= 400. ;
kpress='CENTREP1' ;
prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 1) ;

  rv= eqex 'OMEGA' 0.3  'NITER' 3
  'OPTI' 'EF' 'IMPL' KPRESS 'SUPG' 'DIV2'
  ZONE  $mt       OPER    KBBT (1.)         INCO 'UN' 'PRES'
  ZONE  $mt       OPER    NS 1. 'UN'  (MU/RO)       INCO 'UN'        ;

  rv=eqex rv
  CLIM
  UN VIMP (F3 ET F1) 0. UN WIMP (F1 et F2 et F3)  0.
  UN UIMP PSUP 1.  UN UIMP (PAROI ET F1) 0. PRES TIMP bcp 0.  ;


rv.inco= table inco ;
rv.inco.tn =  kcht $mt  scal sommet 5. ;
rv.inco.un =  kcht $mt  vect sommet (1.e-5 1.e-5 1.e-5) ;

rv.inco.pres =  kcht $mt  scal KPRESS 1.e-5  ;

 rv.'METHINV'.TYPINV=3 ;
 rv.'METHINV'.IMPINV=1 ;
 rv.'METHINV'.NITMAX=1000;
 rv.'METHINV'.PRECOND=3 ;
 rv.'METHINV'.RESID  =1.e-10 ;

 exec rv ;

 srti=doma $AB 'MAILLAGE'        ;
 evolV = EVOL 'CHPO' (rv.'INCO'.'UN') UX (srti  ) ;
 evx=extr evolV 'ORDO' ;
 list evx ;
 evy=extr evolV 'ABSC' ;
 evolV= evol 'MANU' 'Vitesse' evx 'Hauteur' evy ;

 si ('EGA' graph 'O' );
 TAB1=TABLE;
 TAB1.'TITRE'=TABLE ;
 TAB1 . 1 ='MARQ      REGU '     ;

 TAB1.'TITRE' . 1 = mot 'Composante_UX '     ;
 DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1 ;

 c1=vect (rv.inco.'UN') 0.3 ux uy uz jaune ;
 trace c1 mt ;
 pn=elno $mt (rv.inco.'PRES') kpress;
 trace pn mt ;
 finsi ;
 evx=abs evx ;
 lrr='PROG'  3.43000E-06 +3.74144E-02 +6.40402E-02 +6.91359E-02
 +6.26670E-02 +3.46496E-02  5.43639E-02   .17102       .65700  ;
 ER=SOMM( (((abs (evx - lrr))*0.125)**2.) / 1. ) ;
 ER=SOMM( abs (evx - lrr) ) ;
 mess ' Ecart sur CENTREP1 TRI7 QUADR ' er ;
 si ( er > err1 ) ; erreur 5 ; finsi ;
 FIN ;