Télécharger tens2.dgibi

Retour à la liste

* fichier :  tens2.dgibi
'OPTI' echo 0 ;
************************************************************************
************************************************************************
************************************************************************
* NOM         : TENS2
* DESCRIPTION : Test de la décomposition principale, en particulier
*               les valeurs propres et la recomposition.
*               On teste les cas où les valeurs propres sont proches.
*               TENS PRIN et TENS RECO
*               D'autres options de TENS sont egalement testees
*               'INVE' 'DET' pour les matrices de rotation
*               (orthonormales non symetriques)
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
*               mel : stephane.gounand@cea.fr
**********************************************************************
* VERSION    : v1, 09/10/2025, version initiale
* HISTORIQUE : v1, 09/10/2025, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
*
interact=faux ;
*
* erreur relative entre deux champs normee point par point
*
'DEBP' errchpor ch1*'CHPOINT' ch2*'CHPOINT' ;
   nc1 = 'EXTR' ch1 'COMP' ;
   nc2 = 'EXTR' ch2 'COMP' ;
   nct = 'UNIQ' (nc1 'ET' nc2) ;
   first = vrai ;
   'REPE' iict ('DIME' nct) ;
      ict = &iict ;
      nocomp = 'EXTR' nct ict ;
      ach1i = 'ABS' ('EXCO' nocomp ch1 'SCAL') ;
      'SI' first ;
         chmaxa = ach1i ;
         first = faux ;
      'SINO' ;
         chg  = 'MASQ' chmaxa 'EGSUPE' ach1i ;
         chd  = '-' 1. chg ;
         chmaxa = (chg '*' chmaxa) '+' (chd '*' ach1i) ;
      'FINS' ;
      ach2i = 'ABS' ('EXCO' nocomp ch2 'SCAL') ;
      chg  = 'MASQ' chmaxa 'EGSUPE' ach2i ;
      chd  = '-' 1. chg ;
      chmaxa = (chg '*' chmaxa) '+' (chd '*' ach2i) ;
   'FIN' iict ;
* Seuillage
   vpetit = 'VALE' 'PETIT' ;
   chg  = 'MASQ' chmaxa 'EGSUPE' vpetit ;
   chd  = '-' 1. chg ;
   chmaxa = (chg '*' chmaxa) '+' (chd '*' vpetit) ;
   dc12 = '/' ('ABS' ('-' ch1 ch2)) chmaxa ;
'FINP' dc12 ;
*
* Test si les vecteurs propres sont bien des vecteurs propres
*
'DEBP' tvectp ;
'ARGU' tensym*'CHPOINT' ;
'ARGU' chvecp*'CHPOINT' ;
'ARGU' tolrel*'FLOTTANT' ;
lok = vrai ;
amax = 'TENS' 'NORMINF' tensym ;
amax = 'BORN' amax 'MINIMUM' ('VALE' 'PETIT') ;
tol = amax '*' tolrel ;
* Passage CHPOINT -> TABLE Tenseur symetrique
ttens = 'TABL' ;
vdim = 'VALE' 'DIME' ;
'REPE' ii vdim ;
   i = &ii ;
   ttens . i = 'TABL' ;
   'REPE' ij i ;
      j = &ij ;
      ngij = 'CHAI' 'G' i j ;
      ttens . i . j = 'EXCO' ngij tensym ;
      'SI' ('NEG' i j);
         ttens . j . i = ttens . i . j ;
      'FINS' ;
   'FIN' ij ;
'FIN' ii ;
* Passage CHPOINT -> TABLE Vecteurs propres
tvecp = 'TABL' ;
vdim = 'VALE' 'DIME' ;
'REPE' ii vdim ;
   i = &ii ;
   tvecp . i = 'TABL' ;
   'REPE' ij vdim ;
      j = &ij ;
      ngij = 'CHAI' 'CO' j i ;
      tvecp . i . j = 'EXCO' ngij chvecp ;
   'FIN' ij ;
'FIN' ii ;
* Passage CHPOINT -> TABLE Valeurs propres
tvalp = 'TABL' ;
vdim = 'VALE' 'DIME' ;
'REPE' ii vdim ;
   i = &ii ;
   tvalp . i = 'TABL' ;
   ngii = 'CHAI' 'SI' i i ;
   tvalp . i = 'EXCO' ngii chvecp ;
