* fichier :  test_fsur.dgibi
************************************************************************
* Section : Mathematiques Fonctions
* Section : Mecanique Elastique
************************************************************************
*----------------------------------------------------------------------*
* PETIT TEST DE VERIFICATION DE L'OPERATEUR FSUR
*----------------------------------------------------------------------*
'OPTION' 'ECHO' 1 ;
* Mettre IGRAPH a VRAI pour avoir les quelques traces
IGRAPH = FAUX;

*----------------------------------------------------------------------*
* VERIFICATION EN DIMENSION 2
*----------------------------------------------------------------------*
'OPTION' 'DIME' 2 ;

P0 =  0.  0. ;
P1 = 10.  0. ;
P2 =  0. 10. ;

'OPTION' 'ELEM' 'SEG2' ;
L1 = P0 'DROITE' 8 P1 ;

* TEST 1 : COQ2 sur DROITE INCLINEE A 45 degres
* ======== Comparaison de la solution avec PRESS
L2 = P0 'DROITE' 7 (P1 'PLUS' P2) ;

MOD1 = 'MODELISER' L2 'MECANIQUE' 'ELASTIQUE' 'COQ2' ;
MAT1 = 'MATERIAU' MOD1 'YOUN' 200.E+9 'NU' 0.3 'EPAI' 0.1 ;

* Solution avec MODE PRES pression constante = p
MOP = 'MODE' L2 'CHARGEMENT' 'PRESSION' 'COQ2' ;
MAP = 'PRES' MOP 'PRES' 1. ;
PR1 = 'BSIG' MOP MAP ;

* Solution avec FSUR et vecteur constant = - p VECT_N
V1 = -1. 1. ; V1 = V1 / ('NORM' V1) ;
FS1 = 'FSUR' 'COQU' MOD1 V1 ;
* Solution avec FSUR et chpoint constant de composantes = - p VECT_N
ch_fy2 = 'MANU' 'CHPO' L2 2 'FX' (-1. * (SIN 45.)) 'FY' (COS 45.) ;
FS2 = 'FSUR' 'COQU' MOD1 ch_fy2 ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fy2 'FORC' 1. ;
  'TRAC' v_zz L2
     'TITR' 'TEST 1 - COQ2 - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr1 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) L2
     'TITR' 'TEST 1 - COQ2 - Comparaison PRESS (rouge) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr1 = PR1 + 1.E-15 ;
