Aster/Cata/Utilitai/veri_matr_tang.py

   1 #@ MODIF veri_matr_tang Utilitai  DATE 18/09/2007   AUTEUR DURAND C.DURAND
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2007  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20 from Cata.cata import *
  21
  22 # ===========================================================================
  23 #              MACRO "VERI_MATR_TANG"
  24 #           -------------------------------------
  25 import Numeric
  26 from Numeric import *
  27 import cPickle,string
  28 import LinearAlgebra as LA
  29 import aster
  30 class TANGENT :
  31
  32      """
  33        VERIFICATION SUR LES MATRICES TANGENTES
  34
  35        Attributs publics :
  36
  37          mat       : matrice tangente
  38          ddl       : nom des degres de liberte
  39          nddl      : nombre de ddl
  40          norme     : norme de la matrice tangente
  41          prec_zero : en-dessous de prec_zero, on ne compare pas les matrices
  42
  43
  44        Methodes publiques
  45
  46          Save       : sauvegarde la matrice dans un fichier
  47          Load       : lit la matrice depuis un fichier
  48          Aster      : lit la matrice depuis l'espace Aster
  49          Matrice    : range la matrice
  50          Difference : comparaison entre la matrice tangente et une autre matrice
  51          Symetrie   : verification du caractere symetrique de la matrice tangente
  52
  53      """
  54
  55
  56      def __init__(self,ddl='',prec_zero=1.E-12) :
  57
  58        """
  59          ddl       : chaine de caracteres designant les ddl (ex: 'UUP')
  60          prec_zero : en-dessous de prec_zero, on ne compare pas les matrices
  61        """
  62
  63        self.ddl = ddl
  64        self.prec_zero = prec_zero
  65
  66
  67      def Load(self,nom_fichier) :
  68
  69        fichier = file(nom_fichier,'r')
  70        self.__dict__ = cPickle.load(fichier)
  71
  72
  73      def Save(self,nom_fichier) :
  74
  75        fichier = file(nom_fichier,'w')
  76        cPickle.dump(self.__dict__,fichier)
  77
  78
  79      def Aster(self,suffixe = 'MATA') :
  80
  81        """
  82          nom : suffixe de l'objet jeveux
  83        """
  84
  85        nom_obj_jeveux = string.ljust('PYTHON.TANGENT.'+suffixe,24)
  86        obj_jeveux = aster.getvectjev(nom_obj_jeveux)
  87        if not obj_jeveux :
  88          raise 'TANGENT : OBJET JEVEUX DE SUFFIXE '+suffixe+' INEXISTANT'
  89        self.Matrice(obj_jeveux)
  90
  91
  92      def Eigen(self) :
  93
  94        self.vp = sort(LA.eigenvalues(self.mat))
  95
  96
  97      def Matrice(self,matrice) :
  98
  99        """
 100          matrice   : la matrice tangente (rangement carre)
 101        """
 102
 103        if type(matrice) == type((1,)) :
 104          matrice = array(list(matrice))
 105        elif type(matrice) == type([]) :
 106          matrice = array(matrice)
 107        matrice = matrice.astype(Float)
 108
 109        nddl = int(len(matrice)**0.5+0.5)
 110        matrice.shape = (nddl,nddl)
 111
 112        self.mat = matrice
 113        self.nddl = nddl
 114
 115        if not self.ddl :
 116          self.ddl = 'D'*nddl
 117        elif len(self.ddl) <> nddl :
 118          raise 'Nommage des DDL incoherents avec la taille de la matrice'
 119
 120        self.norme = trace(dot(transpose(self.mat),self.mat))
 121
 122
 123      def Difference(self,matp,affi_ok=0,prec_diff = 1.E-4) :
 124
 125        """
 126          COMPARAISON RELATIVE DE LA MATRICE TANGENTE AVEC UNE AUTRE MATRICE
 127
 128          matp      : matrice avec laquelle self.mat est comparee
 129          affi_ok   : si oui, on affiche egalement les valeurs qui collent bien
 130          prec_diff : ecart au-dessus duquel on considere que ce n'est pas OK
 131        """
 132
 133        if type(matp) == type((1,)) :
 134          matp = array(list(matp))
 135        elif type(matp) == type([]) :
 136          matp = array(matp)
 137        elif type(matp) == type(self) :
 138          matp = matp.mat
 139        elif type(matp) == type(array([1])) :
 140          pass
 141        else :
 142          raise '1er argument doit etre une matrice (tuple,liste,TANGENT ou Numeric Array)'
 143        matp = ravel(matp)
 144        matp = matp.astype(Float)
 145
 146        if len(matp) <> self.nddl*self.nddl :
 147          raise 'Matrices de tailles differentes'
 148        matp.shape = (self.nddl,self.nddl)
 149
 150        refe = abs(self.mat) + abs(matp)
 151        diff = where(refe > self.prec_zero,abs(self.mat-matp)/(refe+self.prec_zero),0)
 152        nook = nonzero(ravel(diff)  > prec_diff)
 153        ok   = nonzero(ravel(diff) <= prec_diff)
 154
 155        if affi_ok :
 156          affi = [ok,nook]
 157        else :
 158          affi = [nook]
 159
 160        liste_i=[]
 161        liste_j=[]
 162        liste_matt=[]
 163        liste_matp=[]
 164        liste_diff=[]
 165        for ind in affi :
 166 #         print '-'*80
 167          for pos in ind :
 168            i = pos / self.nddl
 169            j = pos % self.nddl
 170 #           print self.ddl[i],self.ddl[j],'  ',(i+1,j+1),'  ',self.mat[i,j],' ',matp[i,j]
 171            liste_i.append(i+1)
 172            liste_j.append(j+1)
 173            liste_matt.append(self.mat[i,j])
 174            liste_matp.append(matp[i,j])
 175            liste_diff.append( abs(self.mat[i,j]-matp[i,j])/ ( abs(self.mat[i,j]) + abs(matp[i,j]) + self.prec_zero))
 176 #       print '-'*80
 177        if self.norme > self.prec_zero :
 178          ecart = (self.mat - matp)/2.
 179          nor_ecart = trace(dot(transpose(ecart),ecart))
 180          nor_diff= nor_ecart / self.norme
 181        else :
 182          nor_diff= 0.
 183        return liste_i,liste_j,liste_matt,liste_matp, liste_diff,nor_diff
 184
 185
 186      def Symetrie(self,prec_diff = 1.E-4) :
 187
 188        """
 189          VERIFICATION QUE LA MATRICE TANGENTE EST SYMETRIQUE
 190
 191          On retourne la norme relative de l'ecart a la symetrie : || (A-At)/2|| / ||A||
 192          On affiche les termes qui s'ecartent de la symetrie
 193
 194          prec_diff : ecart au-dessus duquel on considere que ce n'est pas OK
 195        """
 196
 197        tran = transpose(self.mat)
 198        liste_i,liste_j,liste_matt,liste_matp,liste_diff,nor_diff=self.Difference(tran,affi_ok=0,prec_diff=prec_diff)
 199
 200
 201 #       if self.norme > self.prec_zero :
 202 #         ecart = (self.mat - tran)/2.
 203 #         nor_ecart = trace(dot(transpose(ecart),ecart))
 204 #         return nor_ecart / self.norme
 205 #       else :
 206 #         return 0.
 207        return liste_i,liste_j,liste_matt,liste_matp, liste_diff,nor_diff
 208
 209      def Sauve(self,nom_fichier) :
 210        cPickler.dump(self.__dict__)
 211
 212 def veri_matr_tang_ops(self,SYMETRIE,DIFFERENCE,PRECISION,**args):
 213
 214    """
 215       Ecriture de la macro verif_matrice_tangente_ops
 216    """
 217    import os
 218    from Accas import _F
 219    from Noyau.N_utils import AsType
 220    from Utilitai.Utmess     import UTMESS
 221    from Utilitai.UniteAster import UniteAster
 222
 223    # On importe les definitions des commandes a utiliser dans la macro
 224    CREA_TABLE  = self.get_cmd('CREA_TABLE')
 225
 226   # La macro compte pour 1 dans la numerotation des commandes
 227    self.set_icmd(1)
 228   # Le concept sortant (de type fonction) est nomme ROTGD dans
 229   # le contexte de la macro
 230
 231    self.DeclareOut('TAB_MAT',self.sd)
 232    ier=0
 233    # La macro compte pour 1 dans la numerotation des commandes
 234    self.set_icmd(1)
 235    tgt=TANGENT()
 236    tgt.Aster(suffixe='MATA')
 237    matp=TANGENT()
 238    matp.Aster(suffixe='MATC')
 239    prec_diff = PRECISION
 240    if (SYMETRIE=='OUI') :
 241       list1_i,list1_j,list1_matt,list1_matp,list1_diff,symetgt=tgt.Symetrie(prec_diff)
 242       list2_i,list2_j,list2_matt,list2_matp,list2_diff,symeper=matp.Symetrie(prec_diff)
 243       print 'Symetrie de la matrice tangente',symetgt
 244       print 'Symetrie de la matrice pr pertubation',symeper
 245       aster.affiche('MESSAGE',str(tgt.Difference(matp,prec_diff) ))
 246    if (DIFFERENCE=='OUI'):
 247        liste_i,liste_j,liste_matt,liste_matp,liste_diff,nor_diff=tgt.Difference(matp,prec_diff)
 248        print 'différence entre matrice tangente et matrice par pertubation',nor_diff
 249        TAB_MAT=CREA_TABLE(LISTE=(
 250                      _F(PARA     ='I',LISTE_I = liste_i),
 251                      _F(PARA     ='J',LISTE_I = liste_j),
 252                      _F(PARA     ='MAT_TGTE',LISTE_R = liste_matt),
 253                      _F(PARA     ='MAT_PERT',LISTE_R = liste_matp),
 254                      _F(PARA     ='MAT_DIFF',LISTE_R = liste_diff),
 255                      ))
 256    return
 257
 258
 259 VERI_MATR_TANG=MACRO(nom="VERI_MATR_TANG",op=veri_matr_tang_ops,sd_prod=table_sdaster,
 260                        docu="",reentrant='n',
 261 fr="verification de la matrice tangente : symétrie et différence par rapport a la matrice calculée par perturbation",
 262          regles=(AU_MOINS_UN('SYMETRIE','DIFFERENCE')),
 263          SYMETRIE        =SIMP(statut='f',typ='TXM',defaut="NON",into=("OUI","NON") ),
 264          DIFFERENCE      =SIMP(statut='f',typ='TXM',defaut="OUI",into=("OUI","NON") ),
 265          PRECISION       =SIMP(statut='f',typ='R',defaut=1.E-4 ),
 266 )  ;
 267
 268
 269