'FIN' ii ;
* Verif Vecteurs propres
tprod = 'TABL' ;
'REPE' ii vdim ;
   i = &ii ;
   tprod . i = 'TABL' ;
   'REPE' ij vdim ;
      j = &ij ;
      chprod = 'VIDE' 'CHPOINT'/'DIFFUS' ;
      'REPE' ik vdim ;
         k = &ik ;
         'SI' ('EGA' j k) ;
            inc = '*' ('-' ttens . k . j tvalp . i) tvecp . i . k ;
         'SINO' ;
            inc = '*' ttens . k . j tvecp . i . k ;
         'FINS' ;
         chprod = chprod '+' inc ;
      'FIN' ik ;
      tprod . i . j = chprod ;
   'FIN' ij ;
'FIN' ii ;
'REPE' ii vdim ;
   i = &ii ;
   'REPE' ij vdim ;
      j= &ij ;
      chprod = tprod . i . j ;
      vprod  = 'MASQ' ('ABS' chprod) 'INFERIEUR' tol ;
      nprod = 'MINI' vprod ;
      tnprod = '>' nprod 0. ;
      lok = lok 'ET' tnprod ;
      'SI' ('NON' tnprod) ;
         'MESS' '!!! Le vecteur i=' i ' nest pas un vecteur propre' ;
*         'MESS' 'par sa composante j=' j ' dont la norme' ' ' nprod  ' >' ' ' tol ;
         'MESS' 'par sa composante j=' j ;
         nocor = 'POIN' vprod 'MINI' ;
         'MESS' '!!! Nb points a problemes = ' ('NBEL' nocor) ;
      'FINS' ;
   'FIN' ij ;
'FIN' ii ;
'SI' ('NON' lok) ; 'ERRE' stop ; 'FINS' ;
'FINP' lok ;
*
* Test si les colonnes d'une matrice sont bien orthonormales
* (matrice de rotation et matrice de vecteurs propres)
*
'DEBP' ttensrot ;
'ARGU' rottot*'CHPOINT' ;
'ARGU' tol*'FLOTTANT' ;
lok = vrai ;
*  Identite
rott1   = 'TENS' 'IDEN' rottot ;
drott1  = '-' rott1 rottot ;
ndrott1 = 'TENS' drott1 'NORMINF' ;
nndrott1 = 'MAXI' ndrott1 'ABS' ;
tdrott1 = '<' nndrott1 tol ;
lok = lok 'ET' tdrott1 ;
'SI' ('NON' tdrott1) ; 'MESS' '!!! Calcul identite' ' ' nndrott1 ' >' ' ' tol ; 'FINS' ;
* Determinant unite
drott2   = '-' ('TENS' 'DET' rottot) 1. ;
nndrott2 = 'MAXI' drott2 'ABS' ;
tdrott2 = '<' nndrott2 tol ;
lok = lok 'ET' tdrott2 ;
'SI' ('NON' tdrott2) ; 'MESS' '!!! Determinant non unite' ' ' nndrott2  ' >' ' ' tol ; 'FINS' ;
* Verif orthonormalite
trot = 'TABL' ;
vdim = 'VALE' 'DIME' ;
'REPE' ii vdim ;
   i = &ii ;
   trot . i = 'TABL' ;
   'REPE' ij vdim ;
      j = &ij ;
      ngij = 'CHAI' 'CO' i j ;
      trot . i . j = 'EXCO' ngij rottot ;
   'FIN' ij ;
'FIN' ii ;
* Pscal des colonnes de trot
tpsca = 'TABL' ;
'REPE' ij vdim ;
   j = &ij ;
   tpsca . j = 'TABL' ;
   'REPE' ik vdim ;
      k = &ik ;
      pjk = 'VIDE' 'CHPOINT'/'DISCRET' ;
      'REPE' ii vdim ;
         i = &ii ;
         pjk = pjk '+' (trot . i . j '*' trot . i . k) ;
      'FIN' ii ;
      tpsca . j . k = pjk ;
   'FIN' ik ;
'FIN' ij ;
* Verif trot orthonormale
'REPE' ij vdim ;
   j = &ij ;
   'REPE' ik vdim ;
      k = &ik ;
      'SI' ('EGA' j k) ;
         dpsca  = (tpsca . j . k) '-' 1.d0 ;
      'SINO' ;
         dpsca  = tpsca . j . k ;
      'FINS' ;
      ddpsca = 'MAXI' dpsca 'ABS' ;
      tdpsca = '<' ddpsca tol ;
      lok = lok 'ET' tdpsca ;
      'SI' ('NON' tdpsca) ;
         'MESS' '!!! Matrice non orthonormale' ' ' ddpsca  ' >' ' ' tol ;
         'MESS' 'colonne' ' ' j '.pscal.colonne' ' '  k '=' ddpsca ;
      'FINS' ;
   'FIN' ik ;