ERR1_1 = 'MAXI' 'ABS' (fS1-pr1) ;
'MESS' 'TEST 1-1 : ECART.ABS = ' ERR1_1 ;
ERR1_1 = 'MAXI' 'ABS' ( (fs1-pr1) / Er_pr1 ) ;
'MESS' 'TEST 1-1 : ECART.REL = ' ERR1_1 ;
'SI' (ERR1_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
ERR1_2 = 'MAXI' (fS2-pr1) 'ABS' ;
'MESS' 'TEST 1-2 : ECART.ABS = ' ERR1_2 ;
ERR1_2 = 'MAXI' 'ABS' ( (fs2-pr1) / Er_pr1 ) ;
'MESS' 'TEST 1-2 : ECART.REL = ' ERR1_2 ;
'SI' (ERR1_2 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;

* TEST 2 : MASSIF avec 1 type EF a bord droit (horizontal/vertical)
* ======== Comparaison de la solution avec PRESS
'OPTION' 'ELEM' 'QUA4' 'MODE' 'PLAN' 'CONT' ;
S1 = L1 'TRANS' 8 P2 ;

MOD1 = 'MODELISER' S1 'MECANIQUE' 'ELASTIQUE' ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' L1 'CHARGEMENT' 'PRESSION' ;
MAP = 'PRES' MOP 'PRES' 2. ;
PR1 = 'BSIG' MOP MAP ;
* Solution avec FSUR et vecteur constant = - p VECT_N (VECT_N = -y)
FS1 = 'FSUR' 'MASS' MOD1 L1 (0. 2.) ;
* Solution avec FSUR et chpoint constant de composantes = + p y
** ch_fy2 = 'MANU' 'CHPO' L1 1 'FX' 0. 'FY' 2. ;
ch_fy2 = 'MANU' 'CHPO' L1 1 'FY' 2. ;
FS2 = 'FSUR' 'MASS' MOD1 ch_fy2 ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fy2 'FORC' 1. ;
  'TRAC' v_zz S1
     'TITR' 'TEST 2 - MASS - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr1 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) S1
     'TITR' 'TEST 2 - MASS - Comparaison PRESS (roug) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr1 = PR1 + 1.E-15 ;
ERR2_1 = 'MAXI' 'ABS' (fs1-pr1) ;
'MESS' 'TEST 2-1 : ECART.ABS = ' ERR2_1 ;
ERR2_1 = 'MAXI' 'ABS' ( (fs1-pr1) / Er_pr1 ) ;
'MESS' 'TEST 2-1 : ECART.REL = ' ERR2_1 ;
ERR2_2 = 'MAXI' (fS2-pr1) 'ABS' ;
'MESS' 'TEST 2-2 : ECART.ABS = ' ERR2_2 ;
ERR2_2 = 'MAXI' 'ABS' ( (fs2-pr1) / Er_pr1 ) ;
'MESS' 'TEST 2-2 : ECART.REL = ' ERR2_2 ;
'SI' (ERR2_2 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;

* TEST 3 : MASSIF avec 2 types EF a bord courbe (arc de cercle)
* ======== Comparaison de la solution avec PRESS
* Nota : Il faut que les elements des bords generent 2 types d'element
*        dans le maillage suite au SURF !
'OPTION' 'ELEM' 'QUA4' 'MODE' 'PLAN' 'GENE' ;
C1 = 'CERC' 25 P1 P0 P2 ;
S2 = 'SURF' (L1 'ET' C1 'ET' ('DROITE' 8 P2 P0)) ;
'SI' IGRAPH ;
  'TRACER' S2 'QUAL' 'TITR' 'TEST3 - DISQUE PLAN - MASSIF' ;
'FINSI' ;

MOD2 = 'MODELISER' S2 'MECANIQUE' 'ELASTIQUE' 'DPGE' P0 ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' C1 'CHARGEMENT' 'PRESSION' ;
MAP = 'PRES' MOP 'PRES' 1. ;
PR2 = 'BSIG' MOP MAP ;

* Solution avec FSUR et chpoint variable decrivant - p VECT_N
x2 y2 = 'COORD' C1 ;
r = (x2*x2) + (y2*y2) ; r = r**(-0.5) ;
x2 = 'NOMC' 'FX' x2 ; y2 = 'NOMC' 'FY' y2 ;
**ch_fs = -1. * r * (x2 'ET' y2) ;
ch_fs = -0.1 * (x2 'ET' y2) ;
FS2 = 'FSUR' 'MASS' MOD2 ch_fs ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fs 'FORC' 1. ;
  'TRAC' v_zz S2
     'TITR' 'TEST 3 - MASS - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr2 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) S2
     'TITR' 'TEST 3 - MASS - Comparaison PRESS (roug) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr2 = PR2 + 1.E-15 ;
ERR3_1 = 'MAXI' (fS2-pr2) 'ABS' ;
'MESS' 'TEST 3-1 : ECART.ABS = ' ERR3_1 ;
ERR3_1 = 'MAXI' 'ABS' ( (fs2-pr2) / Er_pr2 ) ;
'MESS' 'TEST 3-1 : ECART.REL = ' ERR3_1 ;
'SI' (ERR3_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
*Nota : L'ecart absolu ERR3_1 "important" s'explique tres bien.
*       A l'utilisateur de reflechir !

* TEST 4 : Axisymetrie : disque horizontal
* ======== Comparaison de la solution avec PRESS
'OPTI' 'DIME' 2 'MODE' 'AXIS' ;
 
MOD1 = 'MODELISER' L1 'MECANIQUE' 'ELASTIQUE' 'COQ2' ;
MAT1 = 'MATERIAU' MOD1 'YOUN' 200.E+9 'NU' 0.3 'EPAI' 0.1 ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' L1 'CHARGEMENT' 'PRESSION' 'COQ2' ;
MAP = 'PRES' MOP 'PRES' 1. ;
PR1 = 'BSIG' MOP MAP ;
* Solution avec FSUR et vecteur constant = - p VECT_N (VECT_N = -z)
FS1 = 'FSUR' 'COQU' MOD1 (0. 1.) ;
* Solution avec FSUR et chpoint constant de composantes = + p z
ch_fz1 = 'MANU' 'CHPO' L1 2 'FZ' 1. 'FR' 0. ;
FS2 = 'FSUR' 'COQU' MOD1 ch_fz1 ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fz1 'FORC' 1. ;
  'TRAC' v_zz L1
     'TITR' 'TEST 4 - COQ2 - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr1 'FORC' 0.1 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 0.1 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) L1
     'TITR' 'TEST 4 - COQ2 - Comparaison PRESS (roug) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr1 = PR1 + 1.E-15 ;
ERR4_1 = 'MAXI' (fS1-pr1) 'ABS' ;
'MESS' 'TEST 4-1 : ECART.ABS = ' ERR4_1 ;
ERR4_1 = 'MAXI' 'ABS' ( (fs1-pr1) / Er_pr1 ) ;
'MESS' 'TEST 4-1 : ECART.REL = ' ERR4_1 ;
'SI' (ERR4_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
ERR4_2 = 'MAXI' (fS2-pr1) 'ABS' ;
'MESS' 'TEST 4-2 : ECART.ABS = ' ERR4_2 ;
ERR4_2 = 'MAXI' 'ABS' ( (fs2-pr1) / Er_pr1 ) ;
'MESS' 'TEST 4-2 : ECART.REL = ' ERR4_2 ;
'SI' (ERR4_2 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;

*----------------------------------------------------------------------*
* VERIFICATION EN DIMENSION 3
*----------------------------------------------------------------------*
'OPTI' 'DIME' 3 ;

* TEST 5 : COQ4 plaque plane
* ======== Comparaison de la solution avec PRESS
MOD1 = 'MODELISER' S1 'MECANIQUE' 'ELASTIQUE' 'COQ4' ;
MAT1 = 'MATERIAU' MOD1 'YOUN' 200.E+9 'NU' 0.3 'EPAI' 0.1 ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' S1 'CHARGEMENT' 'PRESSION' 'COQ4' ;
MAP = 'PRES' MOP 'PRES' 1. ;
PR1 = 'BSIG' MOP MAP ;
* Solution avec FSUR et vecteur constant = - p VECT_N (VECT_N = -z)
FS1 = 'FSUR' 'COQU' MOD1 (0. 0. 1.) MAT1 ;
* Solution avec FSUR et chpoint constant de composantes = + p z
ch_fz1 = 'MANU' 'CHPO' S1 3 'FX' 0. 'FY' 0. 'FZ' 1. ;
FS2 = 'FSUR' 'COQU' MOD1 ch_fz1 MAT1 ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fz1 'FORC' 1. ;
  'TRAC' v_zz S1
     'TITR' 'TEST 5 - COQ4 - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr1 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) S1
     'TITR' 'TEST 5 - COQ4 - Comparaison PRESS (rouge) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr1 = PR1 + 1.E-15 ;
ERR5_1 = 'MAXI' (fS1-pr1) 'ABS' ;
'MESS' 'TEST 5-1 : ECART.ABS = ' ERR5_1 ;
ERR5_1 = 'MAXI' 'ABS' ( (fs1-pr1) / Er_pr1 ) ;
'MESS' 'TEST 5-1 : ECART.REL = ' ERR5_1 ;
'SI' (ERR5_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
ERR5_2 = 'MAXI' (fS2-pr1) 'ABS' ;
'MESS' 'TEST 5-2 : ECART.ABS = ' ERR5_2 ;
ERR5_2 = 'MAXI' 'ABS' ( (fs2-pr1) / Er_pr1 ) ;
'MESS' 'TEST 5-2 : ECART.REL = ' ERR5_2 ;
'SI' (ERR5_2 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;

* TEST 6 : MASSIF 3D - CUBE
* ======== Comparaison de la solution avec PRESS
'OPTION' 'ELEM' 'CUB8' ;
V1 = 'VOLU' S1 'TRANS' 3 (0. 0. 20.) ;
F3 = 'FACE' 2 V1 ;

MOD1 = 'MODELISER' V1 'MECANIQUE' 'ELASTIQUE' ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' F3 'CHARGEMENT' 'PRESSION' ;
MAP = 'PRES' MOP 'PRES' -1. ;
PR1 = 'BSIG' MOP MAP ;
Er_pr1 = PR1 + 1.E-15 ;
* Solution avec FSUR et vecteur constant = - p VECT_N (VECT_N = +z)
FS1 = 'FSUR' 'MASS' MOD1 F3 (0. 0. -1.) ;
* Solution avec FSUR et chpoint constant de composantes = - p z
ch_fy1 = 'MANU' 'CHPO' F3 1 'FZ' -1. ;
FS2 = 'FSUR' 'MASS' MOD1 ch_fy1 ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fy1 'FORC' 1. ;
  'TRAC' v_zz V1
     'TITR' 'TEST 6 - MASS - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr1 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) V1
     'TITR' 'TEST 6 - MASS - Comparaison PRESS (rouge) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
ERR6_1 = 'MAXI' (fS1-pr1) 'ABS' ;
'MESS' 'TEST 6-1 : ECART.ABS = ' ERR6_1 ;
ERR6_1 = 'MAXI' 'ABS' ( (fs1-pr1) / Er_pr1 ) ;
'MESS' 'TEST 6-1 : ECART.REL = ' ERR6_1 ;
'SI' (ERR6_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
ERR6_2 = 'MAXI' (fS2-pr1) 'ABS' ;
'MESS' 'TEST 6-2 : ECART.ABS = ' ERR6_2 ;
ERR6_2 = 'MAXI' 'ABS' ( (fs2-pr1) / Er_pr1 ) ;
'MESS' 'TEST 6-2 : ECART.REL = ' ERR6_2 ;
'SI' (ERR6_2 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;

* TEST 7 : MASSIF avec 2 types EF a bord courbe (portion de cylindre)
* ======== Comparaison de la solution avec PRESS
vecz = 0. 0. 17. ;
V2 = 'VOLU' S2 'TRANS' 8 vecz ;
PtV2 = 'POINT' V2 'CYLI' P0 (P0 'PLUS' vecz) P2 1.E-8 ;
S2b = 'ELEM' ('ENVE' V2) 'APPU' 'STRI' PtV2 ;

MOD2 = 'MODELISER' V2 'MECANIQUE' 'ELASTIQUE' ;

* Solution avec PRESS pression constante = p
MOP = 'MODE' S2b 'CHARGEMENT' 'PRESSION' ;
MAP = 'PRES' MOP 'PRES' 1. ;
PR2 = 'BSIG' MOP MAP ;
* Solution avec FSUR et chpoint variable decrivant - p VECT_N
x2 y2 z2 = 'COORD' S2b ;
r = (x2*x2) + (y2*y2) ; r = r**(-0.5) ;
x2 = 'NOMC' 'FX' x2 ; y2 = 'NOMC' 'FY' y2 ; z2 = 'NOMC' 'FZ' z2 ;
ch_fs = -1. * r * (x2 'ET' y2) ;
**ch_fs = -0.1 * (x2 'ET' y2) ;
FS2 = 'FSUR' 'MASS' MOD2 ch_fs ;
'SI' IGRAPH ;
  v_zz = 'VECT' ch_fs 'FORC' 1. ;
  'TRAC' v_zz ('ARETE' V2)
     'TITR' 'TEST 7 - MASS - CHPOINT DE Forces SURfaciques' ;
  v_z1 = 'VECT' pr2 'FORC' 1. 'ROUG' ;
  v_z2 = 'VECT' fs2 'FORC' 1. 'BLEU' ;
  'TRAC' (v_z1 'ET' v_z2) ('ARETE' V2)
     'TITR' 'TEST 7 - MASS - Comparaison PRESS (roug) / FSUR (bleu)' ;
'FINSI' ;
* Tests simples de bon fonctionnement
Er_pr2 = PR2 + 1.E-15 ;
ERR7_1 = 'MAXI' (fS2-pr2) 'ABS' ;
'MESS' 'TEST 7-1 : ECART.ABS = ' ERR7_1 ;
ERR7_1 = 'MAXI' 'ABS' ( (fs2-pr2) / Er_pr2 ) ;
'MESS' 'TEST 7-1 : ECART.REL = ' ERR7_1 ;
'SI' (ERR7_1 '>EG' 1.E-11) ; 'ERREUR' 5 ; 'FINSI' ;
*Nota : L'ecart absolu ERR7_1 "important" s'explique tres bien.
*       A l'utilisateur de reflechir !

'FIN' ;