* fichier :  cc3d1.dgibi
************************************************************************
* Section : Fluides Transitoire
************************************************************************
**
**  --- 08 Novembre 1999 ---
**
**   TEST CAVITE CUBIQUE
**
**   teste KBBT NS en 3D + le Bi CG cc2d1.dgibi
**   teste les elements a pression continue sur des tetrahedres
**   et pyramides P1-P1 et Q1-Q1 algorithme de projection
**   stabilises en prenant DT > DT limite
**   Teste aussi la creation de volumes automatique
**   Si on ne veut que des tetrahedres  activer option elem tri3

GRAPH = 'N'                   ;
err1=5.e-2;
DISCR='LINE' ;
kpress='MSOMMET';
COMPLET = FAUX ;

SI(COMPLET) ;
 ds1=0.05 ; ds2=0.08 ;
ITMAX=20 ;
SINON ;
 ds1=0.1  ; ds2=0.1 ;
ITMAX=6;
FINSI ;

 option dime 2 elem qua4 ;
*option dime 2 elem tri3 ;

opti isov suli ;
p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;

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

prof=0.1;
mth=mt plus (0 0 prof) ;
ab= ab12 plus (0 0 prof) ;
oeil = 10 10 100 ;
nph=2 ;
cav=mt volu  nph mth ;
*trace oeil cav;

f1=face 1 cav ;
f2=face 2 cav ;
f3=face 3 cav ;

 psup=cd trans nph (0 0 prof) ;
 psup=(CHAN TRI3 psup) ;
 ff3= (CHAN TRI3 f3) ;
 ELIM (PSUP ET FF3) 1.e-4 ;
 env3= f1 et f2 et ff3 ;

 cav = volu env3 ;

cav= chan cav quaf ;
ab=chan ab quaf ;

elim (cav et psup et ab et ff3 ) 1.e-5 ;

 $mt=mode cav 'NAVIER_STOKES' DISCR ;
 $AB=mode AB 'NAVIER_STOKES'  DISCR ;
  DOMA $mt IMPR ;

DT=2. ;
MU=1. ;
RO= 400. ;

prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 1 ) ;

  rv= eqex 'OMEGA' 1.   'NITER' 1  ITMA ITMAX
  'OPTI' 'EF' 'IMPL' KPRESS 'SUPG'
  ZONE  $mt       OPER    NS 1. 'UN' (MU/RO)   INCO 'UN'
  'OPTI' 'EF' 'BDF2'
  ZONE  $mt    OPER  DFDT 1.  'UNM' 'UNMM' DT  'UN' MU INCO 'UN'
      ;

  rv=eqex rv
  CLIM
  UN UIMP ff3 0. UN VIMP ff3 0. UN WIMP ff3  0.
  UN WIMP (f1 et f2) 0.
  UN UIMP PSUP 1.
  ;


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.'UNM' =  kcht $mt  vect sommet (1.e-5 1.e-5 1.e-5) ;
rv.inco.'UNMM' =  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 ;


RVP= EQEX
    'OPTI' 'EF' KPRESS
    ZONE  $mt    OPER  KBBT  (-1.)   INCO 'UN' 'PRES'
    CLIM PRES TIMP bcp 0.
    ;

    rvp.'METHINV'.TYPINV=2 ;
    rvp.'METHINV'.IMPINV=1 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;

    RV.'PROJ'= RVP ;

 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 ;
 trace c1 cav ;
 trace (rv.inco.'PRESSION') mt ;
 finsi ;

lrr=prog
 3.76000E-38  2.84554E-02  4.47889E-02  6.01943E-02  7.05196E-02
 6.69563E-02  4.10891E-02  4.04425E-03  5.94186E-02  0.26219
   1.0000 ;

 evx = ABS evx ;
 list evx ;
 ER=SOMM( abs (evx - lrr) ) *0.1;
 mess ' Ecart sur MSOMMET TET4 PYR5 ' er ;
 Si ( er > 5.e-3) ; erreur 5 ; finsi ;
 FIN ;