'FIN' ij ;
* Transpose = inverse
trottot = 'TENS' 'TRANS' rottot ;
irottot = 'TENS' 'INVE' rottot ;
drott3  = '-' trottot irottot ;
ndrott3 = 'TENS' drott3 'NORMINF' ;
nndrott3 = 'MAXI' ndrott3 'ABS' ;
tdrott3 = '<' nndrott3 tol ;
lok = lok 'ET' tdrott3 ;
'SI' ('NON' tdrott3) ; 'MESS' '!!! Inverse.NE.Transpose' ' ' nndrott3 ' >' ' ' tol ; 'FINS' ;
* Inverse correcte
tirot = 'TABL' ;
vdim = 'VALE' 'DIME' ;
'REPE' ii vdim ;
   i = &ii ;
   tirot . i = 'TABL' ;
   'REPE' ij vdim ;
      j = &ij ;
      ngij = 'CHAI' 'CO' i j ;
      tirot . i . j = 'EXCO' ngij irottot ;
   'FIN' ij ;
'FIN' ii ;
* Produit
vdim = 'VALE' 'DIME' ;
*tprod = 'TABL'  ;
'REPE' ij vdim ;
   j = &ij ;
*   tprod . j = 'TABL' ;
   'REPE' ik vdim ;
      k = &ik ;
      ngjk = 'CHAI' 'CO' j k ;
      gjk = 'VIDE' 'CHPOINT'/'DIFFUS' ;
      'REPE' ii vdim ;
         i = &ii ;
         gjk = gjk '+' (trot . j . i  '*' trot . i . k) ;
      'FIN' ii ;
*      tprod . j .k = gjk ;
      'SI' ('EGA' j k) ;
         vref = 1. ;
      'SINO' ;
         vref = 0. ;
      'FINS' ;
      drott4   = '-' gjk vref ;
      nndrott4 = 'MAXI' drott2 'ABS' ;
      tdrott4 = '<' nndrott4 tol ;
      lok = lok 'ET' tdrott4 ;
      'SI' ('NON' tdrott4) ; 'MESS' '!!! M*M^-1 non unite' ' j=' j ' k =' k '   ' nndrott4  ' >' ' ' tol ; 'FINS' ;
   'FIN' ik ;
'FIN' ij ;
'SI' ('NON' lok) ; 'ERRE' stop ; 'FINS' ;
'FINP' lok ;
*
lok = vrai ;
pet = 'VALE' 'PETIT' ; lpet = 'LOG' pet ;
pre = 'VALE' 'PREC' ;  lpre = 'LOG' pre ;
*
tol      = '*' pre 10. ;
tolrot   = '*' pre 10. ;
tolvecp  = '*' pre 30. ;
tolrec3d = '*' pre 30. ;
*
'OPTION' 'DIME' 1 'ELEM' 'SEG2' ;
'OPTI' 'MODE' unid plan ;
'MESS' ' ********** ' ;
'MESS' ' *   1D   * ' ;
'MESS' ' ********** ' ;
*
nx = 100000 ;
*nx = 3000000 ;
*nx = 1 ;
p1 = 'POIN' 0. ; p2 = 'POIN' 1. ;
el = 'DROI' nx p1 p2 ;
mod = 'MODE' el 'MECANIQUE' ;
* Scaling global des termes de la matrice
puig = 'BRUI' 'BLAN' 'UNIF' 0. ('/' lpet -2.) el ;
epsg  = 'EXP' puig ;
*
d1 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
d1 = 'NOMC' 'SI11' d1 'NATURE' 'DIFFUS' ;
v1 = 'MANU' 'CHPO' el 1 'CO11' 1. 'NATURE' 'DIFFUS' ;
ltens = ttensrot v1 tolrot ;
lok = lok 'ET' ltens ;
'SI' ('NON' ltens) ;
   'MESS' '!!! Test tenseur de rotation not ok' ;
'SINO' ;
   'MESS' 'Test tenseur de rotation OK' ;
