Accueil | Cast3M

* fichier :  cc3d2.dgibi
************************************************************************
************************************************************************
**
**  --- 08 Novembre 1999 ---
**
**   TEST CAVITE CUBIQUE
**
**   teste KBBT NS en 3D + le Bi CG 
**   teste les elements a pression continue sur des hexahedres
**   et prisme P1-P1 et Q1-Q1 algorithme de projection
**   ces elements ne verifient pas LBB mais si DT pas trop petit
**   ca marche quand meme (pression en damier)
 
GRAPH = 'N'                   ;
DISCR='LINE' ;
 
 ds1=0.02 ; ds2=0.1  ;
 ITMAX=20 ;
 err1=1.e-6;
 
 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.) ;
mt0 = ab trans dini ds1 dfin ds2 (0 0.5) trans dini ds2
dfin ds1 (0 0.5) ;
 
mt1= ab trans dini ds1 dfin ds2 (0 0.5);
 
opti elem tri3 ;
 
mt2= (cote 3 mt1) trans dini ds2
dfin ds1 (0 0.5) ;
 
*cnt2=cont mt2 ;
*mt2 = cnt2 surf ;
mt = mt1 et mt2 ;
 
elim (mt0 et mt) 1.e-5 ;
 
si ('EGA' graph 'O' );
trace mt ;
finsi ;
 
 bc=cote 2 mt0 ;
 cd=cote 3 mt0 ;
 da=cote 4 mt0 ;
ct=ab et bc et cd et da ;
elim ct 1.e-3 ;
 cnt= ab et bc et cd et da ;
* mt=ab bc cd da daller ;
* mt=cnt surf  ;
 
mt=orie mt ;
 
opti dime 3 elem cub8;
 
prof=0.1 ;
mth=mt plus (0 0 prof) ;
ab= ab plus (0 0 prof) ;
oeil = 10 10 100 ;
nph=1 ;
cav=mt volu  nph mth ;
 
f1=face 1 cav ;
f2=face 2 cav ;
f3=face 3 cav ;
 
 psup=(ab plus (0 1 0)) trans nph (0 0 ((-1.)*prof)) ;
*psup=(CHAN TRI3 psup) ;
*ff3= (CHAN TRI3 f3) ;
 ff3 = f3 ;
 
 ELIM (PSUP ET FF3) 1.e-4 ;
 
 
cav= chan cav quaf ;
ab=chan ab quaf ;
 srti=da   plus (0.5 0. 0.) ;
 srti = inve srti ;
 elim (srti et cav) 1.e-5 ;
 srti =chan srti quaf ;
 
 
elim (cav et psup et ab et ff3 et srti) 1.e-5 ;
 
 $mt=mode cav 'NAVIER_STOKES' DISCR ;
 $AB=mode AB 'NAVIER_STOKES'  DISCR ;
 $srti=mode srti 'NAVIER_STOKES'  DISCR ;
  DOMA $mt IMPR ;
 srti = doma $srti maillage ;
 
DT=1.;
MU=1. ;
RO= 400. ;
kpress='MSOMMET';
prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 1 ) ;
 
  rv= eqex 'OMEGA' 1.   'NITER' 1  ITMA ITMAX
  'OPTI' 'EF' 'IMPL' KPRESS 'SUPG' 'NODIV'
  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'
  ZONE  $mt    OPER  DFDT 1.  'UNM'        DT  'UN' MU INCO 'UN'
      ;
 
  rv=eqex rv
  CLIM
  UN UIMP ff3 0. UN VIMP ff3 0. UN WIMP ff3  0.
  UN UIMP PSUP 1. UN WIMP (f2 et f1) 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.'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=0 ;
 rv.'METHINV'.NITMAX=100;
 rv.'METHINV'.PRECOND=3 ;
 rv.'METHINV'.RESID  =1.e-8 ;
 
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=0 ;
    rvp.'METHINV'.NITMAX=100;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-10;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;
 
    RV.'PROJ'= RVP ;
 
 exec rv ;
 
 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_UY '     ;
 DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1 ;
 
 c1=vect (rv.inco.'UN') 0.3 ux uy uz jaune ;
 menvf= doma $mt 'ENVELOPPE' ;
 trace c1 menvf ;
 trace (rv.inco.'PRESSION') menvf cnt ;
 finsi ;
 
  list evx ;
lrr=prog
-2.85788E-38 -1.61082E-02 -3.22813E-02 -4.99973E-02 -7.04952E-02
-9.56495E-02  -.12703      -.16245      -.18934      -.18725
 -.12259     -4.13788E-02  3.10240E-02  8.66209E-02   .12808
  .16266       .21239       .31615       .50353       .75483
 1.0000;
 
 
 lr = evx ;
 
 Dgsrti = doma $srti 'XXDIAGSI' ;
 li= 'EXTR' (evol chpo Dgsrti srti) 'ORDO' ;
 
 ER=( SOMM( ((lr - lrr)*(lr - lrr))*li) )** 0.5 ;
 mess ' Ecart sur MSOMMET PRI6 CUB8  ' er ;
 Si ( er > err1 ) ; erreur 5 ; finsi ;
 FIN ;