* fichier :  lignecourant.dgibi

GRAPH = faux ;

***********************************************************************
* CALCUL DES LIGNES DE COURANT : ON TESTE LES OPERATEURS LAPN ET FIMP *
* ON VERIFIE QUE L'ERREUR ENTRE PSI_NUM ET PSI_EXACT EVOLUE DE FACON  *
* QUADRATIQUE                                                         *
* H. PAILLERE/TTMF/JUILLET 1997                                       *
* ON CONSIDERE LE DOMAINE [0,1]x[0,1] ET LA FONCTION DE COURANT       *
* PSI_EXACT = x(x-1)y^2(y-1)                                          *
* ASSOCIEE AU CHAMP DE VITESSE (u,v) DONNE PAR :                      *
*   u_ex = dPSI/dy et v_ex = -dPSI/dx                                 *
* ON RESOUD L'EQUATION DE POISSON :                                   *
*              d2PSI/dx2 + d2PSI/dy2 = du_ex/dy - dv_ex/dx            *       
* AVEC COMME C.L PSI=0 SUR LES BORDS DU DOMAINE                       *
***********************************************************************

OPTI DIME 2 ;
OPTI ELEM QUA4 ;
OPTI ISOV SULI ;

REPETER BLOC2 6 ;

n = (&BLOC2) ;
NX = 5*n  ;
NY = 5*n ;

A1 = 0.0 0.0 ;
A2 = 1.0 0.0 ;
A3 = 1.0 1.0 ;
A4 = 0.0 1.0 ;

DOMTOT = 'DALLER' (A1 'DROI' NX A2) (A2 'DROI' NY A3)
                  (A3 'DROI' NX A4) (A4 'DROI' NX A1) 'PLAN' ;

MDOMTOT = CHAN DOMTOT QUAF ;
$DOMTOT = 'MODE' MDOMTOT 'NAVIER_STOKES' LINE;
DOMA $DOMTOT 'IMPR' ;
CNT = CONT DOMTOT ;

XX YY = 'COOR' DOMTOT ;

******************************************************************
* EVALUATION DE LA FONCTION DE COURANT EXACTE + CHAMP DE VITESSE *
******************************************************************

PSIEX = KOPS XX '*' ((XX-1.0)) ;
PSIEX = KOPS PSIEX '*' (KOPS YY '*' YY) ;
PSIEX = KOPS PSIEX '*' ((YY-1.0)) ;
PSIEX = KCHT $DOMTOT SCAL SOMMET PSIEX ;

SI ( GRAPH ) ;
TRACE PSIEX DOMTOT 14 'TITR'
        'FONCTION DE COURANT EXACTE' ;
FINSI ;

U = KOPS XX '*' ((XX-1.0)) ;
U = KOPS U '*' YY ;
G = KOPS (YY*3.0) '-' 2.0 ;
U = KOPS U '*' G ;
G = KOPS (XX*2.0) '-' 1.0 ;
V = KOPS G '*' ((YY-1.0)) ;
V = KOPS V '*' YY ;
V = KOPS V '*' YY ;
V = KOPS V '*' (-1.0) ;

U = NOMC 'UX' U 'NATU' 'DISCRET' ;
V = NOMC 'UY' V 'NATU' 'DISCRET' ;
UU = KCHT $DOMTOT VECT SOMMET (U ET V) ;
UNCH = VECT UU 1.0 UX UY VERT ;

SI ( GRAPH ) ;
TRACE UNCH DOMTOT 'TITR' 'CHAMP DE VITESSE' ;
FINSI ;

******************************************************************
* ON CALCULE LA FONCTION DE COURANT NUMERIQUEMENT A PARTIR DE LA *
* DONNEE DU CHAMP DE VITESSE                                     *
******************************************************************

SW = KOPS UU 'ROT' $DOMTOT ;

RK = EQEX $DOMTOT 'OPTI' 'EF' 'IMPL'
ZONE $DOMTOT OPER LAPN 1.0   INCO 'PSI'
ZONE $DOMTOT OPER FIMP sw   INCO 'PSI'
'CLIM' 'PSI' 'TIMP' CNT 0. ;

RK.'INCO' = TABLE INCO ;
RK.'INCO'.'PSI' = KCHT $DOMTOT SCAL SOMMET 0. ;
EXEC RK ;

PSI = RK.'INCO'.'PSI' ;

SI ( GRAPH ) ;
TRACE PSI DOMTOT CNT 14 'TITR' 
        'FONCTION DE COURANT CALCULEE' ;
FINSI ;

************************************************************
* CALCUL DE L'ERREUR ENTRE FONCTION DE COURANT CALCULEE ET *
* FONCTION DE COURANT ANALYTIQUE                           *
************************************************************

ER0 = KOPS PSI '-' PSIEX ;
ER2 = KOPS ER0 '*' ER0 ;
ERR = 0. ;
NPTD=NBNO DOMTOT ;
REPETER BLOC1 NPTD ;
P1 = (DOMA $DOMTOT SOMMET) POIN &BLOC1 ;
ERR = ERR + ('EXTR' ER2 'SCAL' P1) ;
FIN BLOC1 ;
ERR = ERR/NPTD ;
ERR = ERR ** 0.5 ;

MESSAGE '======================================================' ;
MESSAGE 'ERREUR EN NORME L2   = ' ERR ;
MESSAGE '======================================================' ;

SI ( (n 'NEG' 1) ) ;
        LISTE n ;
        PERR = PROG (LOG ERR) ;
        ERRPSI = ERRPSI ET PERR ;
        PXXX = PROG (LOG (1./NX)) ;
        ERRXXX = ERRXXX ET PXXX ;
SINON ;
        ERRPSI = PROG (LOG ERR) ;
        ERRXXX = PROG (LOG (1./NX)) ;
FINSI ;

FIN BLOC2 ;

SI ( GRAPH ) ;
EVOL1 = EVOL 'MANU' 'LOG(DX)' ERRXXX 'LOG(ERR(PSI))' ERRPSI ;
TAB1 = TABLE ;
TAB1.1 = 'MARQ CROI ' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE'. 1 = MOT 'MAILLAGE QUA4' ;
DESS EVOL1 LEGE TAB1 'XBOR' -4. 0. 'YBOR' -11.0 -7.0 'GRIL' ;
FINSI ;

SI ( (MINI ERRPSI) > -10.9 ) ;
        ERREUR 5 ;
FINSI ;

FIN ;