'FINS' ;
dec = d1 'ET' v1 ;
*
cha = 'TENS' 'RECO' dec ;
*
* Décomposition spectrale
*
decb = 'TENS' 'PRIN' cha ;
*  Seulement les vp
decbb = 'TENS' 'VALP' cha ;
legvp = '<' ('MAXI' ('-' ('EXCO' ('MOTS' 'SI11') decb) decbb) 'ABS') pet ;
lok = lok 'ET' legvp ;
'SI' ('NON' legvp) ;
   'MESS' '!!! Test valeurs propres egales not ok' ;
'SINO' ;
   'MESS' 'Test valeurs propres egales OK' ;
'FINS' ;
 
decc = 'EXCO' 'CO11' decb ;
ltens = ttensrot decc tolrot ;
lok = lok 'ET' ltens ;
'SI' ('NON' ltens) ;
   'MESS' '!!! Test vecteurs propres orthonormaux not ok' ;
'SINO' ;
   'MESS' 'Test vecteurs propres orthonormaux OK' ;
'FINS' ;
lvect = tvectp cha decb tolvecp ;
'SI' ('NON' lvect) ;
   'MESS' '!!! Test vecteurs propres not ok' ;
'SINO' ;
   'MESS' 'Test vecteurs propres OK' ;
'FINS' ;
*
chab = 'TENS' 'RECO' decb ;
d1b  = 'EXCO' 'SI11' decb 'SI11' ;
* Test valeur propre
dd1 = 'MAXI' ('-' d1 d1b) 'ABS' ;
tdd1 = '<' dd1 tol ;
lok = lok 'ET' tdd1 ;
'SI' ('NON' tdd1) ; 'MESS' '!!! Calcul err valp 1' ' ' dd1 ' >' ' '  tol ; 'FINS' ;
* Test recomposition matrice
dcha = '-' chab cha ;
ndcha = 'MAXI' ('TENS' dcha 'NORMINF') 'ABS' ;
tdcha = '<' ndcha tol ;
lok = lok 'ET' tdcha ;
'SI' ('NON' tdcha) ; 'MESS' '!!! Calcul recomposition' ' ' ndcha  ' >' ' ' tol ; 'FINS' ;
'SI' lok ; 'MESS' 'OK' ; 'FINS' ;
*
'OPTION' 'DIME' 2 'ELEM' 'TRI3' ;
'OPTI' 'MODE' plan defo ;
'MESS' ' ********** ' ;
'MESS' ' *   2D   * ' ;
'MESS' ' ********** ' ;
*
'REPE' iikas 2 ;
   ikas = &iikas ;
* Scaling global des termes de la matrice
   puig = 'BRUI' 'BLAN' 'UNIF' 0. ('/' lpet -2.) el ;
   epsg  = 'EXP' puig ;
   vp1 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
* Valeurs propres distinctes
   'SI' ('EGA' ikas 1) ;
      'MESS' 'Valeurs propres distinctes' ;
      vp2 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
   'FINS' ;
* Valeurs propres proches
   'SI' ('EGA' ikas 2) ;
      'MESS' 'Valeurs propres proches' ;
      puin = 'ENTI' 'PROC' ('BRUI' 'BLAN' 'UNIF' lpre ('*' lpre -1.d0) el) ;
      eps  = 'EXP' puin ;
      vp2 = vp1 '+' (('BRUI' 'BLAN' 'UNIF' 0. 1. el) '*' eps '*' epsg) ;
   'FINS' ;
* Les plus grandes valeurs dans d1, les plus petites dans d2
   m1 = 'MASQ' vp1 'EGSUPE' vp2 ;
   m2 = '-' 1. m1 ;
   d1 = ('*' vp1 m1) '+' ('*' vp2 m2) ;
   d1 = 'NOMC' 'SI11' d1 'NATURE' 'DIFFUS' ;
   d2 = ('*' vp1 m2) '+' ('*' vp2 m1) ;
   d2 = 'NOMC' 'SI22' d2 'NATURE' 'DIFFUS' ;
*
   ang = 'BRUI' 'BLAN' 'UNIF' 0. 180. el ;
   ca  = 'COS' ang ; sa = 'SIN' ang ; msa = sa '*' -1. ;
*
   rot11 = 'NOMC' 'CO11' ca 'NATURE' 'DIFFUS' ;
   rot21 = 'NOMC' 'CO21' sa 'NATURE' 'DIFFUS' ;
   rot12 = 'NOMC' 'CO12' msa 'NATURE' 'DIFFUS' ;
   rot22 = 'NOMC' 'CO22' ca 'NATURE' 'DIFFUS' ;
   rottot = rot11 'ET' rot21 'ET' rot12 'ET' rot22 ;
* Vérif tenseur de rotation
   ltens = ttensrot rottot tolrot ;
   lok = lok 'ET' ltens ;
   'SI' ('NON' ltens) ;
      'MESS' '!!! Test tenseur de rotation not ok' ;
   'SINO' ;
      'MESS' 'Test tenseur de rotation OK' ;
   'FINS' ;
   vprot = d1 'ET' d2 'ET' rottot ;
*
   vdim = 'VALE' 'DIME' ;
   cha = 'VIDE' 'CHPOINT'/'DIFFUS' ;
   'REPE' ij vdim ;
      j = &ij ;
      'REPE' ik j ;
         k = &ik ;
         ngjk = 'CHAI' 'G' j k ;
         gjk = 'VIDE' 'CHPOINT'/'DIFFUS' ;
         'REPE' ii vdim ;
            i = &ii ;
            nqji = 'CHAI' 'CO' i j ;
            nli  = 'CHAI' 'SI' i i ;
            nqik = 'CHAI' 'CO' i k ;
            qji  = 'EXCO' vprot nqji ;
            li   = 'EXCO' vprot nli ;
            qik  = 'EXCO' vprot nqik ;
            gjk  = gjk '+' (qji '*' li '*' qik) ;
         'FIN' ii ;
         gjk = 'NOMC' ngjk gjk 'NATURE' 'DIFFUS' ;
         cha = cha 'ET' gjk ;
      'FIN' ik ;
   'FIN' ij ;
*
* Décomposition spectrale
*
   decb = 'TENS' 'PRIN' cha ;
*  Seulement les vp
   decbb = 'TENS' 'VALP' cha ;
   legvp = '<' ('MAXI' ('-' ('EXCO' ('MOTS' 'SI11' 'SI22') decb) decbb) 'ABS') pet ;
   lok = lok 'ET' legvp ;
   'SI' ('NON' legvp) ;
      'MESS' '!!! Test valeurs propres egales not ok' ;
   'SINO' ;
      'MESS' 'Test valeurs propres egales OK' ;
   'FINS' ;
*
   decc = 'EXCO' ('MOTS' 'CO11' 'CO12' 'CO21' 'CO22') decb ;
   ltens = ttensrot decc tolrot ;
   lok = lok 'ET' ltens ;
   'SI' ('NON' ltens) ;
      'MESS' '!!! Test vecteurs propres orthonormaux not ok' ;
   'SINO' ;
      'MESS' 'Test vecteurs propres orthonormaux OK' ;
   'FINS' ;
   lvect = tvectp cha decb tolvecp ;
   'SI' ('NON' lvect) ;
      'MESS' '!!! Test vecteurs propres not ok' ;
   'SINO' ;
      'MESS' 'Test vecteurs propres OK' ;
   'FINS' ;
*
   chab = 'TENS' 'RECO' decb ;
   d1b  = 'EXCO' 'SI11' decb 'SI11' ;
   d2b  = 'EXCO' 'SI22' decb 'SI22' ;
* Test valeurs propres
   dd1 = 'NOMC' 'SCAL' ('ABS' ('-' d1 d1b)) ;
   mdd1 = 'MASQ' dd1 'INFERIEUR' ('*' tol epsg) ;
   mmdd1 = 'MINI' mdd1 ;
   tdd1 = '>' mmdd1 0. ;
   lok = lok 'ET' tdd1 ;
   'SI' ('NON' tdd1) ; 'MESS' '!!! Calcul valp1' ; 'FINS' ;
   dd2 = 'NOMC' 'SCAL' ('ABS' ('-' d2 d2b)) ;
   mdd2 = 'MASQ' dd2 'INFERIEUR' ('*' tol epsg) ;
   mmdd2 = 'MINI' mdd2 ;
   tdd2 = '>' mmdd2 0. ;
   lok = lok 'ET' tdd2 ;
   'SI' ('NON' tdd2) ; 'MESS' '!!! Calcul valp2' ; 'FINS' ;
* Test recomposition matrice
   dcha  = 'TENS' 'NORMINF' ('-' chab cha) ;
   mdcha = 'MASQ' dcha 'INFERIEUR' ('*' tol epsg) ;
   mmdcha = 'MINI' mdcha ;
   tdcha = '>' mmdcha 0. ;
   lok = lok 'ET' tdcha ;
   'SI' ('NON' tdcha) ; 'MESS' '!!! Calcul recomposition' ; 'FINS' ;
   'SI' lok ; 'MESS' 'OK' ; 'SINO' ; 'ERRE' stop ; 'FINS' ;
'FIN' iikas ;
*
*
*
'OPTION' 'DIME' 3 'ELEM' 'TET4' ;
'OPTI' 'MODE' 'TRID' ;
'MESS' ' ********** ' ;
'MESS' ' *   3D   * ' ;
'MESS' ' ********** ' ;
*
'REPE' iikas 3 ;
   ikas = &iikas ;
* Scaling global des termes de la matrice
   puig = 'BRUI' 'BLAN' 'UNIF' 0. ('/' lpet -2.) el ;
   epsg  = 'EXP' puig ;
   vp1 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
*
   'SI' ('EGA' ikas 1) ;
      'MESS' 'Valeurs propres distinctes' ;
       vp2 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
       vp3 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
   'FINS' ;
   'SI' ('EGA' ikas 2) ;
      'MESS' 'Deux valeurs propres proches' ;
      vp2 = '*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) epsg ;
      puin = 'ENTI' 'PROC' ('BRUI' 'BLAN' 'UNIF' lpre ('*' lpre -1.d0) el) ;
      eps  = 'EXP' puin ;
      vp3 = vp2 '+' ('*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) (eps '*' epsg)) ;
   'FINS' ;
   'SI' ('EGA' ikas 3) ;
      'MESS' 'Trois valeurs propres proches' ;
      puin = 'ENTI' 'PROC' ('BRUI' 'BLAN' 'UNIF' lpre ('*' lpre -1.d0) el) ;
      eps  = 'EXP' puin ;
      vp2 = vp1 '+' ('*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) (eps '*' epsg)) ;
      vp3 = vp1 '+' ('*' ('BRUI' 'BLAN' 'UNIF' 0. 1. el) (eps '*' epsg)) ;
   'FINS' ;
* Les plus grandes valeurs dans d1, les plus petites dans d3, les autres au milieu
   m12 = 'MASQ' vp1 'EGSUPE' vp2 ; m21 = 1. '-' m12 ;
   m23 = 'MASQ' vp2 'EGSUPE' vp3 ; m32 = 1. '-' m23 ;
   m13 = 'MASQ' vp1 'EGSUPE' vp3 ; m31 = 1. '-' m13 ;
* Pour trouver ces formules, on a énuméré les cas et utilisé le fait que le ET logique est la multiplication
* Plus grandes valeurs
   d1vp1 = m12 '*' m13 ; d1vp2 = m21 '*' m23 ; d1vp3 = m31 '*' m32 ;
* Plus petites valeurs
   d3vp1 = m21 '*' m31 ; d3vp2 = m12 '*' m32 ; d3vp3 = m13 '*' m23 ;
*
   d2vp1 = 1. '-' (d1vp1 '+' d3vp1) ;
   d2vp2 = 1. '-' (d1vp2 '+' d3vp2) ;
   d2vp3 = 1. '-' (d1vp3 '+' d3vp3) ;
*
   d1 =  ('*' d1vp1 vp1) '+' ('*' d1vp2 vp2) '+' ('*' d1vp3 vp3) ;
   d1 = 'NOMC' 'SI11' d1 'NATURE' 'DIFFUS' ;
   d2 =  ('*' d2vp1 vp1) '+' ('*' d2vp2 vp2) '+' ('*' d2vp3 vp3) ;
   d2 = 'NOMC' 'SI22' d2 'NATURE' 'DIFFUS' ;
   d3 =  ('*' d3vp1 vp1) '+' ('*' d3vp2 vp2) '+' ('*' d3vp3 vp3) ;
   d3 = 'NOMC' 'SI33' d3 'NATURE' 'DIFFUS' ;
*
* Génération d'une matrice de rotation quelconque
* Algo III.4 Graphics Gems III par James Arvo (1995)
*
   theta = 'BRUI' 'BLAN' 'UNIF' 0. 180. el ;
   phi   = 'BRUI' 'BLAN' 'UNIF' 0. 180. el ;
   zed   = 'BRUI' 'BLAN' 'UNIF' 0.5 0.5 el ;
* Rotation autour de l'axe z
   ct  = 'COS' theta ; st = 'SIN' theta ; mst = st '*' -1. ;
   trz = 'TABL' ;
   trz . 1 = 'TABL' ;
   trz . 1 . 1 = ct  ; trz . 1 . 2 = st ; trz . 1 . 3 = 0. ;
   trz . 2 = 'TABL' ;
   trz . 2 . 1 = mst ; trz . 2 . 2 = ct ; trz . 2 . 3 = 0. ;
   trz . 3 = 'TABL' ;
   trz . 3 . 1 = 0.  ; trz . 3 . 2 = 0. ; trz . 3 . 3 = 1. ;
* Matrice de Householder H = 2 vvt - I
   cp = 'COS' phi ; sp = 'SIN' phi ; sqz = '**' zed 0.5 ; msqz = '**' (1. '-' zed) 0.5 ;
   v1  = cp '*' sqz ; v2 = sp '*' sqz ; v3 = msqz ;
   th = 'TABL' ;
   th . 1 = 'TABL' ;
   th . 1 . 1 = 2. '*' v1 '*' v1 '-' 1. ; th . 1 . 2 = 2. '*' v1 '*' v2        ; th . 1 . 3 = 2. '*' v1 '*' v3 ;
   th . 2 = 'TABL' ;
   th . 2 . 1 = th . 1 . 2              ; th . 2 . 2 = 2. '*' v2 '*' v2 '-' 1. ; th . 2 . 3 = 2. '*' v2 '*' v3 ;
   th . 3 = 'TABL' ;
   th . 3 . 1 = th . 1 . 3              ; th . 3 . 2 = th . 2 . 3              ; th . 3 . 3 = 2. '*' v3 '*' v3 '-' 1. ;
* Le produit de H et R
   vdim = 'VALE' 'DIME' ;
   charot = 'VIDE' 'CHPOINT'/'DIFFUS' ;
   'REPE' ij vdim ;
      j = &ij ;
      'REPE' ik vdim ;
         k = &ik ;
         ngjk = 'CHAI' 'CO' j k ;
         gjk = 'VIDE' 'CHPOINT'/'DIFFUS' ;
         'REPE' ii vdim ;
            i = &ii ;
            gjk = gjk '+' (th . j . i  '*' trz . i . k) ;
         'FIN' ii ;
         gjk = 'NOMC' ngjk gjk 'NATURE' 'DIFFUS' ;
         charot = charot 'ET' gjk ;
      'FIN' ik ;
   'FIN' ij ;
* Vérif tenseur de rotation
   ltens = ttensrot charot tolrot ;
   lok = lok 'ET' ltens ;
   'SI' ('NON' ltens) ;
      'MESS' '!!! Test tenseur de rotation not ok' ;
   'SINO' ;
      'MESS' 'Test tenseur de rotation OK' ;
   'FINS' ;
   vprot = d1 'ET' d2 'ET' d3 'ET' charot ;
*dbg   mpn6 = 'MANU' 'POI1' ('NOEU' 6) ;
*dbg   vprot = 'REDU' vprot mpn6 ;
*dbg   d1 = 'REDU' d1 mpn6 ;
*dbg   d2 = 'REDU' d2 mpn6 ;
*dbg   d3 = 'REDU' d3 mpn6 ;
* Rotation de la matrice diagonale des vps
   cha = 'VIDE' 'CHPOINT'/'DIFFUS' ;
   'REPE' ij vdim ;
      j = &ij ;
      'REPE' ik j ;
         k = &ik ;
         ngjk = 'CHAI' 'G' j k ;
         gjk = 'VIDE' 'CHPOINT'/'DIFFUS' ;
         'REPE' ii vdim ;
            i = &ii ;
            nqji = 'CHAI' 'CO' i j ;
            nli  = 'CHAI' 'SI' i i ;
            nqik = 'CHAI' 'CO' i k ;
            qji  = 'EXCO' vprot nqji ;
            li   = 'EXCO' vprot nli ;
            qik  = 'EXCO' vprot nqik ;
            gjk  = gjk '+' (qji '*' li '*' qik) ;
         'FIN' ii ;
         gjk = 'NOMC' ngjk gjk 'NATURE' 'DIFFUS' ;
         cha = cha 'ET' gjk ;
      'FIN' ik ;
   'FIN' ij ;
*
* Décomposition spectrale
*
   decb = 'TENS' 'PRIN' cha ;
*  Seulement les vp
   decbb = 'TENS' 'VALP' cha ;
   decbr = 'EXCO' ('MOTS' 'SI11' 'SI22' 'SI33') decb ;
   ddecb = ERRCHPOR decbb decbr ;
   legvp = '<' ('MAXI' ddecb 'ABS') ('VALE' 'PREC') ;
   lok = lok 'ET' legvp ;
   'SI' ('NON' legvp) ;
      'MESS' '!!! Test valeurs propres egales not ok' ;
   'SINO' ;
      'MESS' 'Test valeurs propres egales OK' ;
   'FINS' ;
   decc = 'EXCO' ('MOTS' 'CO11' 'CO12' 'CO13' 'CO21' 'CO22' 'CO23' 'CO31' 'CO32' 'CO33') decb ;
*
   ltens = ttensrot decc tolrot ;
   lok = lok 'ET' ltens ;
   'SI' ('NON' ltens) ;
      'MESS' '!!! Test vecteurs propres orthonormaux not ok' ;
   'SINO' ;
      'MESS' 'Test vecteurs propres orthonormaux OK' ;
   'FINS' ;
   lvect = tvectp cha decb tolvecp ;
   'SI' ('NON' lvect) ;
      'MESS' '!!! Test vecteurs propres not ok' ;
   'SINO' ;
      'MESS' 'Test vecteurs propres OK' ;
   'FINS' ;
*
   chab = 'TENS' 'RECO' decb ;
   d1b  = 'EXCO' 'SI11' decb 'SI11' ;
   d2b  = 'EXCO' 'SI22' decb 'SI22' ;
   d3b  = 'EXCO' 'SI33' decb 'SI33' ;
* Test valeurs propres
   dd1 = 'NOMC' 'SCAL' ('ABS' ('-' d1 d1b)) ;
   mdd1 = 'MASQ' dd1 'INFERIEUR' ('*' tol epsg) ;
   mmdd1 = 'MINI' mdd1 ;
   tdd1 = '>' mmdd1 0. ;
   lok = lok 'ET' tdd1 ;
   'SI' ('NON' tdd1) ; 'MESS' '!!! Calcul valp1' ; 'FINS' ;
   dd2 = 'NOMC' 'SCAL' ('ABS' ('-' d2 d2b)) ;
   mdd2 = 'MASQ' dd2 'INFERIEUR' ('*' tol epsg) ;
   mmdd2 = 'MINI' mdd2 ;
   tdd2 = '>' mmdd2 0. ;
   lok = lok 'ET' tdd2 ;
   'SI' ('NON' tdd2) ; 'MESS' '!!! Calcul valp2' ; 'FINS' ;
   dd3 = 'NOMC' 'SCAL' ('ABS' ('-' d3 d3b)) ;
   mdd3 = 'MASQ' dd3 'INFERIEUR' ('*' tol epsg) ;
   mmdd3 = 'MINI' mdd3 ;
   tdd3 = '>' mmdd3 0. ;
   lok = lok 'ET' tdd3 ;
   'SI' ('NON' tdd3) ; 'MESS' '!!! Calcul valp3' ; 'FINS' ;
* Test recomposition matrice
   tolreco = tolrec3d ;
   dcha  = 'TENS' 'NORMINF' ('-' chab cha) ;
   mdcha = 'MASQ' dcha 'INFERIEUR' ('*' tolreco epsg) ;
   mmdcha = 'MINI' mdcha ;
   tdcha = '>' mmdcha 0. ;
   lok = lok 'ET' tdcha ;
   'SI' ('NON' tdcha) ; 'MESS' '!!! Calcul recomposition' ; 'FINS' ;
   'SI' lok ; 'MESS' 'OK' ; 'SINO' ; 'ERRE' stop ; 'FINS' ;
'FIN' iikas ;
*
'SI' ('NON' lok) ;
   'MESSAGE' ('CHAINE' 'Il y a eu des erreurs') ;
   'ERREUR' 5 ;
'SINON' ;
   'MESSAGE' ('CHAINE' 'Tout sest bien passe !') ;
'FINSI' ;
*
'SI' interact ;
   'OPTION' 'ECHO' 1 'DONN' 5 ;
'FINSI' ;
*
* End of dgibi file TENS2
*
'FIN' ;
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales